symbolic_events.jl

using ModelingToolkit, OrdinaryDiffEq, StochasticDiffEq, JumpProcesses, Test
using SciMLStructures: canonicalize, Discrete
using ModelingToolkit: SymbolicContinuousCallback,
                       SymbolicContinuousCallbacks, NULL_AFFECT,
                       get_callback,
                       t_nounits as t,
                       D_nounits as D
using StableRNGs
import SciMLBase
using SymbolicIndexingInterface
rng = StableRNG(12345)

@variables x(t) = 0

eqs = [D(x) ~ 1]
affect = [x ~ 0]
affect_neg = [x ~ 1]

## Test SymbolicContinuousCallback
@testset "SymbolicContinuousCallback constructors" begin
    e = SymbolicContinuousCallback(eqs[])
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT
    @test e.affect_neg == NULL_AFFECT
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT
    @test e.affect_neg == NULL_AFFECT
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs, NULL_AFFECT)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT
    @test e.affect_neg == NULL_AFFECT
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs[], NULL_AFFECT)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT
    @test e.affect_neg == NULL_AFFECT
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs => NULL_AFFECT)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT
    @test e.affect_neg == NULL_AFFECT
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs[] => NULL_AFFECT)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == NULL_AFFECT
    @test e.affect_neg == NULL_AFFECT
    @test e.rootfind == SciMLBase.LeftRootFind

    ## With affect

    e = SymbolicContinuousCallback(eqs[], affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs, affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs, affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs[], affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs => affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs[] => affect)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect
    @test e.rootfind == SciMLBase.LeftRootFind

    # with only positive edge affect

    e = SymbolicContinuousCallback(eqs[], affect, affect_neg = nothing)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test isnothing(e.affect_neg)
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs, affect, affect_neg = nothing)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test isnothing(e.affect_neg)
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs, affect, affect_neg = nothing)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test isnothing(e.affect_neg)
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs[], affect, affect_neg = nothing)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test isnothing(e.affect_neg)
    @test e.rootfind == SciMLBase.LeftRootFind

    # with explicit edge affects

    e = SymbolicContinuousCallback(eqs[], affect, affect_neg = affect_neg)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect_neg
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs, affect, affect_neg = affect_neg)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect_neg
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs, affect, affect_neg = affect_neg)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect_neg
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(eqs[], affect, affect_neg = affect_neg)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect_neg
    @test e.rootfind == SciMLBase.LeftRootFind

    # with different root finding ops

    e = SymbolicContinuousCallback(
        eqs[], affect, affect_neg = affect_neg, rootfind = SciMLBase.LeftRootFind)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect_neg
    @test e.rootfind == SciMLBase.LeftRootFind

    e = SymbolicContinuousCallback(
        eqs[], affect, affect_neg = affect_neg, rootfind = SciMLBase.RightRootFind)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect_neg
    @test e.rootfind == SciMLBase.RightRootFind

    e = SymbolicContinuousCallback(
        eqs[], affect, affect_neg = affect_neg, rootfind = SciMLBase.NoRootFind)
    @test e isa SymbolicContinuousCallback
    @test isequal(e.eqs, eqs)
    @test e.affect == affect
    @test e.affect_neg == affect_neg
    @test e.rootfind == SciMLBase.NoRootFind
    # test plural constructor

    e = SymbolicContinuousCallbacks(eqs[])
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == NULL_AFFECT

    e = SymbolicContinuousCallbacks(eqs)
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == NULL_AFFECT

    e = SymbolicContinuousCallbacks(eqs[] => affect)
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == affect

    e = SymbolicContinuousCallbacks(eqs => affect)
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == affect

    e = SymbolicContinuousCallbacks([eqs[] => affect])
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == affect

    e = SymbolicContinuousCallbacks([eqs => affect])
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == affect

    e = SymbolicContinuousCallbacks(SymbolicContinuousCallbacks([eqs => affect]))
    @test e isa Vector{SymbolicContinuousCallback}
    @test isequal(e[].eqs, eqs)
    @test e[].affect == affect
end

##

@named sys = ODESystem(eqs, t, continuous_events = [x ~ 1])
@test getfield(sys, :continuous_events)[] ==
      SymbolicContinuousCallback(Equation[x ~ 1], NULL_AFFECT)
@test isequal(equations(getfield(sys, :continuous_events))[], x ~ 1)
fsys = flatten(sys)
@test isequal(equations(getfield(fsys, :continuous_events))[], x ~ 1)

@named sys2 = ODESystem([D(x) ~ 1], t, continuous_events = [x ~ 2], systems = [sys])
@test getfield(sys2, :continuous_events)[] ==
      SymbolicContinuousCallback(Equation[x ~ 2], NULL_AFFECT)
@test all(ModelingToolkit.continuous_events(sys2) .== [
    SymbolicContinuousCallback(Equation[x ~ 2], NULL_AFFECT),
    SymbolicContinuousCallback(Equation[sys.x ~ 1], NULL_AFFECT)
])

