Fix support for BigFloat u0

[ErikQQY]

Fix #43
Unblock #178
Possibly need SciML/PreallocationTools.jl#72

I think get_tmp in PreallocationTools.jl doesn't support BigFloat and would cause cannot reinterpret `BigFloat` as `ForwardDiff.Dual`, type `ForwardDiff.Dual` is not a bits type error, similar issue is SciML/PreallocationTools.jl#72

@avik-pal

Full stacktrace

ERROR: ArgumentError: cannot reinterpret `BigFloat` as `ForwardDiff.Dual{ForwardDiff.Tag{BoundaryValueDiffEq.__Fix3{BoundaryValueDiffEq.var"#131#137"{SciMLBase.StandardBVProblem, BoundaryValueDiffEq.MIRKCache{true, BigFloat, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(simplependulum!), typeof(bc!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing, Nothing, Dict{Any, Any}}}, typeof(bc!), BVProblem{Vector{BigFloat}, Tuple{Float64, Float64}, true, Nothing, SciMLBase.NullParameters, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(simplependulum!), typeof(bc!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing, Nothing, Dict{Any, Any}}}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}}, SciMLBase.StandardBVProblem, SciMLBase.NullParameters, MIRK4{Nothing, BVPJacobianAlgorithm{AutoForwardDiff{nothing, Nothing}, AutoSparse{AutoForwardDiff{nothing, Nothing}, ADTypes.NoSparsityDetector, ADTypes.NoColoringAlgorithm}, AutoSparse{AutoForwardDiff{nothing, Nothing}, ADTypes.NoSparsityDetector, ADTypes.NoColoringAlgorithm}}, Float64}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{BigFloat}, Vector{BigFloat}, Vector{BigFloat}, Matrix{BigFloat}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{BigFloat}, Vector{BigFloat}, Vector{BigFloat}, BigFloat}, Vector{BigFloat}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{BigFloat}, Vector{BigFloat}}}, Vector{Matrix{BigFloat}}, Vector{PreallocationTools.DiffCache{Vector{BigFloat}, Vector{BigFloat}}}, Vector{Vector{BigFloat}}, Vector{PreallocationTools.DiffCache{Vector{BigFloat}, Vector{BigFloat}}}, PreallocationTools.DiffCache{Vector{BigFloat}, Vector{BigFloat}}, Vector{BigFloat}, Vector{Vector{BigFloat}}, Vector{Vector{BigFloat}}, Tuple{Int64}, @NamedTuple{abstol::Float64, dt::Float64, adaptive::Bool}}}, SciMLBase.NullParameters}, BigFloat}, BigFloat, 11}`, type `ForwardDiff.Dual{ForwardDiff.Tag{BoundaryValueDiffEq.__Fix3{BoundaryValueDiffEq.var"#131#137"{SciMLBase.StandardBVProblem, BoundaryValueDiffEq.MIRKCache{true, BigFloat, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(simplependulum!), typeof(bc!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing, Nothing, Dict{Any, Any}}}, typeof(bc!), BVProblem{Vector{BigFloat}, Tuple{Float64, Float64}, true, Nothing, SciMLBase.NullParameters, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(simplependulum!), typeof(bc!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing, Nothing, Dict{Any, Any}}}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}}, SciMLBase.StandardBVProblem, SciMLBase.NullParameters, MIRK4{Nothing, BVPJacobianAlgorithm{AutoForwardDiff{nothing, Nothing}, AutoSparse{AutoForwardDiff{nothing, Nothing}, ADTypes.NoSparsityDetector, ADTypes.NoColoringAlgorithm}, AutoSparse{AutoForwardDiff{nothing, Nothing}, ADTypes.NoSparsityDetector, ADTypes.NoColoringAlgorithm}}, Float64}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{BigFloat}, Vector{BigFloat}, Vector{BigFloat}, Matrix{BigFloat}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{BigFloat}, Vector{BigFloat}, Vector{BigFloat}, BigFloat}, Vector{BigFloat}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{BigFloat}, Vector{BigFloat}}}, Vector{Matrix{BigFloat}}, Vector{PreallocationTools.DiffCache{Vector{BigFloat}, Vector{BigFloat}}}, Vector{Vector{BigFloat}}, Vector{PreallocationTools.DiffCache{Vector{BigFloat}, Vector{BigFloat}}}, PreallocationTools.DiffCache{Vector{BigFloat}, Vector{BigFloat}}, Vector{BigFloat}, Vector{Vector{BigFloat}}, Vector{Vector{BigFloat}}, Tuple{Int64}, @NamedTuple{abstol::Float64, dt::Float64, adaptive::Bool}}}, SciMLBase.NullParameters}, BigFloat}, BigFloat, 11}` is not a bits type
Stacktrace:
  [1] (::Base.var"#throwbits#333")(S::Type, T::Type, U::Type)
    @ Base .\reinterpretarray.jl:16
  [2] reinterpret(::Type{ForwardDiff.Dual{ForwardDiff.Tag{…}, BigFloat, 11}}, a::SubArray{BigFloat, 1, Vector{BigFloat}, Tuple{UnitRange{…}}, true})
    @ Base .\reinterpretarray.jl:62
  [3] get_tmp(dc::PreallocationTools.DiffCache{Vector{BigFloat}, Vector{BigFloat}}, u::Vector{ForwardDiff.Dual{ForwardDiff.Tag{…}, BigFloat, 11}})
    @ PreallocationTools D:\Julia\.julia\packages\PreallocationTools\vWjVk\src\PreallocationTools.jl:152
  [4] (::BoundaryValueDiffEq.var"#3#4"{PreallocationTools.DiffCache{Vector{…}, Vector{…}}, Vector{ForwardDiff.Dual{…}}})()
    @ BoundaryValueDiffEq D:\SciML\BVP\bc\BoundaryValueDiffEq.jl\src\types.jl:169
  [5] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging .\logging.jl:515
  [6] with_logger(f::Function, logger::Base.CoreLogging.NullLogger)
    @ Base.CoreLogging .\logging.jl:627
  [7] get_tmp
    @ D:\SciML\BVP\bc\BoundaryValueDiffEq.jl\src\types.jl:168 [inlined]
  [8] _broadcast_getindex_evalf
    @ .\broadcast.jl:709 [inlined]
  [9] _broadcast_getindex
    @ .\broadcast.jl:682 [inlined]
 [10] getindex
    @ .\broadcast.jl:636 [inlined]
 [11] copy
    @ .\broadcast.jl:942 [inlined]
 [12] materialize
    @ .\broadcast.jl:903 [inlined]

