WIP: Make ForwardDiff Optional

```julia
@time using OrdinaryDiffEqTsit5
# 0.678347 seconds (806.13 k allocations: 60.106 MiB, 3.97% gc time, 7.38% compilation time)

@time begin
  using OrdinaryDiffEqTsit5
  function lorenz(du,u,p,t)
      du[1] = 10.0(u[2]-u[1])
      du[2] = u[1]*(28.0-u[3]) - u[2]
      du[3] = u[1]*u[2] - (8/3)*u[3]
  end
  u0 = [1.0;0.0;0.0]; tspan = (0.0,100.0)
  prob = ODEProblem(lorenz,u0,tspan)
  solve(prob,Tsit5())
end

#  0.765026 seconds (1.29 M allocations: 86.179 MiB, 2.16% gc time, 24.20% compilation time: 2% of which was recompilation)

@time_imports using OrdinaryDiffEqTsit5

#=
               ┌ 0.0 ms DocStringExtensions.__init__()
      1.4 ms  DocStringExtensions
      0.5 ms  Reexport
      7.6 ms  Preferences
      0.6 ms  PrecompileTools
      0.5 ms  FastPower
      1.2 ms  ArrayInterface
      0.9 ms  StaticArraysCore
      0.7 ms  ArrayInterface → ArrayInterfaceStaticArraysCoreExt
      5.6 ms  FunctionWrappers
      0.5 ms  MuladdMacro
      3.2 ms  OrderedCollections
      0.5 ms  UnPack
      0.7 ms  Parameters
      0.8 ms  Statistics
      0.4 ms  IfElse
      0.9 ms  CommonWorldInvalidations
      9.8 ms  Static
      0.8 ms  Compat
      0.5 ms  Compat → CompatLinearAlgebraExt
     12.2 ms  StaticArrayInterface
      1.0 ms  ManualMemory
               ┌ 0.0 ms ThreadingUtilities.__init__()
      4.7 ms  ThreadingUtilities
      0.7 ms  SIMDTypes
      2.0 ms  LayoutPointers
      2.9 ms  CloseOpenIntervals
      8.5 ms  StrideArraysCore
      0.6 ms  BitTwiddlingConvenienceFunctions
               ┌ 0.0 ms CPUSummary.__init__()
      1.6 ms  CPUSummary
               ┌ 0.0 ms PolyesterWeave.__init__()
      2.6 ms  PolyesterWeave
      0.8 ms  Polyester
      1.5 ms  FastBroadcast
      9.4 ms  RecipesBase
      0.8 ms  ExprTools
      1.1 ms  Serialization
      0.9 ms  RuntimeGeneratedFunctions
      7.2 ms  MacroTools
               ┌ 0.0 ms InverseFunctions.__init__()
      1.6 ms  InverseFunctions
      0.7 ms  ConstructionBase
      0.6 ms  ConstructionBase → ConstructionBaseLinearAlgebraExt
      0.6 ms  CompositionsBase
      0.6 ms  CompositionsBase → CompositionsBaseInverseFunctionsExt
      0.6 ms  InverseFunctions → InverseFunctionsDatesExt
               ┌ 0.0 ms Accessors.__init__()
     10.6 ms  Accessors
      0.7 ms  Accessors → LinearAlgebraExt
      2.9 ms  SymbolicIndexingInterface
      0.8 ms  Adapt
      0.7 ms  DataValueInterfaces
      1.0 ms  DataAPI
      0.7 ms  IteratorInterfaceExtensions
      0.7 ms  TableTraits
      8.2 ms  Tables
      0.9 ms  GPUArraysCore
      0.6 ms  ArrayInterface → ArrayInterfaceGPUArraysCoreExt
     17.1 ms  RecursiveArrayTools
      0.8 ms  RecursiveArrayTools → RecursiveArrayToolsFastBroadcastExt
               ┌ 0.0 ms TruncatedStacktraces.__init__()
      0.8 ms  TruncatedStacktraces
      6.6 ms  Setfield
     14.8 ms  IrrationalConstants
      1.2 ms  DiffRules
      1.3 ms  DiffResults
               ┌ 2.1 ms OpenLibm_jll.__init__()
      4.9 ms  OpenLibm_jll
      0.9 ms  NaNMath
      0.9 ms  LogExpFunctions
      0.8 ms  JLLWrappers
               ┌ 12.0 ms CompilerSupportLibraries_jll.__init__()
     15.4 ms  CompilerSupportLibraries_jll
               ┌ 1.3 ms OpenSpecFun_jll.__init__()
      2.3 ms  OpenSpecFun_jll
      5.7 ms  SpecialFunctions
      0.8 ms  CommonSubexpressions
     25.3 ms  ForwardDiff
      0.8 ms  LogExpFunctions → LogExpFunctionsInverseFunctionsExt
      0.7 ms  FastPower → FastPowerForwardDiffExt
      0.7 ms  RecursiveArrayTools → RecursiveArrayToolsForwardDiffExt
      0.9 ms  EnumX
      0.8 ms  ConcreteStructs
      0.9 ms  FastClosures
      1.1 ms  PreallocationTools
      0.8 ms  FunctionWrappersWrappers
      0.9 ms  SciMLStructures
               ┌ 0.0 ms Distributed.__init__()
      9.3 ms  Distributed
      0.8 ms  CommonSolve
      3.1 ms  ADTypes
      1.2 ms  ADTypes → ADTypesConstructionBaseExt
     75.6 ms  MLStyle 20.54% compilation time
      4.6 ms  Expronicon
      8.7 ms  SciMLOperators
      1.0 ms  SciMLOperators → SciMLOperatorsStaticArraysCoreExt
               ┌ 0.0 ms SciMLBase.__init__()
     76.9 ms  SciMLBase
      5.1 ms  DiffEqBase
     38.0 ms  FillArrays
      1.3 ms  FillArrays → FillArraysStatisticsExt
      0.9 ms  SimpleUnPack
     21.6 ms  DataStructures
     31.1 ms  OrdinaryDiffEqCore
     39.0 ms  OrdinaryDiffEqTsit5
=#
```

