Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

have update_coeffs(L::ADVecProd,) recursively update L.f #232

Merged
merged 17 commits into from
Apr 8, 2023
Merged

have update_coeffs(L::ADVecProd,) recursively update L.f #232

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
59b5c76
update_coeffs(L::ADVecProd,...) will recursively call update_coeffs(L…
vpuri3 Mar 27, 2023
c20400a
Add tests for recursive updates of f in JacVec etc.
gaurav-arya Mar 27, 2023
82eee1a
Update VecJac to new f signature
gaurav-arya Mar 27, 2023
3f73643
Test VecJac too
gaurav-arya Mar 27, 2023
9036932
Add missing file
gaurav-arya Mar 28, 2023
0b3241a
use new (u, p, t) in autofwdmode
vpuri3 Mar 27, 2023
530e1b2
Fix typo
gaurav-arya Mar 28, 2023
7e4ac73
Remove u-dependence in WrapFunc
gaurav-arya Mar 28, 2023
c71b336
Make sure to reset state of stateful functions when reconstructing op…
gaurav-arya Mar 28, 2023
0c6f6fd
Revert commenting out tests
gaurav-arya Mar 28, 2023
b5e2fe3
Update WrapFunc in vecjac tests
gaurav-arya Mar 28, 2023
6cd7300
Remove update_coefficients_for_test!, test L(v, u, p, t)
gaurav-arya Mar 28, 2023
f3719de
Remove update_coefficients calls in operator evaluation
gaurav-arya Mar 28, 2023
1f8274f
Introduce nonlinearities into test funcs
gaurav-arya Mar 28, 2023
c126bf8
Fix aliasing of x and v in tests, as well as unreverted edits to x in…
gaurav-arya Mar 28, 2023
c8b45c0
Tighten up a few more tolerances
gaurav-arya Mar 29, 2023
57f4a55
Revert extraneous change to out-of-place zygote hessian
gaurav-arya Mar 30, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 5 additions & 1 deletion src/differentiation/jaches_products.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -207,19 +207,23 @@ struct FwdModeAutoDiffVecProd{F,U,C,V,V!} <: AbstractAutoDiffVecProd
end

function update_coefficients(L::FwdModeAutoDiffVecProd, u, p, t)
FwdModeAutoDiffVecProd(L.f, u, L.vecprod, L.vecprod!, L.cache)
f = update_coefficients(L.f, u, p, t)
vpuri3 marked this conversation as resolved.
Show resolved Hide resolved
FwdModeAutoDiffVecProd(f, u, L.cache, L.vecprod, L.vecprod!)
end

function update_coefficients!(L::FwdModeAutoDiffVecProd, u, p, t)
update_coefficients!(L.f, u, p, t)
copy!(L.u, u)
L
end

function (L::FwdModeAutoDiffVecProd)(v, p, t)
L = update_coefficients(L, v, p, t)
vpuri3 marked this conversation as resolved.
Show resolved Hide resolved
L.vecprod(L.f, L.u, v)
end

function (L::FwdModeAutoDiffVecProd)(dv, v, p, t)
update_coefficients!(L, v, p, t)
L.vecprod!(dv, L.f, L.u, v, L.cache...)
end

Expand Down
16 changes: 9 additions & 7 deletions src/differentiation/vecjac_products.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,26 +65,28 @@ struct RevModeAutoDiffVecProd{ad,iip,oop,F,U,C,V,V!} <: AbstractAutoDiffVecProd
end

function update_coefficients(L::RevModeAutoDiffVecProd, u, p, t)
RevModeAutoDiffVecProd(L.f, u, L.vecprod, L.vecprod!, L.cache)
f = update_coefficients(L.f, u, p, t)
RevModeAutoDiffVecProd(f, u, L.vecprod, L.vecprod!, L.cache)
end

function update_coefficients!(L::RevModeAutoDiffVecProd, u, p, t)
update_coefficients!(L.f, u, p, t)
copy!(L.u, u)
L
end

# Interpret the call as df/du' * u
function (L::RevModeAutoDiffVecProd)(v, p, t)
L.vecprod(_u -> L.f(_u, p, t), L.u, v)
L.vecprod(L.f, L.u, v)
vpuri3 marked this conversation as resolved.
Show resolved Hide resolved
end

# prefer non in-place method
function (L::RevModeAutoDiffVecProd{ad,iip,true})(dv, v, p, t) where{ad,iip}
L.vecprod!(dv, _u -> L.f(_u, p, t), L.u, v, L.cache...)
L.vecprod!(dv, L.f, L.u, v, L.cache...)
end

function (L::RevModeAutoDiffVecProd{ad,true,false})(dv, v, p, t) where{ad}
L.vecprod!(dv, (_du, _u) -> L.f(_du, _u, p, t), L.u, v, L.cache...)
L.vecprod!(dv, L.f, L.u, v, L.cache...)
end

function VecJac(f, u::AbstractArray, p = nothing, t = nothing; autodiff = AutoFiniteDiff(),
Expand All @@ -100,11 +102,11 @@ function VecJac(f, u::AbstractArray, p = nothing, t = nothing; autodiff = AutoFi

cache = (similar(u), similar(u),)

outofplace = static_hasmethod(f, typeof((u, p, t)))
isinplace = static_hasmethod(f, typeof((u, u, p, t)))
outofplace = static_hasmethod(f, typeof((u,)))
isinplace = static_hasmethod(f, typeof((u, u,)))

if !(isinplace) & !(outofplace)
error("$f must have signature f(u, p, t), or f(du, u, p, t)")
error("$f must have signature f(u), or f(du, u)")
end

L = RevModeAutoDiffVecProd(f, u, cache, vecprod, vecprod!; autodiff = autodiff,
Expand Down
93 changes: 67 additions & 26 deletions test/test_jaches_products.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,20 +5,31 @@ using Random
Random.seed!(123)
N = 300
const A = rand(N, N)
f(y, x) = mul!(y, A, x)
f(x) = A * x

_f(y, x) = mul!(y, A, x)
_f(x) = A * x

x = rand(N)
v = rand(N)
a, b = rand(2)
dy = similar(x)
g(x) = sum(abs2, x)
function h(x)
FiniteDiff.finite_difference_gradient(g, x)
_g(x) = sum(abs2, x)
function _h(x)
FiniteDiff.finite_difference_gradient(_g, x)
end
function h(dy, x)
FiniteDiff.finite_difference_gradient!(dy, g, x)
function _h(dy, x)
FiniteDiff.finite_difference_gradient!(dy, _g, x)
end

# Define state-dependent functions for operator tests

include("update_coeffs_testutils.jl")
vpuri3 marked this conversation as resolved.
Show resolved Hide resolved
f = WrapFunc(_f, 1.0, 1.0)
g = WrapFunc(_g, 1.0, 1.0)
h = WrapFunc(_h, 1.0, 1.0)

###

cache1 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(), eltype(x))),
eltype(x), 1}.(x, ForwardDiff.Partials.(tuple.(v)))
cache2 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(), eltype(x))), eltype(x), 1}.(x, ForwardDiff.Partials.(tuple.(v)))
Expand Down Expand Up @@ -67,65 +78,86 @@ cache4 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(Nothing, eltype(x))), eltype(x)

