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

issue with Zygote jacobian when accessing ODE solutions on the time dimension #1102

Open
daviehh opened this issue Aug 30, 2024 · 2 comments
Open

issue with Zygote jacobian when accessing ODE solutions on the time dimension #1102

daviehh opened this issue Aug 30, 2024 · 2 comments
Assignees
@DhairyaLGandhi
Labels
bug

Comments

@daviehh
Copy link

daviehh commented Aug 30, 2024 •
edited
Loading

Describe the bug 🐞

Consider the example at doc, but one wishes to take the jacobian for all parameters, at a given time. The example works when taking the 1st element of the solution, at all times

solve(_prob, Tsit5(), reltol = 1e-6, abstol = 1e-6, saveat = 1)[1, :]

which does work for both ForwardDiff and Zygote. However, when one wishes to calculate the jacobian for all elements but at a given time, i.e. sol[:, i] instead of sol[i, :] as in the example, ForwardDiff gives the result but Zygote can give all zeros depending on how it's accessed: full code

using OrdinaryDiffEq
using Zygote
using SciMLSensitivity, ForwardDiff

function lotka_volterra!(du, u, p, t)
    du[1] = dx = p[1] * u[1] - p[2] * u[1] * u[2]
    du[2] = dy = -p[3] * u[2] + p[4] * u[1] * u[2]
end

p = [1.5, 1.0, 3.0, 1.0];
u0 = [1.0; 1.0];
prob = ODEProblem(lotka_volterra!, u0, (0.0, 10.0), p)
sol0 = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6)

function f_sol2(x, i)
    _prob = remake(prob, u0=x[1:2], p=x[3:end])
    sol = solve(_prob, Tsit5(), reltol=1e-6, abstol=1e-6, saveat=1)
    sol[:, i]
end

u0p = [u0; p]
df_zyg = Zygote.jacobian(x -> f_sol2(x, 6), u0p)[1]
df_fwd = ForwardDiff.jacobian(x -> f_sol2(x, 6), u0p)

in this case, df_zyg is all zero, while df_fwd works. Same if one access the interpolating sol(t) interface. However, if one access the solution at given time-step with the sol.u interface, Zygote works

function f_sol_u(x, i)
    _prob = remake(prob, u0=x[1:2], p=x[3:end])
    sol = solve(_prob, Tsit5(), reltol=1e-6, abstol=1e-6, saveat=1)
    sol.u[i]
end

df_zyg = Zygote.jacobian(x -> f_sol_u(x, 6), u0p)[1]
df_fwd = ForwardDiff.jacobian(x -> f_sol_u(x, 6), u0p)

which is a bit strange since sol.u[i] and sol[:, i] is 'hard equal ( === )', i.e. points to the same memory sol0.u[3] === sol0[:, 3] or check with pointer(sol0.u[3])

Expected behavior

Zygote.jacobian work with sol[:, i]

Minimal Reproducible Example 👇

using OrdinaryDiffEq
using Zygote
using SciMLSensitivity, ForwardDiff

function lotka_volterra!(du, u, p, t)
    du[1] = dx = p[1] * u[1] - p[2] * u[1] * u[2]
    du[2] = dy = -p[3] * u[2] + p[4] * u[1] * u[2]
end

p = [1.5, 1.0, 3.0, 1.0];
u0 = [1.0; 1.0];
prob = ODEProblem(lotka_volterra!, u0, (0.0, 10.0), p)
sol0 = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6)

function f_sol2(x, i)
    _prob = remake(prob, u0=x[1:2], p=x[3:end])
    sol = solve(_prob, Tsit5(), reltol=1e-6, abstol=1e-6, saveat=1)
    sol[:, i]
end

function f_sol_u(x, i)
    _prob = remake(prob, u0=x[1:2], p=x[3:end])
    sol = solve(_prob, Tsit5(), reltol=1e-6, abstol=1e-6, saveat=1)
    sol.u[i]
end

u0p = [u0; p]
df_zyg = Zygote.jacobian(x -> f_sol2(x, 6), u0p)[1]
df_fwd = ForwardDiff.jacobian(x -> f_sol2(x, 6), u0p)

df_zyg = Zygote.jacobian(x -> f_sol_u(x, 6), u0p)[1]
df_fwd = ForwardDiff.jacobian(x -> f_sol_u(x, 6), u0p)

sol0.u[3] === sol0[:, 3]

(added note)
works when returning:

@views sol[:, i]

or

sol.u[i]

wrong Zygote.jacobian when returning

sol[:, i]

or

sol(i)

Error & Stacktrace ⚠️

N/A: error is giving 0 for non-zero results

Environment (please complete the following information):

  • Output of using Pkg; Pkg.status()