This shows that ForwardDiff and its dependencies are a large part of the startup time, and that comes through DiffEqBase. With the DifferentiationInterface changes we will hopefully no longer default to ForwardDiff anyways and instead with Enzyme, and so we should make the stack not require ForwardDiff unless the user/solver wants it.

Author: ChrisRackauckas

Project.toml

@@ -13,15 +13,13 @@ EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b"
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
-PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
@@ -40,6 +38,7 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
 GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
@@ -55,6 +54,7 @@ DiffEqBaseCUDAExt = "CUDA"
 DiffEqBaseChainRulesCoreExt = "ChainRulesCore"
 DiffEqBaseDistributionsExt = "Distributions"
 DiffEqBaseEnzymeExt = ["ChainRulesCore", "Enzyme"]
+DiffEqBaseForwardDiffExt = ["ForwardDiff"]
 DiffEqBaseGeneralizedGeneratedExt = "GeneralizedGenerated"
 DiffEqBaseGTPSAExt = "GTPSA"
 DiffEqBaseMPIExt = "MPI"
@@ -92,7 +92,6 @@ Measurements = "2"
 MonteCarloMeasurements = "1"
 MuladdMacro = "0.2.1"
 Parameters = "0.12.0"
-PreallocationTools = "0.4"
 PrecompileTools = "1"
 Printf = "1.9"
 RecursiveArrayTools = "3"

src/forwarddiff.jl ext/DiffEqBaseForwardDiffExt.jlsrc/forwarddiff.jl renamed to ext/DiffEqBaseForwardDiffExt.jl

@@ -1,5 +1,99 @@
+module DiffEqBaseForwardDiffExt
+
+using DiffEqBase, ForwardDiff
+using DiffeqBase: Void, FunctionWrappersWrappers, OrdinaryDiffEqTag
+import DiffEqBase: hasdualpromote, wrapfun_oop, wrapfun_iip, promote_u0
+
 const DUALCHECK_RECURSION_MAX = 10
 
