feat: Introducing vr_aggregator with VRDirectCB and VRFRMODE

[sivasathyaseeelan]

using DiffEqBase, JumpProcesses, OrdinaryDiffEq, StochasticDiffEq, Test, Plots
using Random, LinearSolve
using StableRNGs
rng = StableRNG(12345)

rate = (u, p, t) -> u[1]
affect! = (integrator) -> (integrator.u[1] = integrator.u[1] / 2)
jump = VariableRateJump(rate, affect!, interp_points = 1000)
jump2 = deepcopy(jump)

f = function (du, u, p, t)
    du[1] = u[1]
end

prob = ODEProblem(f, [0.2], (0.0, 10.0))
jump_prob = JumpProblem(prob, Direct(), jump, jump2; rng = rng)

integrator = init(jump_prob, Tsit5())

sol = solve(jump_prob, Tsit5())
plot(sol)

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[isaacsas]

@sivasathyaseeelan I'm traveling the next few days, but will try to take a look Sunday if I have a chance, or next Tuesday.

[ChrisRackauckas]

This should be moved to SciMLBenchmarks.jl

[ChrisRackauckas]

These should get functionwrapped if they aren't tuples?

[ChrisRackauckas]

Should get an issue.

[ChrisRackauckas]

what is this quantity? Because it's not correctly the cumulative rate: it doesn't take into account condition rejections. But it's also not used.

[ChrisRackauckas]

It should probably just be deleted?

[ChrisRackauckas]

Why is this evaluated here instead of in the affect!?

[ChrisRackauckas]

Move this there, it's just unnecessary on most steps.

[ChrisRackauckas]

This has a bug because cache.current_threshold is never set, so you'll never actually accumulate.

[ChrisRackauckas]

Suggested change

	cache.prev_threshold = cache.current_threshold
	cache.current_threshold = randexp(rng)
	cache.current_threshold = randexp(rng)
	cache.prev_threshold = cache.current_threshold

[ChrisRackauckas]

Suggested change

	return cache.prev_threshold - rate_increment
	cache.current_threshold = cache.prev_threshold - rate_increment # cache increment if not zeroed in this round
	return cache.current_threshold

[ChrisRackauckas]

It's impossible for this to generally work without actually updating the current threshold https://github.com/SciML/JumpProcesses.jl/pull/477/files#r2052157167, so make a test case that stresses that and fails and then after bringing in that change it should pass. You should be able to see this if you do something like set dtmax = 0.0001, as you will get zero jumps in that case since the rates never accumulate.