feat: integrate callback saving into integrator init and finalize

Project.toml

@@ -38,7 +38,7 @@ ADTypes = "1"
 Adapt = "3, 4"
 ArrayInterface = "6, 7"
 DataStructures = "0.18, 0.19"
-DiffEqBase = "6.154"
+DiffEqBase = "6.187"
 DiffEqNoiseProcess = "5.13"
 DocStringExtensions = "0.8, 0.9"
 FastPower = "1"
@@ -48,17 +48,18 @@ JumpProcesses = "9"
 LevyArea = "1.0.0"
 LinearAlgebra = "1.6"
 Logging = "1.6"
+ModelingToolkit = "10"
 MuladdMacro = "0.2.1"
 NLsolve = "4"
 OrdinaryDiffEq = "6.87"
-OrdinaryDiffEqCore = "1.12.1"
+OrdinaryDiffEqCore = "1.32.0"
 OrdinaryDiffEqDifferentiation = "1.9"
 OrdinaryDiffEqNonlinearSolve = "1"
 Random = "1.6"
 RandomNumbers = "1.5.3"
 RecursiveArrayTools = "2, 3"
 Reexport = "0.2, 1.0"
-SciMLBase = "2.71"
+SciMLBase = "2.115"
 SciMLOperators = "0.2.9, 0.3, 0.4, 1"
 SparseArrays = "1.6"
 StaticArrays = "0.11, 0.12, 1.0"

src/integrators/integrator_utils.jl

@@ -234,6 +234,7 @@ end
             copyat_or_push!(integrator.sol.u, integrator.saveiter,
                 integrator.u[integrator.opts.save_idxs], Val{false})
         end
+        SciMLBase.save_final_discretes!(integrator, integrator.opts.callback)
     end
     if (!isnothing(integrator.W) && integrator.W.curt != integrator.t) ||
        (!isnothing(integrator.P) && integrator.P.curt != integrator.t)

src/solve.jl

@@ -791,4 +791,8 @@ function initialize_callbacks!(integrator, initialize_save = true)
 
     # reset this as it is now handled so the integrators should proceed as normal
     integrator.u_modified = false
+
+    if initialize_save
+        SciMLBase.save_discretes_if_enabled!(integrator, integrator.opts.callback; skip_duplicates = true)
+    end
 end

test/reversal_tests.jl

@@ -124,10 +124,10 @@ end
     else
         prob_forward = remake(prob, noise = W_forward, u0 = vec(prob.u0))
     end
-    sys = modelingtoolkitize(prob_forward)
+    sys = complete(modelingtoolkitize(prob_forward))
     sys2 = stochastic_integral_transform(sys, -1 // 1)
-    fdrift = eval(generate_rhs(sys2)[i])
-    fdif = eval(generate_diffusion_function(sys2)[i])
+    fdrift = generate_rhs(sys2; expression = Val{false})[i]
+    fdif = generate_diffusion_function(sys2; expression = Val{false})[i]
 
     for solver in Ito_solver
         println("solver: ", solver)