@info "JacVec"

L = JacVec(f, x)
L = JacVec(f, x, 1.0, 1.0)
update_coefficients!(f, x, 1.0, 1.0)
@test L * x ≈ auto_jacvec(f, x, x)
@test L * v ≈ auto_jacvec(f, x, v)
@test mul!(dy, L, v) ≈ auto_jacvec(f, x, v)
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b) ≈ a*auto_jacvec(f,x,v) + b*_dy
update_coefficients!(L, v, nothing, 0.0)
update_coefficients!(L, v, 3.0, 4.0)
update_coefficients!(f, v, 3.0, 4.0)
@test mul!(dy, L, v) ≈ auto_jacvec(f, v, v)
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b) ≈ a*auto_jacvec(f,x,v) + b*_dy
update_coefficients!(f, v, 5.0, 6.0)
@test L(dy, v, 5.0, 6.0) ≈ auto_jacvec(f,x,v)

L = JacVec(f, x, autodiff = AutoFiniteDiff())
L = JacVec(f, x, 1.0, 1.0; autodiff = AutoFiniteDiff())
update_coefficients!(f, x, 1.0, 1.0)
@test L * x ≈ num_jacvec(f, x, x)
@test L * v ≈ num_jacvec(f, x, v)
@test mul!(dy, L, v)≈num_jacvec(f, x, v) rtol=1e-6
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b) ≈ a*num_jacvec(f,x,v) + b*_dy rtol=1e-6
update_coefficients!(L, v, nothing, 0.0)
update_coefficients!(L, v, 3.0, 4.0)
update_coefficients!(f, v, 3.0, 4.0)
@test mul!(dy, L, v)≈num_jacvec(f, v, v) rtol=1e-6
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b) ≈ a*num_jacvec(f,x,v) + b*_dy rtol=1e-6
update_coefficients!(f, v, 5.0, 6.0)
@test L(dy, v, 5.0, 6.0) ≈ num_jacvec(f,x,v) rtol=1e-6