@test isequal(equations(getfield(sys2, :continuous_events))[1], x ~ 2)
@test length(ModelingToolkit.continuous_events(sys2)) == 2
@test isequal(ModelingToolkit.continuous_events(sys2)[1].eqs[], x ~ 2)
@test isequal(ModelingToolkit.continuous_events(sys2)[2].eqs[], sys.x ~ 1)

sys = complete(sys)
sys_nosplit = complete(sys; split = false)
sys2 = complete(sys2)
# Functions should be generated for root-finding equations
prob = ODEProblem(sys, Pair[], (0.0, 2.0))
p0 = 0
t0 = 0
@test get_callback(prob) isa ModelingToolkit.DiffEqCallbacks.ContinuousCallback
cb = ModelingToolkit.generate_rootfinding_callback(sys)
cond = cb.condition
out = [0.0]
cond.rf_ip(out, [0], p0, t0)
@test out[] ≈ -1 # signature is u,p,t
cond.rf_ip(out, [1], p0, t0)
@test out[] ≈ 0  # signature is u,p,t
cond.rf_ip(out, [2], p0, t0)
@test out[] ≈ 1  # signature is u,p,t

prob = ODEProblem(sys, Pair[], (0.0, 2.0))
prob_nosplit = ODEProblem(sys_nosplit, Pair[], (0.0, 2.0))
sol = solve(prob, Tsit5())
sol_nosplit = solve(prob_nosplit, Tsit5())
@test minimum(t -> abs(t - 1), sol.t) < 1e-10 # test that the solver stepped at the root
@test minimum(t -> abs(t - 1), sol_nosplit.t) < 1e-10 # test that the solver stepped at the root

# Test that a user provided callback is respected
test_callback = DiscreteCallback(x -> x, x -> x)
prob = ODEProblem(sys, Pair[], (0.0, 2.0), callback = test_callback)
prob_nosplit = ODEProblem(sys_nosplit, Pair[], (0.0, 2.0), callback = test_callback)
cbs = get_callback(prob)
cbs_nosplit = get_callback(prob_nosplit)
@test cbs isa CallbackSet
@test cbs.discrete_callbacks[1] == test_callback
@test cbs_nosplit isa CallbackSet
@test cbs_nosplit.discrete_callbacks[1] == test_callback

prob = ODEProblem(sys2, Pair[], (0.0, 3.0))
cb = get_callback(prob)
@test cb isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback

cond = cb.condition
out = [0.0, 0.0]
# the root to find is 2
cond.rf_ip(out, [0, 0], p0, t0)
@test out[1] ≈ -2 # signature is u,p,t
cond.rf_ip(out, [1, 0], p0, t0)
@test out[1] ≈ -1  # signature is u,p,t
cond.rf_ip(out, [2, 0], p0, t0) # this should return 0
@test out[1] ≈ 0  # signature is u,p,t

# the root to find is 1
out = [0.0, 0.0]
cond.rf_ip(out, [0, 0], p0, t0)
@test out[2] ≈ -1 # signature is u,p,t
cond.rf_ip(out, [0, 1], p0, t0) # this should return 0
@test out[2] ≈ 0  # signature is u,p,t
cond.rf_ip(out, [0, 2], p0, t0)
@test out[2] ≈ 1  # signature is u,p,t

sol = solve(prob, Tsit5())
@test minimum(t -> abs(t - 1), sol.t) < 1e-10 # test that the solver stepped at the first root
@test minimum(t -> abs(t - 2), sol.t) < 1e-10 # test that the solver stepped at the second root

@named sys = ODESystem(eqs, t, continuous_events = [x ~ 1, x ~ 2]) # two root eqs using the same unknown
sys = complete(sys)
prob = ODEProblem(sys, Pair[], (0.0, 3.0))
@test get_callback(prob) isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback
sol = solve(prob, Tsit5())
@test minimum(t -> abs(t - 1), sol.t) < 1e-10 # test that the solver stepped at the first root
@test minimum(t -> abs(t - 2), sol.t) < 1e-10 # test that the solver stepped at the second root

## Test bouncing ball with equation affect
@variables x(t)=1 v(t)=0

root_eqs = [x ~ 0]
affect = [v ~ -v]

@named ball = ODESystem([D(x) ~ v
                         D(v) ~ -9.8], t, continuous_events = root_eqs => affect)

@test getfield(ball, :continuous_events)[] ==
      SymbolicContinuousCallback(Equation[x ~ 0], Equation[v ~ -v])
ball = structural_simplify(ball)

@test length(ModelingToolkit.continuous_events(ball)) == 1

tspan = (0.0, 5.0)
prob = ODEProblem(ball, Pair[], tspan)
sol = solve(prob, Tsit5())
@test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close
# plot(sol)

## Test bouncing ball in 2D with walls
@variables x(t)=1 y(t)=0 vx(t)=0 vy(t)=1

continuous_events = [[x ~ 0] => [vx ~ -vx]
                     [y ~ -1.5, y ~ 1.5] => [vy ~ -vy]]

