TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Bastin) doesn't work with Shooting

[avik-pal]

Describe the example

A clear and concise description of what the example is.

Minimal Reproducible Example 👇

Without MRE, we would only be able to help you to a limited extent, and attention to the issue would be limited. to know more about MRE refer to wikipedia and stackoverflow.

using BoundaryValueDiffEq, OrdinaryDiffEq, NonlinearSolve

f1(u, p, t) = [u[2], -u[1]]

function bc1(sol, p, t)
    t₁, t₂ = extrema(t)
    solₜ₁ = sol(t₁)
    solₜ₂ = sol(t₂)
    solₜ₃ = sol((t₁ + t₂) / 2)
    # We know that this overconstrained system has a solution
    return [solₜ₁[1], solₜ₂[1] - 1, solₜ₃[1] - 0.51735, solₜ₃[2] + 1.92533]
end

tspan = (0.0, 100.0)
u0 = [0.0, 1.0]

bvp1 = BVProblem(BVPFunction{false}(f1, bc1; bcresid_prototype = zeros(4)), u0, tspan;
    nlls = Val(true))

solver = Shooting(Tsit5(),
    TrustRegion(; autodiff = AutoForwardDiff(; chunksize = 2),
        radius_update_scheme = RadiusUpdateSchemes.Bastin))

solve(bvp1, solver; verbose = false, abstol = 1e-6, reltol = 1e-6,
    odesolve_kwargs = (; abstol = 1e-6, reltol = 1e-6))