Status `~/test_jl/environments/v1.10/Project.toml`
  [f6369f11] ForwardDiff v0.10.36
  [1dea7af3] OrdinaryDiffEq v6.89.0
  [1ed8b502] SciMLSensitivity v7.65.0
  [e88e6eb3] Zygote v0.6.70
  • Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)
Status `~/test_jl/environments/v1.10/Manifest.toml`
  [47edcb42] ADTypes v1.7.1
  [621f4979] AbstractFFTs v1.5.0
  [7d9f7c33] Accessors v0.1.37
  [79e6a3ab] Adapt v4.0.4
  [66dad0bd] AliasTables v1.1.3
  [ec485272] ArnoldiMethod v0.4.0
  [4fba245c] ArrayInterface v7.15.0
  [4c555306] ArrayLayouts v1.10.3
  [a9b6321e] Atomix v0.1.0
  [62783981] BitTwiddlingConvenienceFunctions v0.1.6
  [fa961155] CEnum v0.5.0
  [2a0fbf3d] CPUSummary v0.2.6
  [7057c7e9] Cassette v0.3.13
  [082447d4] ChainRules v1.69.0
  [d360d2e6] ChainRulesCore v1.24.0
  [fb6a15b2] CloseOpenIntervals v0.1.13
  [38540f10] CommonSolve v0.2.4
  [bbf7d656] CommonSubexpressions v0.3.1
  [f70d9fcc] CommonWorldInvalidations v1.0.0
  [34da2185] Compat v4.16.0
  [a33af91c] CompositionsBase v0.1.2
  [2569d6c7] ConcreteStructs v0.2.3
  [187b0558] ConstructionBase v1.5.7
  [adafc99b] CpuId v0.3.1
  [9a962f9c] DataAPI v1.16.0
  [864edb3b] DataStructures v0.18.20
  [e2d170a0] DataValueInterfaces v1.0.0
  [2b5f629d] DiffEqBase v6.154.0
  [459566f4] DiffEqCallbacks v3.8.0
  [77a26b50] DiffEqNoiseProcess v5.23.0
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.15.1
  [a0c0ee7d] DifferentiationInterface v0.5.15
  [b4f34e82] Distances v0.10.11
  [31c24e10] Distributions v0.25.111
  [ffbed154] DocStringExtensions v0.9.3
  [da5c29d0] EllipsisNotation v1.8.0
  [4e289a0a] EnumX v1.0.4
  [7da242da] Enzyme v0.12.33
  [f151be2c] EnzymeCore v0.7.8
  [d4d017d3] ExponentialUtilities v1.26.1
  [e2ba6199] ExprTools v0.1.10
⌅ [6b7a57c9] Expronicon v0.8.5
  [7034ab61] FastBroadcast v0.3.5
  [9aa1b823] FastClosures v0.3.2
  [29a986be] FastLapackInterface v2.0.4
  [1a297f60] FillArrays v1.13.0
  [6a86dc24] FiniteDiff v2.24.0
  [f6369f11] ForwardDiff v0.10.36
  [f62d2435] FunctionProperties v0.1.2
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v0.1.3
  [d9f16b24] Functors v0.4.12
  [0c68f7d7] GPUArrays v10.3.1
  [46192b85] GPUArraysCore v0.1.6
⌃ [61eb1bfa] GPUCompiler v0.27.1
  [c145ed77] GenericSchur v0.5.4
  [86223c79] Graphs v1.11.2
  [3e5b6fbb] HostCPUFeatures v0.1.17
  [34004b35] HypergeometricFunctions v0.3.24
  [7869d1d1] IRTools v0.4.14
  [615f187c] IfElse v0.1.1
  [d25df0c9] Inflate v0.1.5
  [3587e190] InverseFunctions v0.1.16
  [92d709cd] IrrationalConstants v0.2.2
  [82899510] IteratorInterfaceExtensions v1.0.0
  [692b3bcd] JLLWrappers v1.5.0
  [ccbc3e58] JumpProcesses v9.13.4
  [ef3ab10e] KLU v0.6.0
  [63c18a36] KernelAbstractions v0.9.25
  [ba0b0d4f] Krylov v0.9.6