@named ball = ODESystem(
    [D(x) ~ vx
     D(y) ~ vy
     D(vx) ~ -9.8
     D(vy) ~ -0.01vy], t; continuous_events)

ball = structural_simplify(ball)
ball_nosplit = structural_simplify(ball; split = false)

tspan = (0.0, 5.0)
prob = ODEProblem(ball, Pair[], tspan)
prob_nosplit = ODEProblem(ball_nosplit, Pair[], tspan)

cb = get_callback(prob)
@test cb isa ModelingToolkit.DiffEqCallbacks.VectorContinuousCallback
@test getfield(ball, :continuous_events)[1] ==
      SymbolicContinuousCallback(Equation[x ~ 0], Equation[vx ~ -vx])
@test getfield(ball, :continuous_events)[2] ==
      SymbolicContinuousCallback(Equation[y ~ -1.5, y ~ 1.5], Equation[vy ~ -vy])
cond = cb.condition
out = [0.0, 0.0, 0.0]
cond.rf_ip(out, [0, 0, 0, 0], p0, t0)
@test out ≈ [0, 1.5, -1.5]

sol = solve(prob, Tsit5())
sol_nosplit = solve(prob_nosplit, Tsit5())
@test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close
@test minimum(sol[y]) ≈ -1.5 # check wall conditions
@test maximum(sol[y]) ≈ 1.5  # check wall conditions
@test 0 <= minimum(sol_nosplit[x]) <= 1e-10 # the ball never went through the floor but got very close
@test minimum(sol_nosplit[y]) ≈ -1.5 # check wall conditions
@test maximum(sol_nosplit[y]) ≈ 1.5  # check wall conditions

# tv = sort([LinRange(0, 5, 200); sol.t])
# plot(sol(tv)[y], sol(tv)[x], line_z=tv)
# vline!([-1.5, 1.5], l=(:black, 5), primary=false)
# hline!([0], l=(:black, 5), primary=false)

## Test multi-variable affect
# in this test, there are two variables affected by a single event.
continuous_events = [
    [x ~ 0] => [vx ~ -vx, vy ~ -vy]
]

@named ball = ODESystem([D(x) ~ vx
                         D(y) ~ vy
                         D(vx) ~ -1
                         D(vy) ~ 0], t; continuous_events)

ball_nosplit = structural_simplify(ball)
ball = structural_simplify(ball)

tspan = (0.0, 5.0)
prob = ODEProblem(ball, Pair[], tspan)
prob_nosplit = ODEProblem(ball_nosplit, Pair[], tspan)
sol = solve(prob, Tsit5())
sol_nosplit = solve(prob_nosplit, Tsit5())
@test 0 <= minimum(sol[x]) <= 1e-10 # the ball never went through the floor but got very close
@test -minimum(sol[y]) ≈ maximum(sol[y]) ≈ sqrt(2)  # the ball will never go further than √2 in either direction (gravity was changed to 1 to get this particular number)
@test 0 <= minimum(sol_nosplit[x]) <= 1e-10 # the ball never went through the floor but got very close
@test -minimum(sol_nosplit[y]) ≈ maximum(sol_nosplit[y]) ≈ sqrt(2)  # the ball will never go further than √2 in either direction (gravity was changed to 1 to get this particular number)

# tv = sort([LinRange(0, 5, 200); sol.t])
# plot(sol(tv)[y], sol(tv)[x], line_z=tv)
# vline!([-1.5, 1.5], l=(:black, 5), primary=false)
# hline!([0], l=(:black, 5), primary=false)

# issue https://github.com/SciML/ModelingToolkit.jl/issues/1386
# tests that it works for ODAESystem
@variables vs(t) v(t) vmeasured(t)
eq = [vs ~ sin(2pi * t)
      D(v) ~ vs - v
      D(vmeasured) ~ 0.0]
