refactor: use DI.jl for the jacobian stuff

Project.toml

@@ -8,6 +8,7 @@ ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
@@ -58,20 +59,21 @@ NonlinearSolveSymbolicsExt = "Symbolics"
 NonlinearSolveZygoteExt = "Zygote"
 
 [compat]
-ADTypes = "1.1.0"
+ADTypes = "1.9"
 Aqua = "0.8"
-ArrayInterface = "7.9"
+ArrayInterface = "7.16"
 BandedMatrices = "1.5"
 BenchmarkTools = "1.4"
-CUDA = "5.2"
+CUDA = "5.5"
 ConcreteStructs = "0.2.3"
-DiffEqBase = "6.149.0"
-Enzyme = "0.12"
+DiffEqBase = "6.155.3"
+DifferentiationInterface = "0.6.1"
+Enzyme = "0.13.2"
 ExplicitImports = "1.5"
-FastBroadcast = "0.2.8, 0.3"
+FastBroadcast = "0.3.5"
 FastClosures = "0.3.2"
 FastLevenbergMarquardt = "0.1"
-FiniteDiff = "2.22"
+FiniteDiff = "2.24"
 FixedPointAcceleration = "0.3"
 ForwardDiff = "0.10.36"
 Hwloc = "3"
@@ -80,36 +82,36 @@ LazyArrays = "1.8.2, 2"
 LeastSquaresOptim = "0.8.5"
 LineSearches = "7.2"
 LinearAlgebra = "1.10"
-LinearSolve = "2.30"
+LinearSolve = "2.35"
 MINPACK = "1.2"
 MaybeInplace = "0.1.3"
-ModelingToolkit = "9.15.0"
+ModelingToolkit = "9.41.0"
 NLSolvers = "0.5"
 NLsolve = "4.5"
 NaNMath = "1"
 NonlinearProblemLibrary = "0.1.2"
-OrdinaryDiffEq = "6.75"
+OrdinaryDiffEqTsit5 = "1.1.0"
 Pkg = "1.10"
 PrecompileTools = "1.2"
 Preferences = "1.4"
 Printf = "1.10"
 Random = "1.91"
 ReTestItems = "1.24"
-RecursiveArrayTools = "3.8"
+RecursiveArrayTools = "3.27"
 Reexport = "1.2"
 SIAMFANLEquations = "1.0.1"
-SciMLBase = "2.34.0"
+SciMLBase = "2.54.0"
 SciMLJacobianOperators = "0.1"
-SimpleNonlinearSolve = "1.8"
+SimpleNonlinearSolve = "1.12.3"
 SparseArrays = "1.10"
-SparseDiffTools = "2.19"
+SparseDiffTools = "2.22"
 SpeedMapping = "0.3"
 StableRNGs = "1"
 StaticArrays = "1.9"
 StaticArraysCore = "1.4"
 Sundials = "4.23.1"
-SymbolicIndexingInterface = "0.3.15"
-Symbolics = "5.26, 6"
+SymbolicIndexingInterface = "0.3.31"
+Symbolics = "6.12"
 Test = "1.10"
 TimerOutputs = "0.5.23"
 Zygote = "0.6.69"
@@ -133,7 +135,7 @@ NLSolvers = "337daf1e-9722-11e9-073e-8b9effe078ba"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
@@ -147,4 +149,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "BandedMatrices", "BenchmarkTools", "CUDA", "Enzyme", "ExplicitImports", "FastLevenbergMarquardt", "FixedPointAcceleration", "Hwloc", "InteractiveUtils", "LeastSquaresOptim", "MINPACK", "ModelingToolkit", "NLSolvers", "NLsolve", "NaNMath", "NonlinearProblemLibrary", "OrdinaryDiffEq", "Pkg", "Random", "ReTestItems", "SIAMFANLEquations", "SpeedMapping", "StableRNGs", "StaticArrays", "Sundials", "Symbolics", "Test", "Zygote"]
+test = ["Aqua", "BandedMatrices", "BenchmarkTools", "CUDA", "Enzyme", "ExplicitImports", "FastLevenbergMarquardt", "FixedPointAcceleration", "Hwloc", "InteractiveUtils", "LeastSquaresOptim", "MINPACK", "ModelingToolkit", "NLSolvers", "NLsolve", "NaNMath", "NonlinearProblemLibrary", "OrdinaryDiffEqTsit5", "Pkg", "Random", "ReTestItems", "SIAMFANLEquations", "SpeedMapping", "StableRNGs", "StaticArrays", "Sundials", "Symbolics", "Test", "Zygote"]