⌅ [929cbde3] LLVM v8.1.0
  [10f19ff3] LayoutPointers v0.1.17
  [5078a376] LazyArrays v2.2.0
  [2d8b4e74] LevyArea v1.0.0
  [d3d80556] LineSearches v7.3.0
  [7ed4a6bd] LinearSolve v2.33.0
  [2ab3a3ac] LogExpFunctions v0.3.28
  [bdcacae8] LoopVectorization v0.12.171
  [d8e11817] MLStyle v0.4.17
  [1914dd2f] MacroTools v0.5.13
  [d125e4d3] ManualMemory v0.1.8
  [bb5d69b7] MaybeInplace v0.1.4
  [e1d29d7a] Missings v1.2.0
  [46d2c3a1] MuladdMacro v0.2.4
  [d41bc354] NLSolversBase v7.8.3
  [2774e3e8] NLsolve v4.5.1
  [872c559c] NNlib v0.9.22
  [77ba4419] NaNMath v1.0.2
  [8913a72c] NonlinearSolve v3.14.0
  [d8793406] ObjectFile v0.4.2
  [6fe1bfb0] OffsetArrays v1.14.1
  [429524aa] Optim v1.9.4
  [3bd65402] Optimisers v0.3.3
  [bac558e1] OrderedCollections v1.6.3
  [1dea7af3] OrdinaryDiffEq v6.89.0
  [89bda076] OrdinaryDiffEqAdamsBashforthMoulton v1.1.0
  [6ad6398a] OrdinaryDiffEqBDF v1.1.0
  [bbf590c4] OrdinaryDiffEqCore v1.3.0
  [50262376] OrdinaryDiffEqDefault v1.1.0
  [4302a76b] OrdinaryDiffEqDifferentiation v1.1.0
  [9286f039] OrdinaryDiffEqExplicitRK v1.1.0
  [e0540318] OrdinaryDiffEqExponentialRK v1.1.0
  [becaefa8] OrdinaryDiffEqExtrapolation v1.1.0
  [5960d6e9] OrdinaryDiffEqFIRK v1.1.0
  [101fe9f7] OrdinaryDiffEqFeagin v1.1.0
  [d3585ca7] OrdinaryDiffEqFunctionMap v1.1.0
  [d28bc4f8] OrdinaryDiffEqHighOrderRK v1.1.0
  [9f002381] OrdinaryDiffEqIMEXMultistep v1.1.0
  [521117fe] OrdinaryDiffEqLinear v1.1.0
  [1344f307] OrdinaryDiffEqLowOrderRK v1.1.0
  [b0944070] OrdinaryDiffEqLowStorageRK v1.2.0
  [127b3ac7] OrdinaryDiffEqNonlinearSolve v1.1.0
  [c9986a66] OrdinaryDiffEqNordsieck v1.1.0
  [5dd0a6cf] OrdinaryDiffEqPDIRK v1.1.0
  [5b33eab2] OrdinaryDiffEqPRK v1.1.0
  [04162be5] OrdinaryDiffEqQPRK v1.1.0
  [af6ede74] OrdinaryDiffEqRKN v1.1.0
  [43230ef6] OrdinaryDiffEqRosenbrock v1.1.0
  [2d112036] OrdinaryDiffEqSDIRK v1.1.0
  [669c94d9] OrdinaryDiffEqSSPRK v1.2.0
  [e3e12d00] OrdinaryDiffEqStabilizedIRK v1.1.0
  [358294b1] OrdinaryDiffEqStabilizedRK v1.1.0
  [fa646aed] OrdinaryDiffEqSymplecticRK v1.1.0
  [b1df2697] OrdinaryDiffEqTsit5 v1.1.0
  [79d7bb75] OrdinaryDiffEqVerner v1.1.0
  [90014a1f] PDMats v0.11.31
  [65ce6f38] PackageExtensionCompat v1.0.2
  [d96e819e] Parameters v0.12.3
  [e409e4f3] PoissonRandom v0.4.4
  [f517fe37] Polyester v0.7.16
  [1d0040c9] PolyesterWeave v0.2.2
  [85a6dd25] PositiveFactorizations v0.2.4
  [d236fae5] PreallocationTools v0.4.24
  [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [43287f4e] PtrArrays v1.2.1
  [1fd47b50] QuadGK v2.10.1
  [74087812] Random123 v1.7.0
  [e6cf234a] RandomNumbers v1.6.0
  [c1ae055f] RealDot v0.1.0
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v3.27.0
  [f2c3362d] RecursiveFactorization v0.2.23
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.0
  [ae5879a3] ResettableStacks v1.1.1
  [37e2e3b7] ReverseDiff v1.15.3
  [79098fc4] Rmath v0.7.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.13
  [94e857df] SIMDTypes v0.1.0
  [476501e8] SLEEFPirates v0.6.43
  [0bca4576] SciMLBase v2.50.6
  [c0aeaf25] SciMLOperators v0.3.10
  [1ed8b502] SciMLSensitivity v7.65.0
  [53ae85a6] SciMLStructures v1.4.2
  [6c6a2e73] Scratch v1.2.1
  [efcf1570] Setfield v1.1.1
  [727e6d20] SimpleNonlinearSolve v1.12.0
  [699a6c99] SimpleTraits v0.9.4
  [ce78b400] SimpleUnPack v1.1.0
  [a2af1166] SortingAlgorithms v1.2.1
  [47a9eef4] SparseDiffTools v2.20.0
  [dc90abb0] SparseInverseSubset v0.1.2
  [0a514795] SparseMatrixColorings v0.4.0
  [e56a9233] Sparspak v0.3.9
  [276daf66] SpecialFunctions v2.4.0
  [aedffcd0] Static v1.1.1
  [0d7ed370] StaticArrayInterface v1.8.0
  [90137ffa] StaticArrays v1.9.7
  [1e83bf80] StaticArraysCore v1.4.3
  [82ae8749] StatsAPI v1.7.0
  [2913bbd2] StatsBase v0.34.3
  [4c63d2b9] StatsFuns v1.3.1
  [789caeaf] StochasticDiffEq v6.68.0
  [7792a7ef] StrideArraysCore v0.5.7
  [09ab397b] StructArrays v0.6.18
  [53d494c1] StructIO v0.3.1
  [2efcf032] SymbolicIndexingInterface v0.3.28
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.12.0
  [8290d209] ThreadingUtilities v0.5.2
  [a759f4b9] TimerOutputs v0.5.24
  [9f7883ad] Tracker v0.2.35
  [d5829a12] TriangularSolve v0.2.1
  [410a4b4d] Tricks v0.1.9
  [781d530d] TruncatedStacktraces v1.4.0
  [3a884ed6] UnPack v1.0.2
  [013be700] UnsafeAtomics v0.2.1
  [d80eeb9a] UnsafeAtomicsLLVM v0.2.1
  [3d5dd08c] VectorizationBase v0.21.70
  [19fa3120] VertexSafeGraphs v0.2.0
  [e88e6eb3] Zygote v0.6.70
  [700de1a5] ZygoteRules v0.2.5
  [7cc45869] Enzyme_jll v0.0.145+0
  [1d5cc7b8] IntelOpenMP_jll v2024.2.1+0
⌅ [dad2f222] LLVMExtra_jll v0.0.31+0
  [856f044c] MKL_jll v2024.2.0+0
  [efe28fd5] OpenSpecFun_jll v0.5.5+0
⌅ [f50d1b31] Rmath_jll v0.4.3+0
  [1317d2d5] oneTBB_jll v2021.12.0+0
  [0dad84c5] ArgTools v1.1.1
  [56f22d72] Artifacts
  [2a0f44e3] Base64
  [ade2ca70] Dates
  [8ba89e20] Distributed
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching
  [9fa8497b] Future
  [b77e0a4c] InteractiveUtils
  [4af54fe1] LazyArtifacts
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2
  [8f399da3] Libdl
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging
  [d6f4376e] Markdown
  [a63ad114] Mmap
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.10.0
  [de0858da] Printf
  [3fa0cd96] REPL
  [9a3f8284] Random
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization
  [1a1011a3] SharedArrays
  [6462fe0b] Sockets
  [2f01184e] SparseArrays v1.10.0
  [10745b16] Statistics v1.10.0
  [4607b0f0] SuiteSparse
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [8dfed614] Test
  [cf7118a7] UUIDs
  [4ec0a83e] Unicode
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [deac9b47] LibCURL_jll v8.4.0+0
  [e37daf67] LibGit2_jll v1.6.4+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [c8ffd9c3] MbedTLS_jll v2.28.2+1
  [14a3606d] MozillaCACerts_jll v2023.1.10
  [4536629a] OpenBLAS_jll v0.3.23+4
  [05823500] OpenLibm_jll v0.8.1+2
  [bea87d4a] SuiteSparse_jll v7.2.1+1
  [83775a58] Zlib_jll v1.2.13+1
  [8e850b90] libblastrampoline_jll v5.11.0+0
  [8e850ede] nghttp2_jll v1.52.0+1
  [3f19e933] p7zip_jll v17.4.0+2
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`
  • Output of versioninfo()
Julia Version 1.10.5
Commit 6f3fdf7b362 (2024-08-27 14:19 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 10 × Apple M1 Pro
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)

Additional context

Add any other context about the problem here.
@daviehh daviehh added the bug label Aug 30, 2024
@daviehh daviehh changed the title issue with Zygote jacobian when accessing ODE solutions on the 2nd dimension issue with Zygote jacobian when accessing ODE solutions on the time dimension Sep 1, 2024
@ChrisRackauckas
Copy link
Member

@DhairyaLGandhi could you take a look?

@DhairyaLGandhi
Copy link
Member

Yeah I'll take a look at this.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@DhairyaLGandhi DhairyaLGandhi

Labels
bug
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@ChrisRackauckas @daviehh @DhairyaLGandhi