[github-actions]

Benchmark Results

	master	40656db...	master/40656db0437483...
Simple Pendulum/IIP/BoundaryValueDiffEq.MIRK2()	6.43 ± 0.16 ms	6.59 ± 0.18 ms	0.975
Simple Pendulum/IIP/BoundaryValueDiffEq.MIRK3()	2.33 ± 0.076 ms	2.38 ± 0.08 ms	0.978
Simple Pendulum/IIP/BoundaryValueDiffEq.MIRK4()	0.804 ± 0.027 ms	0.822 ± 0.029 ms	0.978
Simple Pendulum/IIP/BoundaryValueDiffEq.MIRK5()	0.849 ± 0.03 ms	0.874 ± 0.037 ms	0.971
Simple Pendulum/IIP/BoundaryValueDiffEq.MIRK6()	0.957 ± 0.026 ms	0.979 ± 0.031 ms	0.978
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false)	1.32 ± 0.1 ms	1.33 ± 0.11 ms	0.99
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true)	2.57 ± 0.18 ms	2.56 ± 0.17 ms	1
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0348 ± 0.005 s	0.0349 ± 0.0049 s	0.997
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.114 ± 0.0039 s	0.117 ± 0.0057 s	0.974
Simple Pendulum/IIP/Shooting(Tsit5())	0.188 ± 0.0057 ms	0.186 ± 0.0054 ms	1.01
Simple Pendulum/OOP/BoundaryValueDiffEq.MIRK2()	0.0398 ± 0.0031 s	0.0396 ± 0.0029 s	1
Simple Pendulum/OOP/BoundaryValueDiffEq.MIRK3()	11.5 ± 0.14 ms	11.6 ± 0.18 ms	0.992
Simple Pendulum/OOP/BoundaryValueDiffEq.MIRK4()	3.31 ± 0.064 ms	3.36 ± 0.072 ms	0.985
Simple Pendulum/OOP/BoundaryValueDiffEq.MIRK5()	3.4 ± 0.067 ms	3.4 ± 0.079 ms	1
Simple Pendulum/OOP/BoundaryValueDiffEq.MIRK6()	3.52 ± 0.072 ms	3.54 ± 0.095 ms	0.996
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false)	3.21 ± 0.76 ms	3.31 ± 0.78 ms	0.971
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true)	5.84 ± 1.1 ms	5.88 ± 1.2 ms	0.992
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.124 ± 0.0089 s	0.124 ± 0.01 s	0.997
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.377 ± 0.0062 s	0.381 ± 0.0065 s	0.988
Simple Pendulum/OOP/Shooting(Tsit5())	0.742 ± 0.038 ms	0.744 ± 0.035 ms	0.998
time_to_load	5.73 ± 0.015 s	5.74 ± 0.056 s	0.998

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.
Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

[ChrisRackauckas]

__similar(x), why are we not just using zero(x)?

[ChrisRackauckas]

I think get_tmp in PreallocationTools.jl doesn't support BigFloat and would cause cannot reinterpret BigFloat as ForwardDiff.Dual, type ForwardDiff.Dual is not a bits type error, similar issue is SciML/PreallocationTools.jl#72

We should just use a LazyBufferCache then?

[ErikQQY]

Implemented a special diff cache for the BigFloat case, which may not be the best solution, and some improvement is still necessary, the BigFloat support for MIRK is working now!!

[ChrisRackauckas]

I don't understand the need for __init_bigfloat_array!!: aren't you always building BigFloat caches now even if you don't use BigFloats? Why isn't this a part of the DualCache? Can't you just use the Vector{Any} in the dualcache?

[ErikQQY]

I don't understand the need for __init_bigfloat_array!!: aren't you always building BigFloat caches now even if you don't use BigFloats? Why isn't this a part of the DualCache? Can't you just use the Vector{Any} in the dualcache?

We only build the BigFloat cache when the input u0 is BigFloat

[ChrisRackauckas]

SciML/PreallocationTools.jl#112 should be all you need?

[ChrisRackauckas]

Hopefully this can close 🤞 and BVP needs to do nothing special.

[ChrisRackauckas]

is there a reason to not always fill with zeros?

[ErikQQY]

zeros is better, just change them all. While NonlinearSolve is doing the similar thing like specialize on bigfloat and do __similar things in SciML/NonlinearSolve.jl#438, maybe we can change to zeros everywhere in NonlinearSolve.jl as well?

[ChrisRackauckas]

We should. https://github.com/SciML/SciMLStyle?tab=readme-ov-file#always-default-to-constructs-which-initialize-data

[ChrisRackauckas]

Note this will also fix issues with non-fully defined functions. For example, if a user doesn't write to du[i], we should still converge now, while before you'd get something random. So a few bugs should be fixed, some bugs avoided, etc. Avoiding uninitialized memory is just a good idea all around and NonlinearSolve needs to do this change too.