ev = [sin(20pi * t) ~ 0.0] => [vmeasured ~ v]
@named sys = ODESystem(eq, t, continuous_events = ev)
sys = structural_simplify(sys)
prob = ODEProblem(sys, zeros(2), (0.0, 5.1))
sol = solve(prob, Tsit5())
@test all(minimum((0:0.1:5) .- sol.t', dims = 2) .< 0.0001) # test that the solver stepped every 0.1s as dictated by event
@test sol([0.25])[vmeasured][] == sol([0.23])[vmeasured][] # test the hold property

##  https://github.com/SciML/ModelingToolkit.jl/issues/1528
Dₜ = D

@parameters u(t) [input = true]  # Indicate that this is a controlled input
@parameters y(t) [output = true] # Indicate that this is a measured output

function Mass(; name, m = 1.0, p = 0, v = 0)
    ps = @parameters m = m
    sts = @variables pos(t)=p vel(t)=v
    eqs = Dₜ(pos) ~ vel
    ODESystem(eqs, t, [pos, vel], ps; name)
end
function Spring(; name, k = 1e4)
    ps = @parameters k = k
    @variables x(t) = 0 # Spring deflection
    ODESystem(Equation[], t, [x], ps; name)
end
function Damper(; name, c = 10)
    ps = @parameters c = c
    @variables vel(t) = 0
    ODESystem(Equation[], t, [vel], ps; name)
end
function SpringDamper(; name, k = false, c = false)
    spring = Spring(; name = :spring, k)
    damper = Damper(; name = :damper, c)
    compose(ODESystem(Equation[], t; name),
        spring, damper)
end
connect_sd(sd, m1, m2) = [sd.spring.x ~ m1.pos - m2.pos, sd.damper.vel ~ m1.vel - m2.vel]
sd_force(sd) = -sd.spring.k * sd.spring.x - sd.damper.c * sd.damper.vel
@named mass1 = Mass(; m = 1)
@named mass2 = Mass(; m = 1)
@named sd = SpringDamper(; k = 1000, c = 10)
function Model(u, d = 0)
    eqs = [connect_sd(sd, mass1, mass2)
           Dₜ(mass1.vel) ~ (sd_force(sd) + u) / mass1.m
           Dₜ(mass2.vel) ~ (-sd_force(sd) + d) / mass2.m]
    @named _model = ODESystem(eqs, t; observed = [y ~ mass2.pos])
    @named model = compose(_model, mass1, mass2, sd)
end
model = Model(sin(30t))
sys = structural_simplify(model)
@test isempty(ModelingToolkit.continuous_events(sys))

let
    function testsol(osys, u0, p, tspan; tstops = Float64[], paramtotest = nothing,
            kwargs...)
        oprob = ODEProblem(complete(osys), u0, tspan, p; kwargs...)
        sol = solve(oprob, Tsit5(); tstops = tstops, abstol = 1e-10, reltol = 1e-10)
        @test isapprox(sol(1.0000000001)[1] - sol(0.999999999)[1], 1.0; rtol = 1e-6)
        paramtotest === nothing || (@test sol.ps[paramtotest] == 1.0)
        @test isapprox(sol(4.0)[1], 2 * exp(-2.0))
        sol
    end

    @parameters k t1 t2
    @variables A(t) B(t)

    cond1 = (t == t1)
    affect1 = [A ~ A + 1]
    cb1 = cond1 => affect1
    cond2 = (t == t2)
    affect2 = [k ~ 1.0]
    cb2 = cond2 => affect2

    ∂ₜ = D
    eqs = [∂ₜ(A) ~ -k * A]
    @named osys = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2])
    u0 = [A => 1.0]
    p = [k => 0.0, t1 => 1.0, t2 => 2.0]
    tspan = (0.0, 4.0)
    testsol(osys, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k)

    cond1a = (t == t1)
    affect1a = [A ~ A + 1, B ~ A]
    cb1a = cond1a => affect1a
    @named osys1 = ODESystem(eqs, t, [A, B], [k, t1, t2], discrete_events = [cb1a, cb2])
    u0′ = [A => 1.0, B => 0.0]
    sol = testsol(
        osys1, u0′, p, tspan; tstops = [1.0, 2.0], check_length = false, paramtotest = k)
    @test sol(1.0000001, idxs = B) == 2.0

    # same as above - but with set-time event syntax
    cb1‵ = [1.0] => affect1 # needs to be a Vector for the event to happen only once
    cb2‵ = [2.0] => affect2
    @named osys‵ = ODESystem(eqs, t, [A], [k], discrete_events = [cb1‵, cb2‵])
    testsol(osys‵, u0, p, tspan; paramtotest = k)

    # mixing discrete affects
    @named osys3 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵])
    testsol(osys3, u0, p, tspan; tstops = [1.0], paramtotest = k)

    # mixing with a func affect
    function affect!(integrator, u, p, ctx)
        integrator.ps[p.k] = 1.0
        nothing
    end
    cb2‵‵ = [2.0] => (affect!, [], [k], [k], nothing)
    @named osys4 = ODESystem(eqs, t, [A], [k, t1], discrete_events = [cb1, cb2‵‵])
    oprob4 = ODEProblem(complete(osys4), u0, tspan, p)
    testsol(osys4, u0, p, tspan; tstops = [1.0], paramtotest = k)

    # mixing with symbolic condition in the func affect
    cb2‵‵‵ = (t == t2) => (affect!, [], [k], [k], nothing)
    @named osys5 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵‵‵])
    testsol(osys5, u0, p, tspan; tstops = [1.0, 2.0])
    @named osys6 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb2‵‵‵, cb1])
    testsol(osys6, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k)

    # mix a continuous event too
    cond3 = A ~ 0.1
    affect3 = [k ~ 0.0]
    cb3 = cond3 => affect3
    @named osys7 = ODESystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵‵‵],
        continuous_events = [cb3])
    sol = testsol(osys7, u0, p, (0.0, 10.0); tstops = [1.0, 2.0])
    @test isapprox(sol(10.0)[1], 0.1; atol = 1e-10, rtol = 1e-10)
