Split the FIRK (Radau) generator to a separate package

[ChrisRackauckas]

It's pretty heavy in terms of dependencies and pretty niche, so it should definitely be an add-on. It's relatively easy to give an informative error about too.

CC @Shreyas-Ekanathan @oscardssmith

[oscardssmith]

would it change your mind if we removed the symbolics stack from it? that's only used for the error estimator, and we think it shouldn't be necessary

[ChrisRackauckas]

Maybe? Though relying on type piracy packages (GenericLinearAlgebra) probably shouldn't be in the main pathway either.

[oscardssmith]

all we need from that is the eimen transform of a matrix of BigFloat. is there a better way to get that?

[ChrisRackauckas]

I don't see the issue with subpackaging it though. With a solution to JuliaLang/julia#55516 we'd probably want to make it a separate feature anyways

[oscardssmith]

@ChrisRackauckas does this actually work? I think this PR just breaks AdaptiveRadau.

[ChrisRackauckas]

No, AdaptiveRadau passes the tests here?

[ChrisRackauckas]

When you say it breaks, what's the break?

[oscardssmith]

using OrdinaryDiffEqFIRK, OrdinaryDiffEqFirkGenerator
julia> solve(prob3, AdaptiveRadau(max_order=17))
ERROR: num_stages choice 9 out of the pre-generated tableau range. For the fully adaptive Radau, please load the extension via `using OrdinaryDiffEqFIRKGenerator`
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] adaptiveRadauTableau(T1::Type, T2::Type, num_stages::Int64)
    @ OrdinaryDiffEqFIRK ~/.julia/dev/OrdinaryDiffEq/lib/OrdinaryDiffEqFIRK/src/firk_tableaus.jl:532
  [3] alg_cache(alg::AdaptiveRadau{…}, u::Vector{…}, rate_prototype::Vector{…}, ::Type{…}, ::Type{…}, ::Type{…}, uprev::Vector{…}, uprev2::Vector{…}, f::ODEFunction{…}, t::Float64, dt::Float64, reltol::Float64, p::SciMLBase.NullParameters, calck::Bool, ::Val{…})
    @ OrdinaryDiffEqFIRK ~/.julia/dev/OrdinaryDiffEq/lib/OrdinaryDiffEqFIRK/src/firk_caches.jl:523
  [4] __init(prob::ODEProblem{…}, alg::AdaptiveRadau{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Rational{…}, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEqCore.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/V1HJl/src/solve.jl:360
  [5] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEqCore/V1HJl/src/solve.jl:11 [inlined]
  [6] __solve(::ODEProblem{…}, ::AdaptiveRadau{…}; kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/V1HJl/src/solve.jl:6
  [7] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/V1HJl/src/solve.jl:1 [inlined]
  [8] #solve_call#44
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:632 [inlined]
  [9] solve_call(_prob::ODEProblem{…}, args::AdaptiveRadau{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:589
 [10] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::AdaptiveRadau{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1120
 [11] solve_up
    @ ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1099 [inlined]
 [12] solve(prob::ODEProblem{…}, args::AdaptiveRadau{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1036
 [13] solve(prob::ODEProblem{…}, args::AdaptiveRadau{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HW4ge/src/solve.jl:1026
 [14] top-level scope
    @ REPL[13]:1
Some type information was truncated. Use `show(err)` to see complete types.

[ChrisRackauckas]

Interesting, the tests in #2539 ran perfectly fine?

[oscardssmith]

yes. We didn't have any tests for higher than 13th order.

[ChrisRackauckas]

It's testing 17 and 21 https://github.com/SciML/OrdinaryDiffEq.jl/pull/2539/files#diff-205da55b83537a4502451960c1b5fb4907be91225b11165c39aed8de61213e3bR9

[oscardssmith]

ah. never mind. I was testing with a slightly broken setup. I do really think that we might want to reconsider separating them though... With #2531 the only 2 extra dependencies needed are GenericSchur and FastGaussQuadrature with import times (from @time_imports) of 5.4 and 0.5ms respectively.

[ChrisRackauckas]

It's worth considering. Though GenericSchur IIRC is one of the type pirating libraries that gives some issues here and there in binary builds, and so one of the goals of the package split was to remove it by default by isolating it to the exponential integrators.