docs/Project.toml

@@ -10,7 +10,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -30,10 +30,11 @@ DiffEqBase = "6.136"
 Documenter = "1"
 DocumenterCitations = "1"
 IncompleteLU = "0.2"
+InteractiveUtils = "<0.0.1, 1"
 LinearSolve = "2"
 ModelingToolkit = "8, 9"
 NonlinearSolve = "3"
-OrdinaryDiffEq = "6"
+OrdinaryDiffEqTsit5 = "1.1.0"
 Plots = "1"
 Random = "<0.0.1, 1"
 SciMLBase = "2.4"

docs/src/tutorials/code_optimization.md

@@ -33,9 +33,7 @@ Take for example a prototypical small nonlinear solver code in its out-of-place
 ```@example small_opt
 using NonlinearSolve
 
-function f(u, p)
-    u .* u .- p
-end
+f(u, p) = u .* u .- p
 u0 = [1.0, 1.0]
 p = 2.0
 prob = NonlinearProblem(f, u0, p)
@@ -53,9 +51,7 @@ using BenchmarkTools
 Note that this way of writing the function is a shorthand for:
 
 ```@example small_opt
-function f(u, p)
-    [u[1] * u[1] - p, u[2] * u[2] - p]
-end
+f(u, p) = [u[1] * u[1] - p, u[2] * u[2] - p]
 ```
 
 where the function `f` returns an array. This is a common pattern from things like MATLAB's
@@ -71,7 +67,7 @@ by hand, this looks like:
 function f(du, u, p)
     du[1] = u[1] * u[1] - p
     du[2] = u[2] * u[2] - p
-    nothing
+    return nothing
 end
 
 prob = NonlinearProblem(f, u0, p)
@@ -84,6 +80,7 @@ the `.=` in-place broadcasting.
 ```@example small_opt
 function f(du, u, p)
     du .= u .* u .- p
+    return nothing
 end
 
 @benchmark sol = solve(prob, NewtonRaphson())
@@ -114,6 +111,7 @@ to normal array expressions, for example:
 
 ```@example small_opt
 using StaticArrays
+
 A = SA[2.0, 3.0, 5.0]
 typeof(A)
 ```
@@ -135,22 +133,20 @@ want to use the out-of-place allocating form, but this time we want to output a
 array. Doing it with broadcasting looks like:
 
 ```@example small_opt
-function f_SA(u, p)
-    u .* u .- p
-end
+f_SA(u, p) = u .* u .- p
+
 u0 = SA[1.0, 1.0]
 p = 2.0
 prob = NonlinearProblem(f_SA, u0, p)
+
 @benchmark solve(prob, NewtonRaphson())
 ```
 
 Note that only change here is that `u0` is made into a StaticArray! If we needed to write
 `f` out for a more complex nonlinear case, then we'd simply do the following:
 
 ```@example small_opt
-function f_SA(u, p)
-    SA[u[1] * u[1] - p, u[2] * u[2] - p]
-end
+f_SA(u, p) = SA[u[1] * u[1] - p, u[2] * u[2] - p]
 
 @benchmark solve(prob, NewtonRaphson())
 ```

docs/src/tutorials/optimizing_parameterized_ode.md

@@ -4,7 +4,7 @@ Let us fit a parameterized ODE to some data. We will use the Lotka-Volterra mode
 example. We will use Single Shooting to fit the parameters.
 
 ```@example parameterized_ode
-using OrdinaryDiffEq, NonlinearSolve, Plots
+using OrdinaryDiffEqTsit5, NonlinearSolve, Plots
 ```
 
 Let us simulate some real data from the Lotka-Volterra model.