end

let
    function testsol(ssys, u0, p, tspan; tstops = Float64[], paramtotest = nothing,
            kwargs...)
        sprob = SDEProblem(complete(ssys), u0, tspan, p; kwargs...)
        sol = solve(sprob, RI5(); tstops = tstops, abstol = 1e-10, reltol = 1e-10)
        @test isapprox(sol(1.0000000001)[1] - sol(0.999999999)[1], 1.0; rtol = 1e-4)
        paramtotest === nothing || (@test sol.ps[paramtotest] == 1.0)
        @test isapprox(sol(4.0)[1], 2 * exp(-2.0), atol = 1e-4)
        sol
    end

    @parameters k t1 t2
    @variables A(t) B(t)

    cond1 = (t == t1)
    affect1 = [A ~ A + 1]
    cb1 = cond1 => affect1
    cond2 = (t == t2)
    affect2 = [k ~ 1.0]
    cb2 = cond2 => affect2

    ∂ₜ = D
    eqs = [∂ₜ(A) ~ -k * A]
    @named ssys = SDESystem(eqs, [0.0], t, [A], [k, t1, t2],
        discrete_events = [cb1, cb2])
    u0 = [A => 1.0]
    p = [k => 0.0, t1 => 1.0, t2 => 2.0]
    tspan = (0.0, 4.0)
    testsol(ssys, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k)

    cond1a = (t == t1)
    affect1a = [A ~ A + 1, B ~ A]
    cb1a = cond1a => affect1a
    @named ssys1 = SDESystem(eqs, [0.0], t, [A, B], [k, t1, t2],
        discrete_events = [cb1a, cb2])
    u0′ = [A => 1.0, B => 0.0]
    sol = testsol(
        ssys1, u0′, p, tspan; tstops = [1.0, 2.0], check_length = false, paramtotest = k)
    @test sol(1.0000001, idxs = 2) == 2.0

    # same as above - but with set-time event syntax
    cb1‵ = [1.0] => affect1 # needs to be a Vector for the event to happen only once
    cb2‵ = [2.0] => affect2
    @named ssys‵ = SDESystem(eqs, [0.0], t, [A], [k], discrete_events = [cb1‵, cb2‵])
    testsol(ssys‵, u0, p, tspan; paramtotest = k)

    # mixing discrete affects
    @named ssys3 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2],
        discrete_events = [cb1, cb2‵])
    testsol(ssys3, u0, p, tspan; tstops = [1.0], paramtotest = k)

    # mixing with a func affect
    function affect!(integrator, u, p, ctx)
        setp(integrator, p.k)(integrator, 1.0)
        nothing
    end
    cb2‵‵ = [2.0] => (affect!, [], [k], [k], nothing)
    @named ssys4 = SDESystem(eqs, [0.0], t, [A], [k, t1],
        discrete_events = [cb1, cb2‵‵])
    testsol(ssys4, u0, p, tspan; tstops = [1.0], paramtotest = k)

    # mixing with symbolic condition in the func affect
    cb2‵‵‵ = (t == t2) => (affect!, [], [k], [k], nothing)
    @named ssys5 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2],
        discrete_events = [cb1, cb2‵‵‵])
    testsol(ssys5, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k)
    @named ssys6 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2],
        discrete_events = [cb2‵‵‵, cb1])
    testsol(ssys6, u0, p, tspan; tstops = [1.0, 2.0], paramtotest = k)

    # mix a continuous event too
    cond3 = A ~ 0.1
    affect3 = [k ~ 0.0]
    cb3 = cond3 => affect3
    @named ssys7 = SDESystem(eqs, [0.0], t, [A], [k, t1, t2],
        discrete_events = [cb1, cb2‵‵‵],
        continuous_events = [cb3])
    sol = testsol(ssys7, u0, p, (0.0, 10.0); tstops = [1.0, 2.0])
    @test isapprox(sol(10.0)[1], 0.1; atol = 1e-10, rtol = 1e-10)
end

