SciML / DiffEqBase.jl

The lightweight Base library for shared types and functionality for defining differential equation and scientific machine learning (SciML) problems
Other
311 stars 111 forks source link

Norm of nested `ForwardDiff.Dual`s is broken #1023

Closed devmotion closed 5 months ago

devmotion commented 5 months ago

Describe the bug 🐞

https://github.com/SciML/DiffEqBase.jl/pull/981 broke ODE_DEFAULT_NORM for arrays with nested ForwardDiff.Duals which since then returns a ForwardDiff.Dual (one layer of Duals removed). This is problematic e.g. in OrdinaryDiffEq where integrator.EEst is integrator.eigen_est is initialized as a non-Dual number and updated with output of the time-dependent norm which defaults to ODE_DEFAULT_NORM.

More concretely, the main problem in #981 is that the sums are initialized with T(0) where T is the type of the value of the Dual number which for nested Dual numbers is again a Dual. Probably a better initialization would be something like DiffEqBase.value(T(0)) or DiffEqBase.value(zero(T)) - or, consistent with the sum(sse, ...) sse(T(0))/sse(zero(T)).

Expected behavior

It should return a non-Dual number.

Minimal Reproducible Example 👇

julia> using DiffEqBase, ForwardDiff

julia> DiffEqBase.ODE_DEFAULT_NORM(ForwardDiff.Dual{:b}.(ForwardDiff.Dual{:a}.([1.0, 2.0, 3.0], true), true), ForwardDiff.Dual{:b}(4.0, true))

On DiffEqBase 6.144.2 and master (and probably all versions in-between) the call returns

Dual{:a}(1.2909944487358056,0.0)

whereas in 6.144.1 (prior to the linked PR) it returned

1.2909944487358056

Environment (please complete the following information):

julia> versioninfo()
Julia Version 1.10.2
Commit bd47eca2c8a (2024-03-01 10:14 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 10 × Apple M2 Pro
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 6 virtual cores)
Environment:
  JULIA_PKG_USE_CLI_GIT = true
ChrisRackauckas commented 5 months ago

Change it from T(0) to value(T)(0)?

devmotion commented 5 months ago

I think actually that won't work 😄 https://github.com/SciML/DiffEqBase.jl/blob/ec7561abbe0185d0ed7903f5e0b850478fa875fe/src/forwarddiff.jl#L306-L307

The suggestions in the issue should work though, I'm about to prepare a PR.