lib/SciMLJacobianOperators/README.md

@@ -56,4 +56,4 @@ snormal_form * v  # Computes JᵀJ * v
 #  8.0
 #  8.0
 #  8.0
-```
+```

src/NonlinearSolve.jl

@@ -7,8 +7,8 @@ end
 using Reexport: @reexport
 using PrecompileTools: @compile_workload, @setup_workload
 
-using ADTypes: ADTypes, AutoFiniteDiff, AutoForwardDiff, AutoPolyesterForwardDiff,
-               AutoZygote, AutoEnzyme, AutoSparse
+using ADTypes: ADTypes, AbstractADType, AutoFiniteDiff, AutoForwardDiff,
+               AutoPolyesterForwardDiff, AutoZygote, AutoEnzyme, AutoSparse
 # FIXME: deprecated, remove in future
 using ADTypes: AutoSparseFiniteDiff, AutoSparseForwardDiff, AutoSparsePolyesterForwardDiff,
                AutoSparseZygote
@@ -22,6 +22,7 @@ using DiffEqBase: DiffEqBase, AbstractNonlinearTerminationMode,
                   NormTerminationMode, RelNormTerminationMode, RelSafeBestTerminationMode,
                   RelSafeTerminationMode, RelTerminationMode,
                   SimpleNonlinearSolveTerminationMode, SteadyStateDiffEqTerminationMode
+using DifferentiationInterface: DifferentiationInterface, Constant
 using FastBroadcast: @..
 using FastClosures: @closure
 using FiniteDiff: FiniteDiff
@@ -39,7 +40,7 @@ using Printf: @printf
 using Preferences: Preferences, @load_preference, @set_preferences!
 using RecursiveArrayTools: recursivecopy!, recursivefill!
 using SciMLBase: AbstractNonlinearAlgorithm, JacobianWrapper, AbstractNonlinearProblem,
-                 AbstractSciMLOperator, _unwrap_val, has_jac, isinplace, NLStats
+                 AbstractSciMLOperator, _unwrap_val, isinplace, NLStats
 using SciMLJacobianOperators: AbstractJacobianOperator, JacobianOperator, VecJacOperator,
                               JacVecOperator, StatefulJacobianOperator
 using SparseArrays: AbstractSparseMatrix, SparseMatrixCSC
@@ -55,6 +56,8 @@ using SymbolicIndexingInterface: SymbolicIndexingInterface, ParameterIndexingPro
 
 @reexport using SciMLBase, SimpleNonlinearSolve
 
+const DI = DifferentiationInterface
+
 # Type-Inference Friendly Check for Extension Loading
 is_extension_loaded(::Val) = false
 

src/globalization/line_search.jl

@@ -101,7 +101,8 @@ function __internal_init(
         args...; stats, internalnorm::IN = DEFAULT_NORM, kwargs...) where {F, IN}
     T = promote_type(eltype(fu), eltype(u))
     if u isa Number
-        autodiff = get_concrete_forward_ad(alg.autodiff, prob; check_forward_mode = true)
+        autodiff = get_dense_ad(get_concrete_forward_ad(
+            alg.autodiff, prob; check_forward_mode = true))
         if !(autodiff isa AutoForwardDiff ||
              autodiff isa AutoPolyesterForwardDiff ||
              autodiff isa AutoFiniteDiff)
@@ -110,10 +111,10 @@ function __internal_init(
             @warn "Scalar AD is supported only for AutoForwardDiff and AutoFiniteDiff. \
                    Detected $(autodiff). Falling back to AutoFiniteDiff."
         end
-        deriv_op = @closure (du, u, fu, p) -> last(__value_derivative(
-                                                  autodiff, Base.Fix2(f, p), u)) *
-                                              fu *
-                                              du
+        deriv_op = @closure (du, u, fu, p) -> begin
+            # Temporary solution, we are anyways moving to LineSearch.jl
+            return DI.derivative(f, autodiff, u, Constant(p)) * fu * du
+        end
     else
         # Both forward and reverse AD can be used for line-search.
         # We prefer forward AD for better performance, however, reverse AD is also supported if user explicitly requests it.