Replace internal ITP implementation with SimpleNonlinearSolve.jl's ITP

This PR removes the internal `InternalITP` implementation from DiffEqBase
and replaces it with the more robust and feature-complete `ITP` algorithm
from SimpleNonlinearSolve.jl.

## Changes:

1. **Removed internal ITP implementation:**
   - Deleted `src/internal_itp.jl` containing `InternalITP` struct and solver
   - Removed helper functions like `prevfloat_tdir`, `nextfloat_tdir`, `max_tdir`

2. **Added SimpleNonlinearSolve dependency:**
   - Added SimpleNonlinearSolve as a dependency in Project.toml
   - Set minimum version bound to 2.7 to ensure compatibility

3. **Updated usage throughout codebase:**
   - Modified `src/callbacks.jl` to use `ITP()` instead of `InternalITP()`
   - Updated `src/DiffEqBase.jl` to import `ITP` from SimpleNonlinearSolve
   - Updated `ext/DiffEqBaseForwardDiffExt.jl` to handle the new ITP algorithm
   - Updated `test/internal_rootfinder.jl` to test with the new ITP implementation

## Benefits:

- **Better maintenance:** Eliminates duplicate code and reduces maintenance burden
- **Improved reliability:** Uses the well-tested and optimized ITP implementation from SimpleNonlinearSolve
- **Better performance:** SimpleNonlinearSolve's ITP has additional optimizations and tuning parameters
- **Consistency:** Aligns with the broader SciML ecosystem's approach of using specialized packages

## Testing:

All existing tests pass, including the specific internal rootfinder tests that verify
the ITP algorithm works correctly for various edge cases.

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <[email protected]>

Author: ChrisRackauckas

Project.toml

@@ -24,6 +24,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
+SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
@@ -97,6 +98,7 @@ ReverseDiff = "1"
 SciMLBase = "2.115.0"
 SciMLOperators = "1"
 SciMLStructures = "1.5"
+SimpleNonlinearSolve = "2.7"
 Setfield = "1"
 SparseArrays = "1.9"
 Static = "1"

ext/DiffEqBaseForwardDiffExt.jl

@@ -1,13 +1,13 @@
 module DiffEqBaseForwardDiffExt
 
 using DiffEqBase, ForwardDiff
+using SimpleNonlinearSolve: ITP
 using DiffEqBase.ArrayInterface
 using DiffEqBase: Void, FunctionWrappersWrappers, OrdinaryDiffEqTag,
                   AbstractTimeseriesSolution,
                   RecursiveArrayTools, reduce_tup, _promote_tspan, has_continuous_callback
 import DiffEqBase: hasdualpromote, wrapfun_oop, wrapfun_iip, prob2dtmin,
-                   promote_tspan, ODE_DEFAULT_NORM,
-                   InternalITP, nextfloat_tdir
+                   promote_tspan, ODE_DEFAULT_NORM
 import SciMLBase: isdualtype, DualEltypeChecker, sse, __sum
 
 const dualT = ForwardDiff.Dual{ForwardDiff.Tag{OrdinaryDiffEqTag, Float64}, Float64, 1}
@@ -153,7 +153,7 @@ end
 
 # Differentiation of internal solver
 
-function scalar_nlsolve_ad(prob, alg::InternalITP, args...; kwargs...)
+function scalar_nlsolve_ad(prob, alg::ITP, args...; kwargs...)
     f = prob.f
     p = value(prob.p)
 
@@ -186,7 +186,7 @@ end
 function SciMLBase.solve(
         prob::IntervalNonlinearProblem{uType, iip,
             <:ForwardDiff.Dual{T, V, P}},
-        alg::InternalITP, args...;
+        alg::ITP, 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),
@@ -202,7 +202,7 @@ function SciMLBase.solve(
                 V,
                 P},
             }},
-        alg::InternalITP, args...;
+        alg::ITP, args...;
         kwargs...) where {uType, iip, T, V, P}
     sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
 

src/DiffEqBase.jl

@@ -44,6 +44,8 @@ using SciMLBase
 
 using SciMLOperators: AbstractSciMLOperator, AbstractSciMLScalarOperator
 
+using SimpleNonlinearSolve: ITP
+
 using SciMLBase: @def, DEIntegrator, AbstractDEProblem,
                  AbstractDiffEqInterpolation,
                  DECallback, AbstractDEOptions, DECache, AbstractContinuousCallback,
@@ -140,7 +142,6 @@ include("utils.jl")
 include("stats.jl")
 include("calculate_residuals.jl")
 include("tableaus.jl")
-include("internal_itp.jl")
 
 include("callbacks.jl")
 include("common_defaults.jl")

src/callbacks.jl

@@ -358,17 +358,17 @@ end
 # rough implementation, needs multiple type handling
 # always ensures that if r = bisection(f, (x0, x1))
 # then either f(nextfloat(r)) == 0 or f(nextfloat(r)) * f(r) < 0
-# note: not really using bisection - uses the ITP method
+# note: not really using bisection - uses the ITP method  
 function bisection(
         f, tup, t_forward::Bool, rootfind::SciMLBase.RootfindOpt, abstol, reltol;
         maxiters = 1000)
     if rootfind == SciMLBase.LeftRootFind
         solve(IntervalNonlinearProblem{false}(f, tup),
-            InternalITP(), abstol = abstol,
+            ITP(), abstol = abstol,
             reltol = reltol).left
     else
         solve(IntervalNonlinearProblem{false}(f, tup),
-            InternalITP(), abstol = abstol,
+            ITP(), abstol = abstol,
             reltol = reltol).right
     end
 end

src/internal_itp.jl

@@ -1,8 +1,9 @@
 using DiffEqBase
-using DiffEqBase: InternalITP, IntervalNonlinearProblem
+using DiffEqBase: IntervalNonlinearProblem
+using SimpleNonlinearSolve: ITP  
 using ForwardDiff
 
-for Rootfinder in (InternalITP,)
+for Rootfinder in (ITP,)
     rf = Rootfinder()
     # From SimpleNonlinearSolve
     f = (u, p) -> u * u - p