let rng = rng
    function testsol(jsys, u0, p, tspan; tstops = Float64[], paramtotest = nothing,
            N = 40000, kwargs...)
        jsys = complete(jsys)
        dprob = DiscreteProblem(jsys, u0, tspan, p)
        jprob = JumpProblem(jsys, dprob, Direct(); kwargs...)
        sol = solve(jprob, SSAStepper(); tstops = tstops)
        @test (sol(1.000000000001)[1] - sol(0.99999999999)[1]) == 1
        paramtotest === nothing || (@test sol.ps[paramtotest] == 1.0)
        @test sol(40.0)[1] == 0
        sol
    end

    @parameters k t1 t2
    @variables A(t) B(t)

    cond1 = (t == t1)
    affect1 = [A ~ A + 1]
    cb1 = cond1 => affect1
    cond2 = (t == t2)
    affect2 = [k ~ 1.0]
    cb2 = cond2 => affect2

    eqs = [MassActionJump(k, [A => 1], [A => -1])]
    @named jsys = JumpSystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2])
    u0 = [A => 1]
    p = [k => 0.0, t1 => 1.0, t2 => 2.0]
    tspan = (0.0, 40.0)
    testsol(jsys, u0, p, tspan; tstops = [1.0, 2.0], rng, paramtotest = k)

    cond1a = (t == t1)
    affect1a = [A ~ A + 1, B ~ A]
    cb1a = cond1a => affect1a
    @named jsys1 = JumpSystem(eqs, t, [A, B], [k, t1, t2], discrete_events = [cb1a, cb2])
    u0′ = [A => 1, B => 0]
    sol = testsol(jsys1, u0′, p, tspan; tstops = [1.0, 2.0],
        check_length = false, rng, paramtotest = k)
    @test sol(1.000000001, idxs = B) == 2

    # same as above - but with set-time event syntax
    cb1‵ = [1.0] => affect1 # needs to be a Vector for the event to happen only once
    cb2‵ = [2.0] => affect2
    @named jsys‵ = JumpSystem(eqs, t, [A], [k], discrete_events = [cb1‵, cb2‵])
    testsol(jsys‵, u0, [p[1]], tspan; rng, paramtotest = k)

    # mixing discrete affects
    @named jsys3 = JumpSystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵])
    testsol(jsys3, u0, p, tspan; tstops = [1.0], rng, paramtotest = k)

    # mixing with a func affect
    function affect!(integrator, u, p, ctx)
        integrator.ps[p.k] = 1.0
        reset_aggregated_jumps!(integrator)
        nothing
    end
    cb2‵‵ = [2.0] => (affect!, [], [k], [k], nothing)
    @named jsys4 = JumpSystem(eqs, t, [A], [k, t1], discrete_events = [cb1, cb2‵‵])
    testsol(jsys4, u0, p, tspan; tstops = [1.0], rng, paramtotest = k)

    # mixing with symbolic condition in the func affect
    cb2‵‵‵ = (t == t2) => (affect!, [], [k], [k], nothing)
    @named jsys5 = JumpSystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb1, cb2‵‵‵])
    testsol(jsys5, u0, p, tspan; tstops = [1.0, 2.0], rng, paramtotest = k)
    @named jsys6 = JumpSystem(eqs, t, [A], [k, t1, t2], discrete_events = [cb2‵‵‵, cb1])
    testsol(jsys6, u0, p, tspan; tstops = [1.0, 2.0], rng, paramtotest = k)
end

let
    function oscillator_ce(k = 1.0; name)
        sts = @variables x(t)=1.0 v(t)=0.0 F(t)
        ps = @parameters k=k Θ=0.5
        eqs = [D(x) ~ v, D(v) ~ -k * x + F]
        ev = [x ~ Θ] => [x ~ 1.0, v ~ 0.0]
        ODESystem(eqs, t, sts, ps, continuous_events = [ev]; name)
    end

    @named oscce = oscillator_ce()
    eqs = [oscce.F ~ 0]
    @named eqs_sys = ODESystem(eqs, t)
    @named oneosc_ce = compose(eqs_sys, oscce)
    oneosc_ce_simpl = structural_simplify(oneosc_ce)

    prob = ODEProblem(oneosc_ce_simpl, [], (0.0, 2.0), [])
    sol = solve(prob, Tsit5(), saveat = 0.1)

    @test typeof(oneosc_ce_simpl) == ODESystem
    @test sol[1, 6] < 1.0 # test whether x(t) decreases over time
    @test sol[1, 18] > 0.5 # test whether event happened
end