out = similar(v)
gmres!(out, L, v)
@test_nowarn gmres!(out, L, v)

@info "HesVec"

x = rand(N)
v = rand(N)
L = HesVec(g, x, autodiff = AutoFiniteDiff())
@test L * x ≈ num_hesvec(g, x, x)
@test L * v ≈ num_hesvec(g, x, v)
L = HesVec(g, x, 1.0, 1.0, autodiff = AutoFiniteDiff())
update_coefficients!(g, x, 1.0, 1.0)
@test L * x ≈ num_hesvec(g, x, x) rtol=1e-2
num_hesvec(g, x, x)
@test L * v ≈ num_hesvec(g, x, v) rtol=1e-2
@test mul!(dy, L, v)≈num_hesvec(g, x, v) rtol=1e-2
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b) ≈ a*num_hesvec(g,x,v) + b*_dy rtol=1e-2
update_coefficients!(L, v, nothing, 0.0)
update_coefficients!(L, v, 3.0, 4.0)
update_coefficients!(g, v, 3.0, 4.0)
@test mul!(dy, L, v)≈num_hesvec(g, v, v) rtol=1e-2
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b) ≈ a*num_hesvec(g,x,v) + b*_dy rtol=1e-2
update_coefficients!(g, v, 5.0, 6.0)
@test L(dy, v, 5.0, 6.0) ≈ num_hesvec(g,x,v) rtol=1e-2

L = HesVec(g, x)
L = HesVec(g, x, 1.0, 1.0)
update_coefficients!(g, x, 1.0, 1.0)
numauto_hesvec(g, x, x)
num_hesvec(g, x, x)
@test L * x ≈ numauto_hesvec(g, x, x)
@test L * v ≈ numauto_hesvec(g, x, v)
@test mul!(dy, L, v)≈numauto_hesvec(g, x, v) rtol=1e-8
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*numauto_hesvec(g,x,v)+b*_dy rtol=1e-8
update_coefficients!(L, v, nothing, 0.0)
update_coefficients!(L, v, 3.0, 4.0)
update_coefficients!(g, v, 3.0, 4.0)
@test mul!(dy, L, v)≈numauto_hesvec(g, v, v) rtol=1e-8
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*numauto_hesvec(g,x,v)+b*_dy rtol=1e-8
update_coefficients!(g, v, 5.0, 6.0)
@test L(dy, v, 5.0, 6.0) ≈ numauto_hesvec(g,x,v) rtol=1e-2

out = similar(v)
gmres!(out, L, v)

using Zygote

x = rand(N)
v = rand(N)

L = HesVec(g, x, autodiff = AutoZygote())
L = HesVec(g, x, 1.0, 1.0; autodiff = AutoZygote())
update_coefficients!(g, x, 1.0, 1.0)
@test L * x ≈ autoback_hesvec(g, x, x)
@test L * v ≈ autoback_hesvec(g, x, v)
@test mul!(dy, L, v)≈autoback_hesvec(g, x, v) rtol=1e-8
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*autoback_hesvec(g,x,v)+b*_dy rtol=1e-8
update_coefficients!(L, v, nothing, 0.0)
update_coefficients!(L, v, 3.0, 4.0)
update_coefficients!(g, v, 3.0, 4.0)
@test mul!(dy, L, v)≈autoback_hesvec(g, v, v) rtol=1e-8
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*autoback_hesvec(g,x,v)+b*_dy rtol=1e-8
update_coefficients!(g, v, 5.0, 6.0)
@test L(dy, v, 5.0, 6.0) ≈ autoback_hesvec(g,x,v) rtol=1e-2

out = similar(v)
gmres!(out, L, v)
Expand All @@ -134,23 +166,32 @@ gmres!(out, L, v)