Error & Stacktrace ⚠️

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{SparseDiffTools.DeivVecTag, Float64}, Float64, 1})

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  Float64(::Int16)
   @ Base float.jl:159
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{SparseDiffTools.DeivVecTag, Float64}, Float64, 1})
    @ Base ./number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 1}, i1::Int64)
    @ Base ./array.jl:1021
  [3] _unsafe_copyto!(dest::Vector{Float64}, doffs::Int64, src::Vector{ForwardDiff.Dual{…}}, soffs::Int64, n::Int64)
    @ Base ./array.jl:299
  [4] unsafe_copyto!
    @ ./array.jl:353 [inlined]
  [5] _copyto_impl!(dest::Array, doffs::Integer, src::Array, soffs::Integer, n::Integer)
    @ Base ./array.jl:376 [inlined]
  [6] copyto!
    @ ./array.jl:363 [inlined]
  [7] copyto!
    @ ./array.jl:385 [inlined]
  [8] copyto_axcheck!(dest::Any, src::Any)
    @ Base ./abstractarray.jl:1174 [inlined]
  [9] Vector{Float64}(x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{SparseDiffTools.DeivVecTag, Float64}, Float64, 1}})
    @ Base ./array.jl:673
 [10] convert
    @ ./array.jl:665 [inlined]
 [11] setproperty!(x::Any, f::Symbol, v::Any)
    @ Base ./Base.jl:40 [inlined]
 [12] reinit!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, u0::Vector{…}; t0::Float64, tf::Float64, erase_sol::Bool, tstops::Tuple{}, saveat::Tuple{}, d_discontinuities::Tuple{}, reset_dt::Bool, reinit_callbacks::Bool, initialize_save::Bool, reinit_cache::Bool, reinit_retcode::Bool)
    @ OrdinaryDiffEq /mnt/research/nlsolve/OrdinaryDiffEq.jl/src/integrators/integrator_interface.jl:366
 [13] reinit!
    @ /mnt/research/nlsolve/OrdinaryDiffEq.jl/src/integrators/integrator_interface.jl:350 [inlined]
 [14] __single_shooting_loss(u::Any, p::Any, cache::Any, bc::BC, u0_size::Any, pt::Any) where BC
    @ BoundaryValueDiffEq /mnt/research/nlsolve/BoundaryValueDiffEq.jl/src/solve/single_shooting.jl:114 [inlined]
 [15] #225
    @ /mnt/research/nlsolve/BoundaryValueDiffEq.jl/src/solve/single_shooting.jl:33 [inlined]
 [16] NonlinearFunction
    @ /mnt/research/nlsolve/SciMLBase.jl/src/scimlfunctions.jl:2369 [inlined]
 [17] JacobianWrapper
    @ /mnt/research/nlsolve/SciMLBase.jl/src/function_wrappers.jl:97 [inlined]
 [18] auto_jacvec(f::SciMLBase.JacobianWrapper{…}, x::Vector{…}, v::Vector{…})
    @ SparseDiffTools /mnt/research/nlsolve/SparseDiffTools.jl/src/differentiation/jaches_products.jl:32
 [19] #95
    @ /mnt/research/nlsolve/NonlinearSolve.jl/src/internal/operators.jl:114 [inlined]
 [20] JacobianOperator
    @ /mnt/research/nlsolve/NonlinearSolve.jl/src/internal/operators.jl:198 [inlined]
 [21] mul!
    @ /mnt/research/nlsolve/NonlinearSolve.jl/src/internal/operators.jl:234 [inlined]
 [22] __mul!
    @ /mnt/julia/packages/MaybeInplace/DQx9o/src/MaybeInplace.jl:304 [inlined]
 [23] macro expansion
    @ /mnt/julia/packages/MaybeInplace/DQx9o/src/MaybeInplace.jl:171 [inlined]
 [24] __internal_solve!(cache::NonlinearSolve.GenericTrustRegionSchemeCache{…}, J::Matrix{…}, fu::Vector{…}, u::Vector{…}, δu::Vector{…}, descent_stats::@NamedTuple{…})
    @ NonlinearSolve /mnt/research/nlsolve/NonlinearSolve.jl/src/globalization/trust_region.jl:519
 [25] __step!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…}; recompute_jacobian::Nothing, kwargs::@Kwargs{})
    @ NonlinearSolve /mnt/research/nlsolve/NonlinearSolve.jl/src/core/generalized_first_order.jl:244
 [26] __step!
    @ /mnt/research/nlsolve/NonlinearSolve.jl/src/core/generalized_first_order.jl:203 [inlined]
 [27] #step!#204
    @ /mnt/research/nlsolve/NonlinearSolve.jl/src/core/generic.jl:52 [inlined]
 [28] step!(cache::NonlinearSolve.AbstractNonlinearSolveCache{iip, timeit}, args::Vararg{Any}) where {iip, timeit}
    @ NonlinearSolve /mnt/research/nlsolve/NonlinearSolve.jl/src/core/generic.jl:48 [inlined]
 [29] solve!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…})
    @ NonlinearSolve /mnt/research/nlsolve/NonlinearSolve.jl/src/core/generic.jl:13
 [30] __solve(prob::Union{…}, alg::NonlinearSolve.AbstractNonlinearSolveAlgorithm, args::Vararg{…}; kwargs...)
    @ NonlinearSolve /mnt/research/nlsolve/NonlinearSolve.jl/src/core/generic.jl:4 [inlined]
 [31] __solve
    @ /mnt/research/nlsolve/NonlinearSolve.jl/src/core/generic.jl:1 [inlined]
 [32] __solve(prob::BVProblem{…}, alg_::Shooting{…}; odesolve_kwargs::@NamedTuple{…}, nlsolve_kwargs::@NamedTuple{}, verbose::Bool, kwargs::@Kwargs{…})
    @ BoundaryValueDiffEq /mnt/research/nlsolve/BoundaryValueDiffEq.jl/src/solve/single_shooting.jl:75
 [33] __solve(prob::BVProblem, alg_::Shooting)
    @ BoundaryValueDiffEq /mnt/research/nlsolve/BoundaryValueDiffEq.jl/src/solve/single_shooting.jl:1 [inlined]
 [34] solve_call(_prob::BVProblem{…}, args::Shooting{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase /mnt/julia/packages/DiffEqBase/eLhx9/src/solve.jl:609
 [35] solve_call
    @ DiffEqBase /mnt/julia/packages/DiffEqBase/eLhx9/src/solve.jl:567 [inlined]
 [36] #solve_up#42
    @ DiffEqBase /mnt/julia/packages/DiffEqBase/eLhx9/src/solve.jl:1058 [inlined]
 [37] solve_up
    @ DiffEqBase /mnt/julia/packages/DiffEqBase/eLhx9/src/solve.jl:1044 [inlined]
 [38] #solve#40
    @ DiffEqBase /mnt/julia/packages/DiffEqBase/eLhx9/src/solve.jl:981 [inlined]
 [39] top-level scope
    @ REPL[69]:1
Some type information was truncated. Use `show(err)` to see complete types.

Not Working Environment (please complete the following information):

NonlinearSolve -- 3.5+

Working Environment (please complete the following information):

NonlinearSolve -- <3.5

Additional context

This change was introduced because in the latest updates we made Bastin and Yuan more efficient by directly using JVP and VJPs instead of materializing the full jacobian as before. This causes a problem with the cached versions which try to reuse the ODE solver cache.