diff --git a/src/jacobians.jl b/src/jacobians.jl index e1b742c..108a60d 100644 --- a/src/jacobians.jl +++ b/src/jacobians.jl @@ -287,6 +287,13 @@ function finite_difference_jacobian( dir = true) where {T1, T2, T3, T4, cType, sType, fdtype, returntype} x1, fx, fx1 = cache.x1, cache.fx, cache.fx1 + # Issue #213: cache.x1 may have been initialized via `similar(x)` (e.g. by + # DifferentiationInterface) or built around a previous x. Synchronize it + # with the current x so the cache is observably consistent after the call. + if x1 isa AbstractArray && ArrayInterface.ismutable(x1) && x1 !== x + copyto!(x1, x) + end + if !(f_in isa Nothing) vecfx = _vec(f_in) elseif fdtype == Val(:forward) @@ -297,7 +304,6 @@ function finite_difference_jacobian( vecfx = _vec(fx) end vecx = _vec(x) - vecx1 = _vec(x1) J = jac_prototype isa Nothing ? (sparsity isa Nothing ? Array{eltype(x), 2}(undef, length(vecfx), 0) : zeros(eltype(x), size(sparsity))) : zero(jac_prototype) @@ -343,11 +349,14 @@ function finite_difference_jacobian( end end elseif fdtype == Val(:central) + # Both halves of the central difference must perturb around the *current* + # x. Reading the unperturbed components from `cache.x1` (issue #213) is + # unsafe — the cache may have been built via `similar(x)` or reused at a + # different x — so we always perturb around `vecx` directly. function calculate_Ji_central(i) - x1_save = ArrayInterface.allowed_getindex(vecx1, i) x_save = ArrayInterface.allowed_getindex(vecx, i) - epsilon = compute_epsilon(Val(:forward), x1_save, relstep, absstep, dir) - _vecx1 = setindex(vecx1, x1_save+epsilon, i) + epsilon = compute_epsilon(Val(:forward), x_save, relstep, absstep, dir) + _vecx1 = setindex(vecx, x_save+epsilon, i) _vecx = setindex(vecx, x_save-epsilon, i) _x1 = reshape(_vecx1, axes(x)) _x = reshape(_vecx, axes(x)) @@ -366,10 +375,10 @@ function finite_difference_jacobian( dx = calculate_Ji_central(color_i) J = J + _make_Ji(J, eltype(x), dx, color_i, nrows, ncols) else - tmp = norm(vecx1 .* (colorvec .== color_i)) + tmp = norm(vecx .* (colorvec .== color_i)) epsilon = compute_epsilon( Val(:forward), sqrt(tmp), relstep, absstep, dir) - _vecx1 = @. vecx1 + epsilon * (colorvec == color_i) + _vecx1 = @. vecx + epsilon * (colorvec == color_i) _vecx = @. vecx - epsilon * (colorvec == color_i) _x1 = reshape(_vecx1, axes(x)) _x = reshape(_vecx, axes(x)) diff --git a/test/cache_reuse_tests.jl b/test/cache_reuse_tests.jl new file mode 100644 index 0000000..53a0807 --- /dev/null +++ b/test/cache_reuse_tests.jl @@ -0,0 +1,162 @@ +using FiniteDiff, LinearAlgebra, SparseArrays, StaticArrays, Test + +# Tests for issue #213: caches must be safe to reuse at a new x, regardless of +# how their internal scratch fields (e.g. JacobianCache.x1) were initialized. +# The original symptom was DI building a JacobianCache from `similar(x)` +# (uninitialized) and getting junk Jacobians in :central mode. + +const J_REF = [2.0 0.0; 0.0 3.0; 4.0 0.0] +foo_oop(x) = [2x[1], 3x[2], 4x[1]] +foo_iip!(y, x) = (y[1] = 2x[1]; y[2] = 3x[2]; y[3] = 4x[1]; y) + +# A non-zero point where the bug becomes obvious — at zeros(2) the junk in x1 +# cancels by symmetry on this affine function and hides the issue. +const X_TEST = [1.0, 2.0] + +"""Build a JacobianCache whose scratch fields are explicitly poisoned with a +huge value, mimicking what happens when a caller hands FiniteDiff a cache +allocated via `similar(x)` (which gives uninitialized memory).""" +function poisoned_jcache(fdtype; x_template = X_TEST, y_template = foo_oop(X_TEST), poison = 1.0e10) + x1 = fill(poison, length(x_template)) + fx = fill(poison, length(y_template)) + if fdtype === Val(:complex) + FiniteDiff.JacobianCache(x1, fx, nothing, fdtype) + else + fx1 = fill(poison, length(y_template)) + FiniteDiff.JacobianCache(x1, fx, fx1, fdtype) + end +end + +@testset "Cache reuse safety (issue #213)" begin + +@testset "JacobianCache out-of-place reuse" begin + @testset "fresh cache reused at new x" for fdtype in (Val(:forward), Val(:central), Val(:complex)) + cache = FiniteDiff.JacobianCache(zeros(2), zeros(3), fdtype) + # Exercise the cache once at x_old, then reuse at X_TEST. + FiniteDiff.finite_difference_jacobian(foo_oop, zeros(2), cache) + J = FiniteDiff.finite_difference_jacobian(foo_oop, X_TEST, cache) + @test J ≈ J_REF atol=1e-6 + end + + @testset "cache built with garbage x1/fx ($(fdtype))" for fdtype in (Val(:forward), Val(:central), Val(:complex)) + cache = poisoned_jcache(fdtype) + J = FiniteDiff.finite_difference_jacobian(foo_oop, X_TEST, cache) + @test J ≈ J_REF atol=1e-6 + end + + @testset "cache built with garbage x1/fx + sparsity ($(fdtype))" for fdtype in (Val(:forward), Val(:central)) + spJ = sparse(J_REF) + cache = poisoned_jcache(fdtype) + J = FiniteDiff.finite_difference_jacobian(foo_oop, X_TEST, cache; + sparsity = spJ, jac_prototype = spJ) + @test Matrix(J) ≈ J_REF atol=1e-6 + end +end + +@testset "JacobianCache in-place reuse" begin + @testset "fresh cache reused at new x ($(fdtype))" for fdtype in (Val(:forward), Val(:central), Val(:complex)) + cache = FiniteDiff.JacobianCache(zeros(2), zeros(3), fdtype) + J = zeros(3, 2) + FiniteDiff.finite_difference_jacobian!(J, foo_iip!, zeros(2), cache) + fill!(J, 0) + FiniteDiff.finite_difference_jacobian!(J, foo_iip!, X_TEST, cache) + @test J ≈ J_REF atol=1e-6 + end + + @testset "cache built with garbage x1/fx ($(fdtype))" for fdtype in (Val(:forward), Val(:central), Val(:complex)) + cache = poisoned_jcache(fdtype) + J = zeros(3, 2) + FiniteDiff.finite_difference_jacobian!(J, foo_iip!, X_TEST, cache) + @test J ≈ J_REF atol=1e-6 + end + + @testset "in-place :central must not mutate x (sparse path)" begin + spJ = sparse(J_REF) + cache = poisoned_jcache(Val(:central)) + J = zeros(3, 2) + x = copy(X_TEST) + x_orig = copy(x) + FiniteDiff.finite_difference_jacobian!(J, foo_iip!, x, cache; + sparsity = spJ, colorvec = 1:2) + @test Matrix(J) ≈ J_REF atol=1e-6 + @test x == x_orig # x should be restored / unmutated + end +end + +@testset "GradientCache reuse" begin + # `:complex` requires an analytic function, so don't use abs2 here. + g(x) = x[1]^2 + x[2]^2 + x[3]^2 + grad_ref = [2.0, 4.0, 6.0] + x = [1.0, 2.0, 3.0] + + # Use the allocating constructor so buffer types are correct, then poison + # any non-`nothing` buffer to simulate a stale cache. + @testset "vector → scalar with poisoned cache ($(fdtype))" for fdtype in (Val(:forward), Val(:central), Val(:complex)) + df = zeros(3) + cache = FiniteDiff.GradientCache(df, x, fdtype, Float64, Val(false)) + for fld in (:c1, :c2, :c3) + buf = getfield(cache, fld) + buf isa AbstractArray && fill!(buf, 1e10) + end + grad = zeros(3) + FiniteDiff.finite_difference_gradient!(grad, g, x, cache) + @test grad ≈ grad_ref atol=1e-5 + end + + @testset "fresh cache reused at new x ($(fdtype))" for fdtype in (Val(:forward), Val(:central)) + cache = FiniteDiff.GradientCache(zeros(3), zeros(3), fdtype) + grad = zeros(3) + FiniteDiff.finite_difference_gradient!(grad, g, zeros(3), cache) + FiniteDiff.finite_difference_gradient!(grad, g, x, cache) + @test grad ≈ grad_ref atol=1e-5 + end +end + +@testset "JVPCache reuse" begin + foo_iip!_3 = (y, x) -> (y[1] = 2x[1]; y[2] = 3x[2]; y[3] = 4x[1]; y) + v = [1.0, 0.0] + jvp_ref = J_REF * v + + @testset "garbage cache ($(fdtype))" for fdtype in (Val(:forward), Val(:central)) + x1 = fill(1e10, 2) + fx1 = fill(1e10, 3) + cache = FiniteDiff.JVPCache(x1, fx1, fdtype) + jvp = zeros(3) + FiniteDiff.finite_difference_jvp!(jvp, foo_iip!_3, X_TEST, v, cache) + @test jvp ≈ jvp_ref atol=1e-6 + end +end + +@testset "HessianCache reuse" begin + h(x) = x[1]^2 + 2 * x[2]^2 + H_ref = [2.0 0.0; 0.0 4.0] + + xpp = fill(1e10, 2); xpm = fill(1e10, 2); xmp = fill(1e10, 2); xmm = fill(1e10, 2) + cache = FiniteDiff.HessianCache(xpp, xpm, xmp, xmm, Val(:hcentral), Val(true)) + H = zeros(2, 2) + FiniteDiff.finite_difference_hessian!(H, h, X_TEST, cache) + @test H ≈ H_ref atol=1e-3 +end + +# Mirrors the failure mode from JuliaDiff/DifferentiationInterface.jl#983: a +# caller building a cache with `similar(x)` fields and then asking for a +# Jacobian via the non-allocating entry point. +@testset "DI-style similar() cache (issue #983 reproduction)" begin + foo(x) = [2x[1], 3x[2], 4x[1]] + y = foo(X_TEST) + + @testset "$(fdtype)" for fdtype in (Val(:forward), Val(:central), Val(:complex)) + x1 = similar(X_TEST) + fx = similar(y) + if fdtype === Val(:complex) + cache = FiniteDiff.JacobianCache(x1, fx, nothing, fdtype) + else + fx1 = similar(y) + cache = FiniteDiff.JacobianCache(x1, fx, fx1, fdtype) + end + J = FiniteDiff.finite_difference_jacobian(foo, X_TEST, cache) + @test J ≈ J_REF atol=1e-6 + end +end + +end # outer testset diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml index 71f4223..00a6714 100644 --- a/test/downstream/Project.toml +++ b/test/downstream/Project.toml @@ -1,4 +1,6 @@ [deps] +ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce" diff --git a/test/downstream/ordinarydiffeq_tridiagonal_solve.jl b/test/downstream/ordinarydiffeq_tridiagonal_solve.jl index 9032a92..ec44ef0 100644 --- a/test/downstream/ordinarydiffeq_tridiagonal_solve.jl +++ b/test/downstream/ordinarydiffeq_tridiagonal_solve.jl @@ -1,4 +1,4 @@ -using OrdinaryDiffEq, ForwardDiff, LinearAlgebra, Test +using OrdinaryDiffEq, OrdinaryDiffEqRosenbrock, ADTypes, ForwardDiff, LinearAlgebra, Test const nknots = 10 const h = 1.0/(nknots+1) @@ -21,8 +21,11 @@ sol_true = solve(prob, Rodas4P(), saveat=0.1) function loss(p) _prob = remake(prob, p=p) - sol = solve(_prob, Rodas4P(autodiff=false), saveat=0.1) + sol = solve(_prob, Rodas4P(autodiff=AutoFiniteDiff()), saveat=0.1) sum((sol .- sol_true).^2) end -@test ForwardDiff.gradient(loss, [1.0])[1] ≈ 0.6645766813735486 +# Loose tolerance: this is a smoke test that FiniteDiff works through a +# Rosenbrock solver with a Tridiagonal jacobian prototype; the exact value +# drifts with solver internals across OrdinaryDiffEq releases. +@test ForwardDiff.gradient(loss, [1.0])[1] ≈ 0.665 atol=1e-2 diff --git a/test/runtests.jl b/test/runtests.jl index ba2b432..fc66d8a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,7 @@ if GROUP == "All" || GROUP == "Core" @time @safetestset "FiniteDiff Standard Tests" begin include("finitedifftests.jl") end @time @safetestset "Color Differentiation Tests" begin include("coloring_tests.jl") end @time @safetestset "Out of Place Tests" begin include("out_of_place_tests.jl") end + @time @safetestset "Cache Reuse Safety Tests" begin include("cache_reuse_tests.jl") end end if GROUP == "All" || GROUP == "Downstream"