x = rand(N)
v = rand(N)
L = HesVecGrad(h, x, autodiff = AutoFiniteDiff())
L = HesVecGrad(h, x, 1.0, 1.0; autodiff = AutoFiniteDiff())
update_coefficients!(h, x, 1.0, 1.0)
update_coefficients!(g, x, 1.0, 1.0)
@test L * x ≈ num_hesvec(g, x, x) rtol=1e-2
@test L * v ≈ num_hesvec(g, x, v) rtol=1e-2
@test mul!(dy, L, v)≈num_hesvec(g, x, v) rtol=1e-2
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*num_hesvec(g,x,v)+b*_dy rtol=1e-2
update_coefficients!(L, v, nothing, 0.0)
for op in (L, g, h) update_coefficients!(op, v, 3.0, 4.0) end
@test mul!(dy, L, v)≈num_hesvec(g, v, v) rtol=1e-2
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*num_hesvec(g,x,v)+b*_dy rtol=1e-2
update_coefficients!(g, v, 5.0, 6.0)
@test L(dy, v, 5.0, 6.0) ≈ num_hesvec(g,x,v) rtol=1e-2

L = HesVecGrad(h, x)
L = HesVecGrad(h, x, 1.0, 1.0)
update_coefficients!(g, x, 1.0, 1.0)
update_coefficients!(h, x, 1.0, 1.0)
@test L * x ≈ autonum_hesvec(g, x, x)
@test L * v ≈ numauto_hesvec(g, x, v)
@test mul!(dy, L, v)≈numauto_hesvec(g, x, v) rtol=1e-8
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*numauto_hesvec(g,x,v)+b*_dy rtol=1e-8
update_coefficients!(L, v, nothing, 0.0)
for op in (L, g, h) update_coefficients!(op, v, 3.0, 4.0) end
@test mul!(dy, L, v)≈numauto_hesvec(g, v, v) rtol=1e-8
dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*numauto_hesvec(g,x,v)+b*_dy rtol=1e-8
update_coefficients!(g, v, 5.0, 6.0)
update_coefficients!(h, v, 5.0, 6.0)
@test L(dy, v, 5.0, 6.0) ≈ num_hesvec(g,x,v) rtol=1e-2

out = similar(v)
gmres!(out, L, v)
Expand Down
28 changes: 20 additions & 8 deletions test/test_vecjac_products.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,16 +9,28 @@ const A = rand(N, N)
x = rand(Float32, N)
v = rand(Float32, N)

f(du,u,p,t) = mul!(du, A, u)
f(u,p,t) = A * u
_f(du,u) = mul!(du, A, u)
_f(u) = A * u

# Define state-dependent functions for operator tests
include("update_coeffs_testutils.jl")
f = WrapFunc(_f, 1.0, 1.0)

@info "VecJac"

L = VecJac(f, x)
actual_vjp = Zygote.jacobian(x -> f(x, nothing, 0.0), x)[1]' * v
update_coefficients!(L, v, nothing, 0.0)
L = VecJac(f, x, 1.0, 1.0)
update_coefficients!(L, v, 3.0, 4.0)
update_coefficients!(f, v, 3.0, 4.0)
actual_vjp = Zygote.jacobian(f, x)[1]' * v
@test L * v ≈ actual_vjp
L = VecJac(f, x; autodiff = AutoFiniteDiff())
update_coefficients!(L, v, nothing, 0.0)
update_coefficients!(f, v, 5.0, 6.0)
actual_vjp2 = Zygote.jacobian(f, x)[1]' * v
@test L(copy(v), v, 5.0, 6.0) ≈ actual_vjp2

L = VecJac(f, x, 1.0, 1.0; autodiff = AutoFiniteDiff())
update_coefficients!(L, v, 3.0, 4.0)
update_coefficients!(f, v, 3.0, 4.0)
@test L * v ≈ actual_vjp
#
update_coefficients!(f, v, 5.0, 6.0)
@test L(copy(v), v, 5.0, 6.0) ≈ actual_vjp2
#
21 changes: 21 additions & 0 deletions test/update_coeffs_testutils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
import SciMLOperators: update_coefficients, update_coefficients!

# Wrapper function for testing update coefficient behaviour with state-dependent (i.e. dependent on u/p/t) functions

mutable struct WrapFunc{F,P,T}
func::F
p::P
t::T
end

(w::WrapFunc)(u) = w.p * w.t * w.func(u)
function (w::WrapFunc)(v, u)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@vpuri3 can you take a look at how this operator is using u as part of its state, and applying it to a different RHS v? Is this something that is invalid in your opinion? Because your recent changes don't work with this -- this also relates to the discussion in SciML/SciMLOperators.jl#57.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(This is a separate, SciMLOperators issue, so resolved for our purposes here)

w.func(v, u)
lmul!(w.p * w.t, v)
end

update_coefficients(w::WrapFunc, u, p, t) = WrapFunc(w.func, p, t)
function update_coefficients!(w::WrapFunc, u, p, t)
w.p = p
w.t = t
end