+eltypedual(x)  = eltype(x) <: ForwardDiff.Dual
+const dualT = ForwardDiff.Dual{ForwardDiff.Tag{OrdinaryDiffEqTag, Float64}, Float64, 1}
+dualgen(::Type{T}) where {T} = ForwardDiff.Dual{ForwardDiff.Tag{OrdinaryDiffEqTag, T}, T, 1}
+
+hasdualpromote(u0,t::Number) = hasmethod(ArrayInterface.promote_eltype,
+                    Tuple{Type{typeof(u0)}, Type{dualgen(eltype(u0))}}) &&
+                hasmethod(promote_rule,
+                    Tuple{Type{eltype(u0)}, Type{dualgen(eltype(u0))}}) &&
+                hasmethod(promote_rule,
+                    Tuple{Type{eltype(u0)}, Type{typeof(t)}}))
+
+const NORECOMPILE_IIP_SUPPORTED_ARGS = (
+    Tuple{Vector{Float64}, Vector{Float64},
+        Vector{Float64}, Float64},
+    Tuple{Vector{Float64}, Vector{Float64},
+        SciMLBase.NullParameters, Float64})
+
+const oop_arglists = (Tuple{Vector{Float64}, Vector{Float64}, Float64},
+    Tuple{Vector{Float64}, SciMLBase.NullParameters, Float64},
+    Tuple{Vector{Float64}, Vector{Float64}, dualT},
+    Tuple{Vector{dualT}, Vector{Float64}, Float64},
+    Tuple{Vector{dualT}, SciMLBase.NullParameters, Float64},
+    Tuple{Vector{Float64}, SciMLBase.NullParameters, dualT})
+
+const NORECOMPILE_OOP_SUPPORTED_ARGS = (Tuple{Vector{Float64},
+        Vector{Float64}, Float64},
+    Tuple{Vector{Float64},
+        SciMLBase.NullParameters, Float64})
+const oop_returnlists = (Vector{Float64}, Vector{Float64},
+    ntuple(x -> Vector{dualT}, length(oop_arglists) - 2)...)
+
+function wrapfun_oop(ff, inputs::Tuple = ())
+    if !isempty(inputs)
+        IT = Tuple{map(typeof, inputs)...}
+        if IT ∉ NORECOMPILE_OOP_SUPPORTED_ARGS
+            throw(NoRecompileArgumentError(IT))
+        end
+    end
+    FunctionWrappersWrappers.FunctionWrappersWrapper(ff, oop_arglists,
+        oop_returnlists)
+end
+
+function wrapfun_iip(ff,
+        inputs::Tuple{T1, T2, T3, T4}) where {T1, T2, T3, T4}
+    T = eltype(T2)
+    dualT = dualgen(T)
+    dualT1 = ArrayInterface.promote_eltype(T1, dualT)
+    dualT2 = ArrayInterface.promote_eltype(T2, dualT)
+    dualT4 = dualgen(promote_type(T, T4))
+
+    iip_arglists = (Tuple{T1, T2, T3, T4},
+        Tuple{dualT1, dualT2, T3, T4},
+        Tuple{dualT1, T2, T3, dualT4},
+        Tuple{dualT1, dualT2, T3, dualT4})
+
+    iip_returnlists = ntuple(x -> Nothing, 4)
+
+    fwt = map(iip_arglists, iip_returnlists) do A, R
+        FunctionWrappersWrappers.FunctionWrappers.FunctionWrapper{R, A}(Void(ff))
+    end
+    FunctionWrappersWrappers.FunctionWrappersWrapper{typeof(fwt), false}(fwt)
+end
+
+const iip_arglists_default = (
+    Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64},
+        Float64},
+    Tuple{Vector{Float64}, Vector{Float64},
+        SciMLBase.NullParameters,
+        Float64
+    },
+    Tuple{Vector{dualT}, Vector{Float64}, Vector{Float64}, dualT},
+    Tuple{Vector{dualT}, Vector{dualT}, Vector{Float64}, dualT},
+    Tuple{Vector{dualT}, Vector{dualT}, Vector{Float64}, Float64},
+    Tuple{Vector{dualT}, Vector{dualT}, SciMLBase.NullParameters,
+        Float64
+    },
+    Tuple{Vector{dualT}, Vector{Float64},
+        SciMLBase.NullParameters, dualT
+    })
+const iip_returnlists_default = ntuple(x -> Nothing, length(iip_arglists_default))
+
+function wrapfun_iip(@nospecialize(ff))
+    fwt = map(iip_arglists_default, iip_returnlists_default) do A, R
+        FunctionWrappersWrappers.FunctionWrappers.FunctionWrapper{R, A}(Void(ff))
+    end
+    FunctionWrappersWrappers.FunctionWrappersWrapper{typeof(fwt), false}(fwt)
+end
+
 """
     promote_dual(::Type{T},::Type{T2})
 
@@ -371,8 +465,6 @@ end
 anyeltypedual(::@Kwargs{}, ::Type{Val{counter}} = Val{0}) where {counter} = Any
 anyeltypedual(::Type{@Kwargs{}}, ::Type{Val{counter}} = Val{0}) where {counter} = Any
 
-@inline promote_u0(::Nothing, p, t0) = nothing
-
 @inline function promote_u0(u0, p, t0)
     if SciMLStructures.isscimlstructure(p)
         _p = SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]
@@ -482,6 +574,69 @@ end
 
 # Static Arrays don't support the `init` keyword argument for `sum`
 @inline __sum(f::F, args...; init, kwargs...) where {F} = sum(f, args...; init, kwargs...)
-@inline function __sum(f::F, a::StaticArraysCore.StaticArray...; init, kwargs...) where {F}
+@inline function __sum(f::F, a::DiffEqBase.StaticArraysCore.StaticArray...; init, kwargs...) where {F}
     return mapreduce(f, +, a...; init, kwargs...)
 end
+
+# Differentiation of internal solver
+
+function scalar_nlsolve_ad(prob, alg::InternalITP, args...; kwargs...)
+    f = prob.f
+    p = value(prob.p)
+
+    if prob isa IntervalNonlinearProblem
+        tspan = value(prob.tspan)
+        newprob = IntervalNonlinearProblem(f, tspan, p; prob.kwargs...)
+    else
+        u0 = value(prob.u0)
+        newprob = NonlinearProblem(f, u0, p; prob.kwargs...)
+    end
+
+    sol = solve(newprob, alg, args...; kwargs...)
+
+    uu = sol.u
+    if p isa Number
+        f_p = ForwardDiff.derivative(Base.Fix1(f, uu), p)
+    else
+        f_p = ForwardDiff.gradient(Base.Fix1(f, uu), p)
+    end
+
+    f_x = ForwardDiff.derivative(Base.Fix2(f, p), uu)
+    pp = prob.p
+    sumfun = let f_x′ = -f_x
+        ((fp, p),) -> (fp / f_x′) * ForwardDiff.partials(p)
+    end
+    partials = sum(sumfun, zip(f_p, pp))
+    return sol, partials
+end
+
+function SciMLBase.solve(
+        prob::IntervalNonlinearProblem{uType, iip,
+            <:ForwardDiff.Dual{T, V, P}},
+        alg::InternalITP, args...;
+        kwargs...) where {uType, iip, T, V, P}
+    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
+    return SciMLBase.build_solution(prob, alg, ForwardDiff.Dual{T, V, P}(sol.u, partials),
+        sol.resid; retcode = sol.retcode,
+        left = ForwardDiff.Dual{T, V, P}(sol.left, partials),
+        right = ForwardDiff.Dual{T, V, P}(sol.right, partials))
+end
+
+function SciMLBase.solve(
+        prob::IntervalNonlinearProblem{uType, iip,
+            <:AbstractArray{
+                <:ForwardDiff.Dual{T,
+                V,
+                P},
+            }},
+        alg::InternalITP, args...;
+        kwargs...) where {uType, iip, T, V, P}
+    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
+
+    return SciMLBase.build_solution(prob, alg, ForwardDiff.Dual{T, V, P}(sol.u, partials),
+        sol.resid; retcode = sol.retcode,
+        left = ForwardDiff.Dual{T, V, P}(sol.left, partials),
+        right = ForwardDiff.Dual{T, V, P}(sol.right, partials))
+end
+
+end

src/DiffEqBase.jl

@@ -32,18 +32,13 @@ import TruncatedStacktraces
 
 using Setfield
 
-using ForwardDiff
-
 using EnumX
 
 using Markdown
 
 using ConcreteStructs: @concrete
 using FastClosures: @closure
 
-# Could be made optional/glue
-import PreallocationTools
-
 import FunctionWrappersWrappers
 
 using SciMLBase
@@ -111,6 +106,13 @@ Reexport.@reexport using SciMLBase
 
 SciMLBase.isfunctionwrapper(x::FunctionWrapper) = true
 
+## Extension Functions
+
+eltypedual(x) = false
+promote_u0(::Nothing, p, t0) = nothing
+
+## Types
+
 """
 $(TYPEDEF)
 """
@@ -132,7 +134,6 @@ include("utils.jl")
 include("stats.jl")
 include("calculate_residuals.jl")
 include("tableaus.jl")
-include("internal_falsi.jl")
 include("internal_itp.jl")
 
 include("callbacks.jl")