@testset "Additional SymbolicContinuousCallback options" begin
    # baseline affect (pos + neg + left root find)
    @variables c1(t)=1.0 c2(t)=1.0 # c1 = cos(t), c2 = cos(3t)
    eqs = [D(c1) ~ -sin(t); D(c2) ~ -3 * sin(3 * t)]
    record_crossings(i, u, _, c) = push!(c, i.t => i.u[u.v])
    cr1 = []
    cr2 = []
    evt1 = ModelingToolkit.SymbolicContinuousCallback(
        [c1 ~ 0], (record_crossings, [c1 => :v], [], [], cr1))
    evt2 = ModelingToolkit.SymbolicContinuousCallback(
        [c2 ~ 0], (record_crossings, [c2 => :v], [], [], cr2))
    @named trigsys = ODESystem(eqs, t; continuous_events = [evt1, evt2])
    trigsys_ss = structural_simplify(trigsys)
    prob = ODEProblem(trigsys_ss, [], (0.0, 2π))
    sol = solve(prob, Tsit5())
    required_crossings_c1 = [π / 2, 3 * π / 2]
    required_crossings_c2 = [π / 6, π / 2, 5 * π / 6, 7 * π / 6, 3 * π / 2, 11 * π / 6]
    @test maximum(abs.(first.(cr1) .- required_crossings_c1)) < 1e-4
    @test maximum(abs.(first.(cr2) .- required_crossings_c2)) < 1e-4
    @test sign.(cos.(required_crossings_c1 .- 1e-6)) == sign.(last.(cr1))
    @test sign.(cos.(3 * (required_crossings_c2 .- 1e-6))) == sign.(last.(cr2))

    # with neg affect (pos * neg + left root find)
    cr1p = []
    cr2p = []
    cr1n = []
    cr2n = []
    evt1 = ModelingToolkit.SymbolicContinuousCallback(
        [c1 ~ 0], (record_crossings, [c1 => :v], [], [], cr1p);
        affect_neg = (record_crossings, [c1 => :v], [], [], cr1n))
    evt2 = ModelingToolkit.SymbolicContinuousCallback(
        [c2 ~ 0], (record_crossings, [c2 => :v], [], [], cr2p);
        affect_neg = (record_crossings, [c2 => :v], [], [], cr2n))
    @named trigsys = ODESystem(eqs, t; continuous_events = [evt1, evt2])
    trigsys_ss = structural_simplify(trigsys)
    prob = ODEProblem(trigsys_ss, [], (0.0, 2π))
    sol = solve(prob, Tsit5(); dtmax = 0.01)
    c1_pc = filter((<=)(0) ∘ sin, required_crossings_c1)
    c1_nc = filter((>=)(0) ∘ sin, required_crossings_c1)
    c2_pc = filter(c -> -sin(3c) > 0, required_crossings_c2)
    c2_nc = filter(c -> -sin(3c) < 0, required_crossings_c2)
    @test maximum(abs.(c1_pc .- first.(cr1p))) < 1e-5
    @test maximum(abs.(c1_nc .- first.(cr1n))) < 1e-5
    @test maximum(abs.(c2_pc .- first.(cr2p))) < 1e-5
    @test maximum(abs.(c2_nc .- first.(cr2n))) < 1e-5
    @test sign.(cos.(c1_pc .- 1e-6)) == sign.(last.(cr1p))
    @test sign.(cos.(c1_nc .- 1e-6)) == sign.(last.(cr1n))
    @test sign.(cos.(3 * (c2_pc .- 1e-6))) == sign.(last.(cr2p))
    @test sign.(cos.(3 * (c2_nc .- 1e-6))) == sign.(last.(cr2n))

    # with nothing neg affect (pos * neg + left root find)
    cr1p = []
    cr2p = []
    evt1 = ModelingToolkit.SymbolicContinuousCallback(
        [c1 ~ 0], (record_crossings, [c1 => :v], [], [], cr1p); affect_neg = nothing)
    evt2 = ModelingToolkit.SymbolicContinuousCallback(
        [c2 ~ 0], (record_crossings, [c2 => :v], [], [], cr2p); affect_neg = nothing)
    @named trigsys = ODESystem(eqs, t; continuous_events = [evt1, evt2])
    trigsys_ss = structural_simplify(trigsys)
    prob = ODEProblem(trigsys_ss, [], (0.0, 2π))
    sol = solve(prob, Tsit5(); dtmax = 0.01)
    @test maximum(abs.(c1_pc .- first.(cr1p))) < 1e-5
    @test maximum(abs.(c2_pc .- first.(cr2p))) < 1e-5
    @test sign.(cos.(c1_pc .- 1e-6)) == sign.(last.(cr1p))
    @test sign.(cos.(3 * (c2_pc .- 1e-6))) == sign.(last.(cr2p))

    #mixed
    cr1p = []
    cr2p = []
    cr1n = []
    cr2n = []
    evt1 = ModelingToolkit.SymbolicContinuousCallback(
        [c1 ~ 0], (record_crossings, [c1 => :v], [], [], cr1p); affect_neg = nothing)
    evt2 = ModelingToolkit.SymbolicContinuousCallback(
        [c2 ~ 0], (record_crossings, [c2 => :v], [], [], cr2p);
        affect_neg = (record_crossings, [c2 => :v], [], [], cr2n))
    @named trigsys = ODESystem(eqs, t; continuous_events = [evt1, evt2])
    trigsys_ss = structural_simplify(trigsys)
    prob = ODEProblem(trigsys_ss, [], (0.0, 2π))
    sol = solve(prob, Tsit5(); dtmax = 0.01)
    c1_pc = filter((<=)(0) ∘ sin, required_crossings_c1)
    c2_pc = filter(c -> -sin(3c) > 0, required_crossings_c2)
    c2_nc = filter(c -> -sin(3c) < 0, required_crossings_c2)
    @test maximum(abs.(c1_pc .- first.(cr1p))) < 1e-5
    @test maximum(abs.(c2_pc .- first.(cr2p))) < 1e-5
    @test maximum(abs.(c2_nc .- first.(cr2n))) < 1e-5
    @test sign.(cos.(c1_pc .- 1e-6)) == sign.(last.(cr1p))
    @test sign.(cos.(3 * (c2_pc .- 1e-6))) == sign.(last.(cr2p))
    @test sign.(cos.(3 * (c2_nc .- 1e-6))) == sign.(last.(cr2n))

    # baseline affect w/ right rootfind (pos + neg + right root find)
    @variables c1(t)=1.0 c2(t)=1.0 # c1 = cos(t), c2 = cos(3t)
    cr1 = []
    cr2 = []
    evt1 = ModelingToolkit.SymbolicContinuousCallback(
        [c1 ~ 0], (record_crossings, [c1 => :v], [], [], cr1);
        rootfind = SciMLBase.RightRootFind)
    evt2 = ModelingToolkit.SymbolicContinuousCallback(
        [c2 ~ 0], (record_crossings, [c2 => :v], [], [], cr2);
        rootfind = SciMLBase.RightRootFind)
    @named trigsys = ODESystem(eqs, t; continuous_events = [evt1, evt2])
    trigsys_ss = structural_simplify(trigsys)
    prob = ODEProblem(trigsys_ss, [], (0.0, 2π))
    sol = solve(prob, Tsit5())
    required_crossings_c1 = [π / 2, 3 * π / 2]
    required_crossings_c2 = [π / 6, π / 2, 5 * π / 6, 7 * π / 6, 3 * π / 2, 11 * π / 6]
    @test maximum(abs.(first.(cr1) .- required_crossings_c1)) < 1e-4
    @test maximum(abs.(first.(cr2) .- required_crossings_c2)) < 1e-4
    @test sign.(cos.(required_crossings_c1 .+ 1e-6)) == sign.(last.(cr1))
    @test sign.(cos.(3 * (required_crossings_c2 .+ 1e-6))) == sign.(last.(cr2))

    # baseline affect w/ mixed rootfind (pos + neg + right root find)
    cr1 = []
    cr2 = []
    evt1 = ModelingToolkit.SymbolicContinuousCallback(
        [c1 ~ 0], (record_crossings, [c1 => :v], [], [], cr1);
        rootfind = SciMLBase.LeftRootFind)
    evt2 = ModelingToolkit.SymbolicContinuousCallback(
        [c2 ~ 0], (record_crossings, [c2 => :v], [], [], cr2);
        rootfind = SciMLBase.RightRootFind)
    @named trigsys = ODESystem(eqs, t; continuous_events = [evt1, evt2])
    trigsys_ss = structural_simplify(trigsys)
    prob = ODEProblem(trigsys_ss, [], (0.0, 2π))
    sol = solve(prob, Tsit5())
    @test maximum(abs.(first.(cr1) .- required_crossings_c1)) < 1e-4
    @test maximum(abs.(first.(cr2) .- required_crossings_c2)) < 1e-4
    @test sign.(cos.(required_crossings_c1 .- 1e-6)) == sign.(last.(cr1))
    @test sign.(cos.(3 * (required_crossings_c2 .+ 1e-6))) == sign.(last.(cr2))

    #flip order and ensure results are okay
    cr1 = []
    cr2 = []
    evt1 = ModelingToolkit.SymbolicContinuousCallback(
        [c1 ~ 0], (record_crossings, [c1 => :v], [], [], cr1);
        rootfind = SciMLBase.LeftRootFind)
    evt2 = ModelingToolkit.SymbolicContinuousCallback(
        [c2 ~ 0], (record_crossings, [c2 => :v], [], [], cr2);
        rootfind = SciMLBase.RightRootFind)
    @named trigsys = ODESystem(eqs, t; continuous_events = [evt2, evt1])
    trigsys_ss = structural_simplify(trigsys)
    prob = ODEProblem(trigsys_ss, [], (0.0, 2π))
    sol = solve(prob, Tsit5())
    @test maximum(abs.(first.(cr1) .- required_crossings_c1)) < 1e-4
    @test maximum(abs.(first.(cr2) .- required_crossings_c2)) < 1e-4
    @test sign.(cos.(required_crossings_c1 .- 1e-6)) == sign.(last.(cr1))
    @test sign.(cos.(3 * (required_crossings_c2 .+ 1e-6))) == sign.(last.(cr2))
end

@testset "Discrete variable timeseries" begin
    @variables x(t)
    @parameters a(t) b(t) c(t)
    cb1 = [x ~ 1.0] => [a ~ -a]
    function save_affect!(integ, u, p, ctx)
        integ.ps[p.b] = 5.0
    end
    cb2 = [x ~ 0.5] => (save_affect!, [], [b], [b], nothing)
    cb3 = 1.0 => [c ~ t]

    @mtkbuild sys = ODESystem(D(x) ~ cos(x), t, [x], [a, b, c]; continuous_events = [cb1, cb2], discrete_events = [cb3])
    prob = ODEProblem(sys, [x => 1.0], (0.0, 2pi), [a => 1.0, b => 2.0, c => 0.0])
    @test sort(canonicalize(Discrete(), prob.p)[1]) == [0.0, 1.0, 2.0]
    sol = solve(prob, Tsit5())

    @test sol[a] == [-1.0]
    @test sol[b] == [5.0, 5.0]
    @test sol[c] == [1.0, 2.0, 3.0]
end