From 59b5c763689ed750a51eaba5a35f0ca143de6dfb Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 27 Mar 2023 14:38:18 -0400 Subject: [PATCH 01/17] update_coeffs(L::ADVecProd,...) will recursively call update_coeffs(L.f,...) --- src/differentiation/jaches_products.jl | 4 +++- src/differentiation/vecjac_products.jl | 4 +++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/differentiation/jaches_products.jl b/src/differentiation/jaches_products.jl index 6ce3c4f8..597ed036 100644 --- a/src/differentiation/jaches_products.jl +++ b/src/differentiation/jaches_products.jl @@ -207,10 +207,12 @@ 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) + FwdModeAutoDiffVecProd(f, u, L.vecprod, L.vecprod!, L.cache) end function update_coefficients!(L::FwdModeAutoDiffVecProd, u, p, t) + update_coefficients!(L.f, u, p, t) copy!(L.u, u) L end diff --git a/src/differentiation/vecjac_products.jl b/src/differentiation/vecjac_products.jl index 003e7671..a7921188 100644 --- a/src/differentiation/vecjac_products.jl +++ b/src/differentiation/vecjac_products.jl @@ -65,10 +65,12 @@ 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 From c20400a31cca40a9d737e7b6b3b78f2db86d3289 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 27 Mar 2023 15:31:04 -0400 Subject: [PATCH 02/17] Add tests for recursive updates of f in JacVec etc. --- test/test_jaches_products.jl | 82 ++++++++++++++++++++++++++---------- 1 file changed, 59 insertions(+), 23 deletions(-) diff --git a/test/test_jaches_products.jl b/test/test_jaches_products.jl index f693df2f..26d1c90f 100644 --- a/test/test_jaches_products.jl +++ b/test/test_jaches_products.jl @@ -1,24 +1,62 @@ using SparseDiffTools, ForwardDiff, FiniteDiff, Zygote, IterativeSolvers using LinearAlgebra, Test +import SciMLOperators: update_coefficients, update_coefficients! 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) +end + +# Define state-dependent (i.e. dependent on u/p/t) functions for tests of operators + +mutable struct WrapFunc{F,U,P,T} + func::F + u::U + p::P + t::T +end + +(w::WrapFunc)(u) = sum(w.u) * w.p * w.t * w.func(u) +function (w::WrapFunc)(v, u) + w.func(v, u) + lmul!(sum(w.u) * w.p * w.t, v) +end + +update_coefficients(w::WrapFunc, u, p, t) = WrapFunc(w.func, u, p, t) +function update_coefficients!(w::WrapFunc, u, p, t) + w.u = u + w.p = p + w.t = t end -function h(dy, x) - FiniteDiff.finite_difference_gradient!(dy, g, x) + +# Helper function for testing correct update coefficients behaviour of operators +function update_coefficients_for_test!(L, u, p, t) + update_coefficients!(L, u, p, t) + # Force function hiding inside L to update. Should be a null-op if previous line works correctly + update_coefficients!(L.op.f, u, p, t) end +f = WrapFunc(_f, ones(N) * 2, 1.0, 1.0) +g = WrapFunc(_g, ones(N), 1.0, 1.0) +h = WrapFunc(_h, ones(N), 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))) @@ -67,21 +105,21 @@ 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) @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_for_test!(L, 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 -L = JacVec(f, x, autodiff = AutoFiniteDiff()) +L = JacVec(f, x, 1.0, 1.0; autodiff = AutoFiniteDiff()) @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_for_test!(L, 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 @@ -92,38 +130,36 @@ gmres!(out, L, v) x = rand(N) v = rand(N) -L = HesVec(g, x, autodiff = AutoFiniteDiff()) +L = HesVec(g, x, 1.0, 1.0, autodiff = AutoFiniteDiff()) @test L * x ≈ num_hesvec(g, x, x) @test L * v ≈ num_hesvec(g, x, v) @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_for_test!(L, 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 -L = HesVec(g, x) +L = HesVec(g, x, 1.0, 1.0) @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_for_test!(L, 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 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()) @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_for_test!(L, 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 @@ -134,21 +170,21 @@ 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()) @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) +update_coefficients_for_test!(L, 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 -L = HesVecGrad(h, x) +L = HesVecGrad(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) +update_coefficients_for_test!(L, 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 From 82eee1a4745422996e5235b49356cc54b0d72c8f Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 27 Mar 2023 17:28:29 -0400 Subject: [PATCH 03/17] Update VecJac to new f signature --- src/differentiation/vecjac_products.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/differentiation/vecjac_products.jl b/src/differentiation/vecjac_products.jl index a7921188..7893d890 100644 --- a/src/differentiation/vecjac_products.jl +++ b/src/differentiation/vecjac_products.jl @@ -77,16 +77,16 @@ 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) 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(), @@ -102,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, From 3f73643a1cf91bc5144f1972d5bf86ff3641be14 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 27 Mar 2023 17:30:22 -0400 Subject: [PATCH 04/17] Test VecJac too --- test/runtests.jl | 18 +++++++++--------- test/test_jaches_products.jl | 32 +++----------------------------- test/test_vecjac_products.jl | 20 ++++++++++++-------- 3 files changed, 24 insertions(+), 46 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 09db608c..08c73bfb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,17 +12,17 @@ function activate_gpu_env() end if GROUP == "All" - @time @safetestset "Exact coloring via contraction" begin include("test_contraction.jl") end - @time @safetestset "Greedy distance-1 coloring" begin include("test_greedy_d1.jl") end - @time @safetestset "Greedy star coloring" begin include("test_greedy_star.jl") end - @time @safetestset "Acyclic coloring" begin include("test_acyclic.jl") end - @time @safetestset "Matrix to graph conversion" begin include("test_matrix2graph.jl") end - @time @safetestset "Hessian colorvecs" begin include("test_sparse_hessian.jl") end - @time @safetestset "Integration test" begin include("test_integration.jl") end - @time @safetestset "Special matrices" begin include("test_specialmatrices.jl") end + # @time @safetestset "Exact coloring via contraction" begin include("test_contraction.jl") end + # @time @safetestset "Greedy distance-1 coloring" begin include("test_greedy_d1.jl") end + # @time @safetestset "Greedy star coloring" begin include("test_greedy_star.jl") end + # @time @safetestset "Acyclic coloring" begin include("test_acyclic.jl") end + # @time @safetestset "Matrix to graph conversion" begin include("test_matrix2graph.jl") end + # @time @safetestset "Hessian colorvecs" begin include("test_sparse_hessian.jl") end + # @time @safetestset "Integration test" begin include("test_integration.jl") end + # @time @safetestset "Special matrices" begin include("test_specialmatrices.jl") end @time @safetestset "Jac Vecs and Hes Vecs" begin include("test_jaches_products.jl") end @time @safetestset "Vec Jac Products" begin include("test_vecjac_products.jl") end - @time @safetestset "AD using colorvec vector" begin include("test_ad.jl") end + # @time @safetestset "AD using colorvec vector" begin include("test_ad.jl") end end if GROUP == "GPU" diff --git a/test/test_jaches_products.jl b/test/test_jaches_products.jl index 26d1c90f..e1b826dd 100644 --- a/test/test_jaches_products.jl +++ b/test/test_jaches_products.jl @@ -22,36 +22,10 @@ function _h(dy, x) FiniteDiff.finite_difference_gradient!(dy, _g, x) end -# Define state-dependent (i.e. dependent on u/p/t) functions for tests of operators +# Define state-dependent functions for operator tests -mutable struct WrapFunc{F,U,P,T} - func::F - u::U - p::P - t::T -end - -(w::WrapFunc)(u) = sum(w.u) * w.p * w.t * w.func(u) -function (w::WrapFunc)(v, u) - w.func(v, u) - lmul!(sum(w.u) * w.p * w.t, v) -end - -update_coefficients(w::WrapFunc, u, p, t) = WrapFunc(w.func, u, p, t) -function update_coefficients!(w::WrapFunc, u, p, t) - w.u = u - w.p = p - w.t = t -end - -# Helper function for testing correct update coefficients behaviour of operators -function update_coefficients_for_test!(L, u, p, t) - update_coefficients!(L, u, p, t) - # Force function hiding inside L to update. Should be a null-op if previous line works correctly - update_coefficients!(L.op.f, u, p, t) -end - -f = WrapFunc(_f, ones(N) * 2, 1.0, 1.0) +include("update_coeffs_testutils.jl") +f = WrapFunc(_f, ones(N), 1.0, 1.0) g = WrapFunc(_g, ones(N), 1.0, 1.0) h = WrapFunc(_h, ones(N), 1.0, 1.0) diff --git a/test/test_vecjac_products.jl b/test/test_vecjac_products.jl index a13a8b57..400fd09f 100644 --- a/test/test_vecjac_products.jl +++ b/test/test_vecjac_products.jl @@ -9,16 +9,20 @@ 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, ones(N), 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) +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) +L = VecJac(f, x, 1.0, 1.0; autodiff = AutoFiniteDiff()) +update_coefficients!(L, v, 3.0, 4.0) @test L * v ≈ actual_vjp -# \ No newline at end of file +# From 9036932e1f99e74baa1a7d0a90feb765d2edfc2a Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 27 Mar 2023 21:57:58 -0400 Subject: [PATCH 05/17] Add missing file --- test/update_coeffs_testutils.jl | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 test/update_coeffs_testutils.jl diff --git a/test/update_coeffs_testutils.jl b/test/update_coeffs_testutils.jl new file mode 100644 index 00000000..91399c25 --- /dev/null +++ b/test/update_coeffs_testutils.jl @@ -0,0 +1,28 @@ +# Utilities for testing update coefficient behaviour with state-dependent (i.e. dependent on u/p/t) functions + +mutable struct WrapFunc{F,U,P,T} + func::F + u::U + p::P + t::T +end + +(w::WrapFunc)(u) = sum(w.u) * w.p * w.t * w.func(u) +function (w::WrapFunc)(v, u) + w.func(v, u) + lmul!(sum(w.u) * w.p * w.t, v) +end + +update_coefficients(w::WrapFunc, u, p, t) = WrapFunc(w.func, u, p, t) +function update_coefficients!(w::WrapFunc, u, p, t) + w.u = u + w.p = p + w.t = t +end + +# Helper function for testing correct update coefficients behaviour of operators +function update_coefficients_for_test!(L, u, p, t) + update_coefficients!(L, u, p, t) + # Force function hiding inside L to update. Should be a null-op if previous line works correctly + update_coefficients!(L.op.f, u, p, t) +end \ No newline at end of file From 0b3241abdc259558cb78cee1bc173bd359c2415a Mon Sep 17 00:00:00 2001 From: Vedant Puri Date: Mon, 27 Mar 2023 19:20:23 -0400 Subject: [PATCH 06/17] use new (u, p, t) in autofwdmode --- src/differentiation/jaches_products.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/differentiation/jaches_products.jl b/src/differentiation/jaches_products.jl index 597ed036..6c38e7bd 100644 --- a/src/differentiation/jaches_products.jl +++ b/src/differentiation/jaches_products.jl @@ -218,10 +218,12 @@ function update_coefficients!(L::FwdModeAutoDiffVecProd, u, p, t) end function (L::FwdModeAutoDiffVecProd)(v, p, t) + L = update_coefficients(L, v, p, t) 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 From 530e1b267837e2bdcd9b6a22c7c3e8b256433b98 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 27 Mar 2023 22:13:21 -0400 Subject: [PATCH 07/17] Fix typo --- src/differentiation/jaches_products.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/differentiation/jaches_products.jl b/src/differentiation/jaches_products.jl index 6c38e7bd..5e1e7489 100644 --- a/src/differentiation/jaches_products.jl +++ b/src/differentiation/jaches_products.jl @@ -208,7 +208,7 @@ end function update_coefficients(L::FwdModeAutoDiffVecProd, u, p, t) f = update_coefficients(L.f, u, p, t) - FwdModeAutoDiffVecProd(f, u, L.vecprod, L.vecprod!, L.cache) + FwdModeAutoDiffVecProd(f, u, L.cache, L.vecprod, L.vecprod!) end function update_coefficients!(L::FwdModeAutoDiffVecProd, u, p, t) From 7e4ac735d3922c5d00ca73bdc9fde2526c5fa990 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 27 Mar 2023 22:29:57 -0400 Subject: [PATCH 08/17] Remove u-dependence in WrapFunc --- test/update_coeffs_testutils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/update_coeffs_testutils.jl b/test/update_coeffs_testutils.jl index 91399c25..669ec291 100644 --- a/test/update_coeffs_testutils.jl +++ b/test/update_coeffs_testutils.jl @@ -7,10 +7,10 @@ mutable struct WrapFunc{F,U,P,T} t::T end -(w::WrapFunc)(u) = sum(w.u) * w.p * w.t * w.func(u) +(w::WrapFunc)(u) = w.p * w.t * w.func(u) function (w::WrapFunc)(v, u) w.func(v, u) - lmul!(sum(w.u) * w.p * w.t, v) + lmul!(w.p * w.t, v) end update_coefficients(w::WrapFunc, u, p, t) = WrapFunc(w.func, u, p, t) From c71b3366161b3eb0f0ae1a8b2d49ff0ea6ed3212 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 27 Mar 2023 22:45:03 -0400 Subject: [PATCH 09/17] Make sure to reset state of stateful functions when reconstructing operator --- test/test_jaches_products.jl | 27 +++++++++++++++++++-------- test/update_coeffs_testutils.jl | 6 ++---- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/test/test_jaches_products.jl b/test/test_jaches_products.jl index e1b826dd..cbff19fc 100644 --- a/test/test_jaches_products.jl +++ b/test/test_jaches_products.jl @@ -25,9 +25,9 @@ end # Define state-dependent functions for operator tests include("update_coeffs_testutils.jl") -f = WrapFunc(_f, ones(N), 1.0, 1.0) -g = WrapFunc(_g, ones(N), 1.0, 1.0) -h = WrapFunc(_h, ones(N), 1.0, 1.0) +f = WrapFunc(_f, 1.0, 1.0) +g = WrapFunc(_g, 1.0, 1.0) +h = WrapFunc(_h, 1.0, 1.0) ### @@ -80,6 +80,7 @@ cache4 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(Nothing, eltype(x))), eltype(x) @info "JacVec" 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) @@ -89,6 +90,7 @@ update_coefficients_for_test!(L, v, 3.0, 4.0) dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b) ≈ a*auto_jacvec(f,x,v) + b*_dy 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 @@ -105,8 +107,10 @@ gmres!(out, L, v) x = rand(N) v = rand(N) L = HesVec(g, x, 1.0, 1.0, autodiff = AutoFiniteDiff()) -@test L * x ≈ num_hesvec(g, x, x) -@test L * v ≈ num_hesvec(g, x, v) +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_for_test!(L, v, 3.0, 4.0) @@ -114,6 +118,9 @@ update_coefficients_for_test!(L, v, 3.0, 4.0) dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b) ≈ a*num_hesvec(g,x,v) + b*_dy rtol=1e-2 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 @@ -129,6 +136,7 @@ x = rand(N) v = rand(N) 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 @@ -145,23 +153,26 @@ gmres!(out, L, v) x = rand(N) v = rand(N) 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_for_test!(L, v, 3.0, 4.0) +update_coefficients!(g, x, 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 L = HesVecGrad(h, x, 1.0, 1.0) +update_coefficients!(h, x, 1.0, 1.0) +update_coefficients!(g, 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_for_test!(L, v, 3.0, 4.0) +update_coefficients!(g, x, 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 - -out = similar(v) -gmres!(out, L, v) # diff --git a/test/update_coeffs_testutils.jl b/test/update_coeffs_testutils.jl index 669ec291..9cb75150 100644 --- a/test/update_coeffs_testutils.jl +++ b/test/update_coeffs_testutils.jl @@ -1,8 +1,7 @@ # Utilities for testing update coefficient behaviour with state-dependent (i.e. dependent on u/p/t) functions -mutable struct WrapFunc{F,U,P,T} +mutable struct WrapFunc{F,P,T} func::F - u::U p::P t::T end @@ -13,9 +12,8 @@ function (w::WrapFunc)(v, u) lmul!(w.p * w.t, v) end -update_coefficients(w::WrapFunc, u, p, t) = WrapFunc(w.func, u, p, t) +update_coefficients(w::WrapFunc, u, p, t) = WrapFunc(w.func, p, t) function update_coefficients!(w::WrapFunc, u, p, t) - w.u = u w.p = p w.t = t end From 0c6f6fdc872d63344dbc911f36307cdd83f0cd99 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 27 Mar 2023 23:00:31 -0400 Subject: [PATCH 10/17] Revert commenting out tests --- test/runtests.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 08c73bfb..09db608c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,17 +12,17 @@ function activate_gpu_env() end if GROUP == "All" - # @time @safetestset "Exact coloring via contraction" begin include("test_contraction.jl") end - # @time @safetestset "Greedy distance-1 coloring" begin include("test_greedy_d1.jl") end - # @time @safetestset "Greedy star coloring" begin include("test_greedy_star.jl") end - # @time @safetestset "Acyclic coloring" begin include("test_acyclic.jl") end - # @time @safetestset "Matrix to graph conversion" begin include("test_matrix2graph.jl") end - # @time @safetestset "Hessian colorvecs" begin include("test_sparse_hessian.jl") end - # @time @safetestset "Integration test" begin include("test_integration.jl") end - # @time @safetestset "Special matrices" begin include("test_specialmatrices.jl") end + @time @safetestset "Exact coloring via contraction" begin include("test_contraction.jl") end + @time @safetestset "Greedy distance-1 coloring" begin include("test_greedy_d1.jl") end + @time @safetestset "Greedy star coloring" begin include("test_greedy_star.jl") end + @time @safetestset "Acyclic coloring" begin include("test_acyclic.jl") end + @time @safetestset "Matrix to graph conversion" begin include("test_matrix2graph.jl") end + @time @safetestset "Hessian colorvecs" begin include("test_sparse_hessian.jl") end + @time @safetestset "Integration test" begin include("test_integration.jl") end + @time @safetestset "Special matrices" begin include("test_specialmatrices.jl") end @time @safetestset "Jac Vecs and Hes Vecs" begin include("test_jaches_products.jl") end @time @safetestset "Vec Jac Products" begin include("test_vecjac_products.jl") end - # @time @safetestset "AD using colorvec vector" begin include("test_ad.jl") end + @time @safetestset "AD using colorvec vector" begin include("test_ad.jl") end end if GROUP == "GPU" From b5e2fe38b55d2cb8dffa4c098876bbf6931867c6 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 27 Mar 2023 23:07:25 -0400 Subject: [PATCH 11/17] Update WrapFunc in vecjac tests --- test/test_vecjac_products.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_vecjac_products.jl b/test/test_vecjac_products.jl index 400fd09f..6aea35de 100644 --- a/test/test_vecjac_products.jl +++ b/test/test_vecjac_products.jl @@ -14,7 +14,7 @@ _f(u) = A * u # Define state-dependent functions for operator tests include("update_coeffs_testutils.jl") -f = WrapFunc(_f, ones(N), 1.0, 1.0) +f = WrapFunc(_f, 1.0, 1.0) @info "VecJac" From 6cd7300a3d0c21ad559a572e9dca7910309065d4 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Mon, 27 Mar 2023 23:12:37 -0400 Subject: [PATCH 12/17] Remove update_coefficients_for_test!, test L(v, u, p, t) --- test/test_jaches_products.jl | 44 ++++++++++++++++++++++++--------- test/test_vecjac_products.jl | 8 ++++++ test/update_coeffs_testutils.jl | 11 +++------ 3 files changed, 43 insertions(+), 20 deletions(-) diff --git a/test/test_jaches_products.jl b/test/test_jaches_products.jl index cbff19fc..f89ec418 100644 --- a/test/test_jaches_products.jl +++ b/test/test_jaches_products.jl @@ -1,6 +1,5 @@ using SparseDiffTools, ForwardDiff, FiniteDiff, Zygote, IterativeSolvers using LinearAlgebra, Test -import SciMLOperators: update_coefficients, update_coefficients! using Random Random.seed!(123) @@ -85,9 +84,12 @@ update_coefficients!(f, x, 1.0, 1.0) @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_for_test!(L, v, 3.0, 4.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, 1.0, 1.0; autodiff = AutoFiniteDiff()) update_coefficients!(f, x, 1.0, 1.0) @@ -95,12 +97,15 @@ update_coefficients!(f, x, 1.0, 1.0) @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_for_test!(L, v, 3.0, 4.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" @@ -113,9 +118,12 @@ 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_for_test!(L, v, 3.0, 4.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, 1.0, 1.0) update_coefficients!(g, x, 1.0, 1.0) @@ -125,9 +133,12 @@ num_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_for_test!(L, v, 3.0, 4.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) @@ -141,9 +152,12 @@ update_coefficients!(g, x, 1.0, 1.0) @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_for_test!(L, v, 3.0, 4.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) @@ -159,20 +173,26 @@ update_coefficients!(g, x, 1.0, 1.0) @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_for_test!(L, v, 3.0, 4.0) -update_coefficients!(g, x, 3.0, 4.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, 1.0, 1.0) -update_coefficients!(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_for_test!(L, v, 3.0, 4.0) -update_coefficients!(g, x, 3.0, 4.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) # diff --git a/test/test_vecjac_products.jl b/test/test_vecjac_products.jl index 6aea35de..eb6511e5 100644 --- a/test/test_vecjac_products.jl +++ b/test/test_vecjac_products.jl @@ -20,9 +20,17 @@ f = WrapFunc(_f, 1.0, 1.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 +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 # diff --git a/test/update_coeffs_testutils.jl b/test/update_coeffs_testutils.jl index 9cb75150..8e82c394 100644 --- a/test/update_coeffs_testutils.jl +++ b/test/update_coeffs_testutils.jl @@ -1,4 +1,6 @@ -# Utilities for testing update coefficient behaviour with state-dependent (i.e. dependent on u/p/t) functions +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 @@ -16,11 +18,4 @@ 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 - -# Helper function for testing correct update coefficients behaviour of operators -function update_coefficients_for_test!(L, u, p, t) - update_coefficients!(L, u, p, t) - # Force function hiding inside L to update. Should be a null-op if previous line works correctly - update_coefficients!(L.op.f, u, p, t) end \ No newline at end of file From f3719defdd17e3652f4191459617f8936c1d09da Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Tue, 28 Mar 2023 18:13:17 -0400 Subject: [PATCH 13/17] Remove update_coefficients calls in operator evaluation --- src/differentiation/jaches_products.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/differentiation/jaches_products.jl b/src/differentiation/jaches_products.jl index 5e1e7489..3c6b6c15 100644 --- a/src/differentiation/jaches_products.jl +++ b/src/differentiation/jaches_products.jl @@ -218,12 +218,10 @@ function update_coefficients!(L::FwdModeAutoDiffVecProd, u, p, t) end function (L::FwdModeAutoDiffVecProd)(v, p, t) - L = update_coefficients(L, v, p, t) 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 From 1f8274fbaec7cd543f6617d37334442597ccdf83 Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Tue, 28 Mar 2023 18:40:33 -0400 Subject: [PATCH 14/17] Introduce nonlinearities into test funcs --- test/test_jaches_products.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_jaches_products.jl b/test/test_jaches_products.jl index f89ec418..284b614c 100644 --- a/test/test_jaches_products.jl +++ b/test/test_jaches_products.jl @@ -6,14 +6,14 @@ 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.^2) +_f(x) = A * (x.^2) x = rand(N) v = rand(N) a, b = rand(2) dy = similar(x) -_g(x) = sum(abs2, x) +_g(x) = sum(abs2, x.^2) function _h(x) FiniteDiff.finite_difference_gradient(_g, x) end From c126bf8cb0576b1cde9559b46be393d82833eb0b Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Tue, 28 Mar 2023 18:37:38 -0400 Subject: [PATCH 15/17] Fix aliasing of x and v in tests, as well as unreverted edits to x in in-place hesvecs, and update tests accordingly --- ext/SparseDiffToolsZygote.jl | 2 + src/differentiation/jaches_products.jl | 3 + test/test_jaches_products.jl | 129 +++++++++++++------------ test/test_vecjac_products.jl | 48 +++++++-- 4 files changed, 109 insertions(+), 73 deletions(-) diff --git a/ext/SparseDiffToolsZygote.jl b/ext/SparseDiffToolsZygote.jl index b91e7d31..9275e012 100644 --- a/ext/SparseDiffToolsZygote.jl +++ b/ext/SparseDiffToolsZygote.jl @@ -25,6 +25,7 @@ function SparseDiffTools.numback_hesvec!(dy, f, x, v, cache1 = similar(v), cache g(cache1, x) @. x -= 2ϵ * v g(cache2, x) + @. x += ϵ * v @. dy = (cache1 - cache2) / (2ϵ) end @@ -37,6 +38,7 @@ function SparseDiffTools.numback_hesvec(f, x, v) gxp = g(x) x -= 2ϵ * v gxm = g(x) + @. x += ϵ * v (gxp - gxm) / (2ϵ) end diff --git a/src/differentiation/jaches_products.jl b/src/differentiation/jaches_products.jl index 3c6b6c15..846ba795 100644 --- a/src/differentiation/jaches_products.jl +++ b/src/differentiation/jaches_products.jl @@ -78,6 +78,7 @@ function num_hesvec!(dy, g(cache2, x) @. x -= 2ϵ * v g(cache3, x) + @. x += ϵ * v @. dy = (cache2 - cache3) / (2ϵ) end @@ -110,6 +111,7 @@ function numauto_hesvec!(dy, g(cache1, x) @. x -= 2ϵ * v g(cache2, x) + @. x += ϵ * v @. dy = (cache1 - cache2) / (2ϵ) end @@ -158,6 +160,7 @@ function num_hesvecgrad!(dy, g, x, v, cache2 = similar(v), cache3 = similar(v)) g(cache2, x) @. x -= 2ϵ * v g(cache3, x) + @. x += ϵ * v @. dy = (cache2 - cache3) / (2ϵ) end diff --git a/test/test_jaches_products.jl b/test/test_jaches_products.jl index 284b614c..09f0a4d1 100644 --- a/test/test_jaches_products.jl +++ b/test/test_jaches_products.jl @@ -4,15 +4,23 @@ using LinearAlgebra, Test using Random Random.seed!(123) N = 300 -const A = rand(N, N) - -_f(y, x) = mul!(y, A, x.^2) -_f(x) = A * (x.^2) x = rand(N) v = rand(N) + +# Save original values of x and v to make sure they are not ever mutated +x0 = copy(x) +v0 = copy(v) + a, b = rand(2) dy = similar(x) + +# Define functions for testing + +A = rand(N, N) +_f(y, x) = mul!(y, A, x.^2) +_f(x) = A * (x.^2) + _g(x) = sum(abs2, x.^2) function _h(x) FiniteDiff.finite_difference_gradient(_g, x) @@ -21,7 +29,7 @@ function _h(dy, x) FiniteDiff.finite_difference_gradient!(dy, _g, x) end -# Define state-dependent functions for operator tests +# Make functions state-dependent for operator tests include("update_coeffs_testutils.jl") f = WrapFunc(_f, 1.0, 1.0) @@ -47,38 +55,38 @@ cache2 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(), e similar(v))≈ForwardDiff.hessian(g, x) * v rtol=1e-2 @test num_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-2 -@test numauto_hesvec!(dy, g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-8 +@test numauto_hesvec!(dy, g, x, v)≈ForwardDiff.hessian(g, x) * v @test numauto_hesvec!(dy, g, x, v, ForwardDiff.GradientConfig(g, x), similar(v), - similar(v))≈ForwardDiff.hessian(g, x) * v rtol=1e-8 -@test numauto_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-8 + similar(v))≈ForwardDiff.hessian(g, x) * v +@test numauto_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v -@test autonum_hesvec!(dy, g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-2 +@test autonum_hesvec!(dy, g, x, v)≈ForwardDiff.hessian(g, x) * v @test autonum_hesvec!(dy, g, x, v, cache1, cache2)≈ForwardDiff.hessian(g, x) * v rtol=1e-2 -@test autonum_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-8 +@test autonum_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v -@test numback_hesvec!(dy, g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-8 -@test numback_hesvec!(dy, g, x, v, similar(v), similar(v))≈ForwardDiff.hessian(g, x) * v rtol=1e-8 -@test numback_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-8 +@test numback_hesvec!(dy, g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-6 +@test numback_hesvec!(dy, g, x, v, similar(v), similar(v))≈ForwardDiff.hessian(g, x) * v rtol=1e-6 +@test numback_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-6 cache3 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(Nothing, eltype(x))), eltype(x), 1 }.(x, ForwardDiff.Partials.(tuple.(v))) cache4 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(Nothing, eltype(x))), eltype(x), 1 }.(x, ForwardDiff.Partials.(tuple.(v))) -@test autoback_hesvec!(dy, g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-8 -@test autoback_hesvec!(dy, g, x, v, cache3, cache4)≈ForwardDiff.hessian(g, x) * v rtol=1e-8 -@test autoback_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-8 +@test autoback_hesvec!(dy, g, x, v)≈ForwardDiff.hessian(g, x) * v +@test autoback_hesvec!(dy, g, x, v, cache3, cache4)≈ForwardDiff.hessian(g, x) * v +@test autoback_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v @test num_hesvecgrad!(dy, h, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-2 @test num_hesvecgrad!(dy, h, x, v, similar(v), similar(v))≈ForwardDiff.hessian(g, x) * v rtol=1e-2 @test num_hesvecgrad(h, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-2 -@test auto_hesvecgrad!(dy, h, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-2 -@test auto_hesvecgrad!(dy, h, x, v, cache1, cache2)≈ForwardDiff.hessian(g, x) * v rtol=1e-2 -@test auto_hesvecgrad(h, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-2 +@test auto_hesvecgrad!(dy, h, x, v)≈ForwardDiff.hessian(g, x) * v +@test auto_hesvecgrad!(dy, h, x, v, cache1, cache2)≈ForwardDiff.hessian(g, x) * v +@test auto_hesvecgrad(h, x, v)≈ForwardDiff.hessian(g, x) * v @info "JacVec" -L = JacVec(f, x, 1.0, 1.0) +L = JacVec(f, copy(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) @@ -86,12 +94,12 @@ update_coefficients!(f, x, 1.0, 1.0) 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, 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 +@test mul!(dy, L, x) ≈ auto_jacvec(f, v, x) +dy=rand(N);_dy=copy(dy);@test mul!(dy,L,x,a,b) ≈ a*auto_jacvec(f,v,x) + b*_dy update_coefficients!(f, v, 5.0, 6.0) -@test L(dy, v, 5.0, 6.0) ≈ auto_jacvec(f,x,v) +@test L(dy, v, 5.0, 6.0) ≈ auto_jacvec(f,v,v) -L = JacVec(f, x, 1.0, 1.0; autodiff = AutoFiniteDiff()) +L = JacVec(f, copy(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) @@ -99,74 +107,63 @@ update_coefficients!(f, x, 1.0, 1.0) 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, 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 +@test mul!(dy, L, x)≈num_jacvec(f, v, x) rtol=1e-6 +dy=rand(N);_dy=copy(dy);@test mul!(dy,L,x,a,b) ≈ a*num_jacvec(f,v,x) + 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 +@test L(dy, v, 5.0, 6.0) ≈ num_jacvec(f,v,v) rtol=1e-6 out = similar(v) @test_nowarn gmres!(out, L, v) @info "HesVec" -x = rand(N) -v = rand(N) -L = HesVec(g, x, 1.0, 1.0, autodiff = AutoFiniteDiff()) +L = HesVec(g, copy(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, 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 +@test mul!(dy, L, x)≈num_hesvec(g, v, x) rtol=1e-2 +dy=rand(N);_dy=copy(dy);@test mul!(dy,L,x,a,b) ≈ a*num_hesvec(g,v,x) + 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 +@test L(dy, v, 5.0, 6.0) ≈ num_hesvec(g,v,v) rtol=1e-2 -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) +L = HesVec(g, copy(x), 1.0, 1.0) @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 +@test mul!(dy, L, v)≈numauto_hesvec(g, x, v) +dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,0)≈a*numauto_hesvec(g,x,v)+0*_dy 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 +@test mul!(dy, L, x)≈numauto_hesvec(g, v, x) +dy=rand(N);_dy=copy(dy);@test mul!(dy,L,x,a,b)≈a*numauto_hesvec(g,v,x)+b*_dy update_coefficients!(g, v, 5.0, 6.0) -@test L(dy, v, 5.0, 6.0) ≈ numauto_hesvec(g,x,v) rtol=1e-2 +@test L(dy, v, 5.0, 6.0) ≈ numauto_hesvec(g,v,v) out = similar(v) gmres!(out, L, v) -x = rand(N) -v = rand(N) - -L = HesVec(g, x, 1.0, 1.0; autodiff = AutoZygote()) +L = HesVec(g, copy(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 +@test mul!(dy, L, v)≈autoback_hesvec(g, x, v) +dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*autoback_hesvec(g,x,v)+b*_dy 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 +@test mul!(dy, L, x)≈autoback_hesvec(g, v, x) +dy=rand(N);_dy=copy(dy);@test mul!(dy,L,x,a,b)≈a*autoback_hesvec(g,v,x)+b*_dy update_coefficients!(g, v, 5.0, 6.0) -@test L(dy, v, 5.0, 6.0) ≈ autoback_hesvec(g,x,v) rtol=1e-2 +@test L(dy, v, 5.0, 6.0) ≈ autoback_hesvec(g,v,v) out = similar(v) gmres!(out, L, v) @info "HesVecGrad" -x = rand(N) -v = rand(N) -L = HesVecGrad(h, x, 1.0, 1.0; autodiff = AutoFiniteDiff()) +L = HesVecGrad(h, copy(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 @@ -174,25 +171,31 @@ update_coefficients!(g, x, 1.0, 1.0) @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 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 +@test mul!(dy, L, x)≈num_hesvec(g, v, x) rtol=1e-2 +dy=rand(N);_dy=copy(dy);@test mul!(dy,L,x,a,b)≈a*num_hesvec(g,v,x)+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 +@test L(dy, v, 5.0, 6.0) ≈ num_hesvec(g,v,v) rtol=1e-2 -L = HesVecGrad(h, x, 1.0, 1.0) +L = HesVecGrad(h, copy(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 +dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*numauto_hesvec(g,x,v)+b*_dy 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 +@test mul!(dy, L, x)≈numauto_hesvec(g, v, x) +dy=rand(N);_dy=copy(dy);@test mul!(dy,L,x,a,b)≈a*numauto_hesvec(g,v,x)+b*_dy 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 +@test L(dy, v, 5.0, 6.0) ≈ numauto_hesvec(g,v,v) out = similar(v) gmres!(out, L, v) + +# Test that x and v were not mutated +# x's rtol can't be too large since it is mutated and then restored in some algorithms +@test x ≈ x0 +@test v ≈ v0 + # diff --git a/test/test_vecjac_products.jl b/test/test_vecjac_products.jl index eb6511e5..2978d8fc 100644 --- a/test/test_vecjac_products.jl +++ b/test/test_vecjac_products.jl @@ -4,33 +4,61 @@ using LinearAlgebra, Test using Random Random.seed!(123) N = 300 -const A = rand(N, N) +# Use Float32 since Zygote defaults to Float32 x = rand(Float32, N) v = rand(Float32, N) +# Save original values of x and v to make sure they are not ever mutated +x0 = copy(x) +v0 = copy(v) + +a, b = rand(2) +dy = similar(x) + +A = rand(Float32, N, N) _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) +f = WrapFunc(_f, 1.0f0, 1.0f0) + +# Compute Jacobian via Zygote @info "VecJac" -L = VecJac(f, x, 1.0, 1.0) +L = VecJac(f, copy(x), 1.0f0, 1.0f0; autodiff = AutoZygote()) +update_coefficients!(f, x, 1.0, 1.0) +actual_jac = Zygote.jacobian(f, x)[1] +@test L * x ≈ actual_jac' * x +@test L * v ≈ actual_jac' * v +@test mul!(dy, L, v) ≈ actual_jac' * v 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 +actual_jac = Zygote.jacobian(f, v)[1] +@test mul!(dy, L, x) ≈ actual_jac' * x +_dy=copy(dy); @test mul!(dy,L,x,a,b) ≈ a*actual_jac'*x + b*_dy 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 +actual_jac = Zygote.jacobian(f, v)[1] +@test L(dy, v, 5.0, 6.0) ≈ actual_jac' * v -L = VecJac(f, x, 1.0, 1.0; autodiff = AutoFiniteDiff()) +L = VecJac(f, copy(x), 1.0f0, 1.0f0; autodiff = AutoFiniteDiff()) +update_coefficients!(f, x, 1.0, 1.0) +actual_jac = Zygote.jacobian(f, x)[1] +@test L * x ≈ actual_jac' * x +@test L * v ≈ actual_jac' * v +@test mul!(dy, L, v) ≈ actual_jac' * v update_coefficients!(L, v, 3.0, 4.0) update_coefficients!(f, v, 3.0, 4.0) -@test L * v ≈ actual_vjp +actual_jac = Zygote.jacobian(f, v)[1] +@test mul!(dy, L, x) ≈ actual_jac' * x +_dy=copy(dy); @test mul!(dy,L,x,a,b) ≈ a*actual_jac'*x + b*_dy update_coefficients!(f, v, 5.0, 6.0) -@test L(copy(v), v, 5.0, 6.0) ≈ actual_vjp2 +actual_jac = Zygote.jacobian(f, v)[1] +@test L(dy, v, 5.0, 6.0) ≈ actual_jac' * v + +# Test that x and v were not mutated +@test x ≈ x0 +@test v ≈ v0 # From c8b45c07a1cf8165d000a7fdc92fa40189b5c3ef Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Tue, 28 Mar 2023 20:15:41 -0400 Subject: [PATCH 16/17] Tighten up a few more tolerances --- test/test_jaches_products.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_jaches_products.jl b/test/test_jaches_products.jl index 09f0a4d1..200b2584 100644 --- a/test/test_jaches_products.jl +++ b/test/test_jaches_products.jl @@ -64,9 +64,9 @@ cache2 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(), e @test autonum_hesvec!(dy, g, x, v, cache1, cache2)≈ForwardDiff.hessian(g, x) * v rtol=1e-2 @test autonum_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v -@test numback_hesvec!(dy, g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-6 -@test numback_hesvec!(dy, g, x, v, similar(v), similar(v))≈ForwardDiff.hessian(g, x) * v rtol=1e-6 -@test numback_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v rtol=1e-6 +@test numback_hesvec!(dy, g, x, v)≈ForwardDiff.hessian(g, x) * v +@test numback_hesvec!(dy, g, x, v, similar(v), similar(v))≈ForwardDiff.hessian(g, x) * v +@test numback_hesvec(g, x, v)≈ForwardDiff.hessian(g, x) * v cache3 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(Nothing, eltype(x))), eltype(x), 1 }.(x, ForwardDiff.Partials.(tuple.(v))) @@ -181,7 +181,7 @@ 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 +@test mul!(dy, L, v)≈numauto_hesvec(g, x, v) dy=rand(N);_dy=copy(dy);@test mul!(dy,L,v,a,b)≈a*numauto_hesvec(g,x,v)+b*_dy for op in (L, g, h) update_coefficients!(op, v, 3.0, 4.0) end @test mul!(dy, L, x)≈numauto_hesvec(g, v, x) From 57f4a5568b39e68bc6e0baffa1b2c09a44add59c Mon Sep 17 00:00:00 2001 From: Gaurav Arya Date: Thu, 30 Mar 2023 17:08:58 -0400 Subject: [PATCH 17/17] Revert extraneous change to out-of-place zygote hessian --- ext/SparseDiffToolsZygote.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/ext/SparseDiffToolsZygote.jl b/ext/SparseDiffToolsZygote.jl index 9275e012..eb6df440 100644 --- a/ext/SparseDiffToolsZygote.jl +++ b/ext/SparseDiffToolsZygote.jl @@ -38,7 +38,6 @@ function SparseDiffTools.numback_hesvec(f, x, v) gxp = g(x) x -= 2ϵ * v gxm = g(x) - @. x += ϵ * v (gxp - gxm) / (2ϵ) end