From 29ffd992535ddcd5427d69905b6905719c747442 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Mon, 4 Nov 2024 14:17:22 -0500
Subject: [PATCH 01/12] add equations

---
 src/systems/jumps/jumpsystem.jl | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 9fa073b4b0..8b2e8e5ded 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -194,7 +194,8 @@ function JumpSystem(eqs, iv, unknowns, ps;
     # this and the treatment of continuous events are the only part 
     # unique to JumpSystems
     eqs = scalarize.(eqs)
-    ap = ArrayPartition(MassActionJump[], ConstantRateJump[], VariableRateJump[])
+    ap = ArrayPartition(
+        MassActionJump[], ConstantRateJump[], VariableRateJump[], Equation[])
     for eq in eqs
         if eq isa MassActionJump
             push!(ap.x[1], eq)
@@ -202,8 +203,10 @@ function JumpSystem(eqs, iv, unknowns, ps;
             push!(ap.x[2], eq)
         elseif eq isa VariableRateJump
             push!(ap.x[3], eq)
+        elseif eq isa Equation
+            push!(ap.x[4], eq)
         else
-            error("JumpSystem equations must contain MassActionJumps, ConstantRateJumps, or VariableRateJumps.")
+            error("JumpSystem equations must contain MassActionJumps, ConstantRateJumps, VariableRateJumps, or Equations.")
         end
     end
 
@@ -245,6 +248,7 @@ end
 has_massactionjumps(js::JumpSystem) = !isempty(equations(js).x[1])
 has_constantratejumps(js::JumpSystem) = !isempty(equations(js).x[2])
 has_variableratejumps(js::JumpSystem) = !isempty(equations(js).x[3])
+has_equations(js::JumpSystem) = !isempty(equations(js).x[4])
 
 function generate_rate_function(js::JumpSystem, rate)
     consts = collect_constants(rate)
@@ -390,6 +394,11 @@ function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple,
     if !iscomplete(sys)
         error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblem`")
     end
+
+    if has_equations(sys)
+        error("The passed in JumpSystem contains `Equations`, please use a problem type that supports equations such as such as ODEProblem.")
+    end
+
     _, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap;
         t = tspan === nothing ? nothing : tspan[1], use_union, tofloat = false, check_length = false)
     f = DiffEqBase.DISCRETE_INPLACE_DEFAULT

From 6a720421d8ebc628b8b9e97f868916954a66b36d Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Mon, 4 Nov 2024 15:06:43 -0500
Subject: [PATCH 02/12] hybrid jump-ODEs

---
 src/systems/jumps/jumpsystem.jl | 32 ++++++++++++++++++++++----------
 test/jumpsystem.jl              | 22 ++++++++++++++++++++++
 2 files changed, 44 insertions(+), 10 deletions(-)

diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 8b2e8e5ded..7950d2cc98 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -482,19 +482,31 @@ function DiffEqBase.ODEProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothi
         use_union = false,
         eval_expression = false,
         eval_module = @__MODULE__,
+        check_length = false,
         kwargs...)
     if !iscomplete(sys)
         error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblem`")
     end
 
-    _, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap;
-        t = tspan === nothing ? nothing : tspan[1], use_union, tofloat = false, check_length = false)
-
-    observedfun = ObservedFunctionCache(sys; eval_expression, eval_module)
-
-    f = (du, u, p, t) -> (du .= 0; nothing)
-    df = ODEFunction(f; sys, observed = observedfun)
-    ODEProblem(df, u0, tspan, p; kwargs...)
+    # forward everything to be an ODESystem but the jumps
+    if has_equations(sys)
+        osys = ODESystem(equations(sys).x[4], get_iv(sys), unknowns(sys), parameters(sys);
+            observed = observed(sys), name = nameof(sys), description = description(sys),
+            systems = get_systems(sys), defaults = defaults(sys), 
+            discrete_events = discrete_events(sys), 
+            parameter_dependencies = parameter_dependencies(sys),
+            metadata = get_metadata(sys), gui_metadata = get_gui_metadata(sys))
+        osys = complete(osys)
+        return ODEProblem(osys, u0map, tspan, parammap; check_length, kwargs...)
+    else
+        _, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap;
+            t = tspan === nothing ? nothing : tspan[1], use_union, tofloat = false, 
+            check_length = false)
+        f = (du, u, p, t) -> (du .= 0; nothing)
+        observedfun = ObservedFunctionCache(sys; eval_expression, eval_module)    
+        df = ODEFunction(f; sys, observed = observedfun)
+        return ODEProblem(df, u0, tspan, p; check_length, kwargs...)
+    end
 end
 
 """
@@ -530,8 +542,8 @@ function JumpProcesses.JumpProblem(js::JumpSystem, prob,
                             for j in eqs.x[2]]
     vrjs = VariableRateJump[assemble_vrj(js, j, unknowntoid; eval_expression, eval_module)
                             for j in eqs.x[3]]
-    ((prob isa DiscreteProblem) && !isempty(vrjs)) &&
-        error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps")
+    ((prob isa DiscreteProblem) && (!isempty(vrjs) || has_equations(js))) &&
+        error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps and/or coupled differential equations.")
     jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, majs)
 
     # dep graphs are only for constant rate jumps
diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl
index 30cb636226..76d11f6205 100644
--- a/test/jumpsystem.jl
+++ b/test/jumpsystem.jl
@@ -422,3 +422,25 @@ let
     @test issetequal(us, [x5])
     @test issetequal(ps, [p5])
 end
+
+# PDMP test
+let
+    @variables X(t) Y(t)
+    @parameters k1 k2
+    rate1 = k1 * X
+    affect1! = [X ~ X - 1]
+    rate2 = k1 
+    affect2! = [Y ~ Y + 1]
+    eqs = [D(X) ~ k2, D(Y) ~ -k2/100*Y]
+    vrj1 = VariableRateJump(rate1, affect1!)
+    vrj2 = VariableRateJump(rate2, affect2!)
+    @named jsys = JumpSystem([vrj1, vrj2, eqs[1], eqs[2]], t, [X, Y], [k1, k2])
+    jsys = complete(jsys)
+    u0 = [X => 0.0, Y => 0.0]
+    p = [k1 => 1.0, k2 => 20.0]
+    tspan = (0.0, 20.0)
+    oprob = ODEProblem(jsys, u0, tspan, p)
+    jprob = JumpProblem(jsys, oprob; rng)
+    sol = solve(jprob, Tsit5())
+
+end
\ No newline at end of file

From 922800bf24c38bc3d63fc6c66547e2734084e288 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Mon, 4 Nov 2024 15:42:56 -0500
Subject: [PATCH 03/12] add continuous events

---
 src/systems/jumps/jumpsystem.jl | 32 +++++++++++++++++++++-----------
 1 file changed, 21 insertions(+), 11 deletions(-)

diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 7950d2cc98..f046b7ed78 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -88,6 +88,11 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
     """
     connector_type::Any
     """
+    A `Vector{SymbolicContinuousCallback}` that model events.
+    The integrator will use root finding to guarantee that it steps at each zero crossing.
+    """
+    continuous_events::Vector{SymbolicContinuousCallback}
+    """
     A `Vector{SymbolicDiscreteCallback}` that models events. Symbolic
     analog to `SciMLBase.DiscreteCallback` that executes an affect when a given condition is
     true at the end of an integration step. Note, one must make sure to call
@@ -120,8 +125,7 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
 
     function JumpSystem{U}(
             tag, ap::U, iv, unknowns, ps, var_to_name, observed, name, description,
-            systems,
-            defaults, connector_type, devents, parameter_dependencies,
+            systems, defaults, connector_type, cevents, devents, parameter_dependencies,
             metadata = nothing, gui_metadata = nothing,
             complete = false, index_cache = nothing, isscheduled = false;
             checks::Union{Bool, Int} = true) where {U <: ArrayPartition}
@@ -136,8 +140,8 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
         end
         new{U}(tag, ap, iv, unknowns, ps, var_to_name,
             observed, name, description, systems, defaults,
-            connector_type, devents, parameter_dependencies, metadata, gui_metadata,
-            complete, index_cache, isscheduled)
+            connector_type, cevents, devents, parameter_dependencies, metadata, 
+            gui_metadata, complete, index_cache, isscheduled)
     end
 end
 function JumpSystem(tag, ap, iv, states, ps, var_to_name, args...; kwargs...)
@@ -210,13 +214,12 @@ function JumpSystem(eqs, iv, unknowns, ps;
         end
     end
 
-    (continuous_events === nothing) ||
-        error("JumpSystems currently only support discrete events.")
+    cont_callbacks = SymbolicContinuousCallbacks(continuous_events)
     disc_callbacks = SymbolicDiscreteCallbacks(discrete_events)
 
     JumpSystem{typeof(ap)}(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)),
         ap, iv′, us′, ps′, var_to_name, observed, name, description, systems,
-        defaults, connector_type, disc_callbacks, parameter_dependencies,
+        defaults, connector_type, cont_callbacks, disc_callbacks, parameter_dependencies,
         metadata, gui_metadata, checks = checks)
 end
 
@@ -399,6 +402,10 @@ function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple,
         error("The passed in JumpSystem contains `Equations`, please use a problem type that supports equations such as such as ODEProblem.")
     end
 
+    if !isempty(continuous_events(sys))
+        error("The passed in JumpSystem contains `ContinuousEvents`, please use a problem type that supports continuous events such as such as ODEProblem.")
+    end
+
     _, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap;
         t = tspan === nothing ? nothing : tspan[1], use_union, tofloat = false, check_length = false)
     f = DiffEqBase.DISCRETE_INPLACE_DEFAULT
@@ -488,12 +495,12 @@ function DiffEqBase.ODEProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothi
         error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblem`")
     end
 
-    # forward everything to be an ODESystem but the jumps
+    # forward everything to be an ODESystem but the jumps and discrete events
     if has_equations(sys)
         osys = ODESystem(equations(sys).x[4], get_iv(sys), unknowns(sys), parameters(sys);
             observed = observed(sys), name = nameof(sys), description = description(sys),
             systems = get_systems(sys), defaults = defaults(sys), 
-            discrete_events = discrete_events(sys), 
+            continuous_events = continuous_events(sys),
             parameter_dependencies = parameter_dependencies(sys),
             metadata = get_metadata(sys), gui_metadata = get_gui_metadata(sys))
         osys = complete(osys)
@@ -542,8 +549,11 @@ function JumpProcesses.JumpProblem(js::JumpSystem, prob,
                             for j in eqs.x[2]]
     vrjs = VariableRateJump[assemble_vrj(js, j, unknowntoid; eval_expression, eval_module)
                             for j in eqs.x[3]]
-    ((prob isa DiscreteProblem) && (!isempty(vrjs) || has_equations(js))) &&
-        error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps and/or coupled differential equations.")
+    if prob isa DiscreteProblem
+        if (!isempty(vrjs) || has_equations(js) || !isempty(continous_events(js)))
+            error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps, coupled differential equations, or continuous events.")
+        end
+    end
     jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, majs)
 
     # dep graphs are only for constant rate jumps

From a6c5c6cf9366686ad078084c2ce92d8a8e58cde2 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Mon, 4 Nov 2024 15:45:58 -0500
Subject: [PATCH 04/12] typos

---
 src/systems/jumps/jumpsystem.jl | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index f046b7ed78..4cf82e4bec 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -398,12 +398,8 @@ function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple,
         error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblem`")
     end
 
-    if has_equations(sys)
-        error("The passed in JumpSystem contains `Equations`, please use a problem type that supports equations such as such as ODEProblem.")
-    end
-
-    if !isempty(continuous_events(sys))
-        error("The passed in JumpSystem contains `ContinuousEvents`, please use a problem type that supports continuous events such as such as ODEProblem.")
+    if has_equations(sys) || (!isempty(continuous_events(sys)))
+        error("The passed in JumpSystem contains `Equation`s or continuous events, please use a problem type that supports these features, such as ODEProblem.")
     end
 
     _, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap;

From 4f1cd262144f152d67fc9ebbaefeee9d59c7751d Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Mon, 4 Nov 2024 17:04:07 -0500
Subject: [PATCH 05/12] spelling

---
 src/systems/jumps/jumpsystem.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 4cf82e4bec..94cf2b67e8 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -546,7 +546,7 @@ function JumpProcesses.JumpProblem(js::JumpSystem, prob,
     vrjs = VariableRateJump[assemble_vrj(js, j, unknowntoid; eval_expression, eval_module)
                             for j in eqs.x[3]]
     if prob isa DiscreteProblem
-        if (!isempty(vrjs) || has_equations(js) || !isempty(continous_events(js)))
+        if (!isempty(vrjs) || has_equations(js) || !isempty(continuous_events(js)))
             error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps, coupled differential equations, or continuous events.")
         end
     end

From 022af77d450f68fd611885afb77cb8de5820f065 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 5 Nov 2024 14:42:41 -0500
Subject: [PATCH 06/12] Add continuous event support

---
 src/systems/callbacks.jl        | 16 ++++++++--------
 src/systems/jumps/jumpsystem.jl |  6 ++----
 2 files changed, 10 insertions(+), 12 deletions(-)

diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl
index c427ae02a1..9ca69bf151 100644
--- a/src/systems/callbacks.jl
+++ b/src/systems/callbacks.jl
@@ -1,6 +1,6 @@
 #################################### system operations #####################################
 get_continuous_events(sys::AbstractSystem) = SymbolicContinuousCallback[]
-get_continuous_events(sys::AbstractODESystem) = getfield(sys, :continuous_events)
+get_continuous_events(sys::AbstractTimeDependentSystem) = getfield(sys, :continuous_events)
 has_continuous_events(sys::AbstractSystem) = isdefined(sys, :continuous_events)
 
 has_discrete_events(sys::AbstractSystem) = isdefined(sys, :discrete_events)
@@ -560,7 +560,7 @@ function compile_affect(eqs::Vector{Equation}, cb, sys, dvs, ps; outputidxs = no
     end
 end
 
-function generate_rootfinding_callback(sys::AbstractODESystem, dvs = unknowns(sys),
+function generate_rootfinding_callback(sys::AbstractTimeDependentSystem, dvs = unknowns(sys),
         ps = parameters(sys); kwargs...)
     cbs = continuous_events(sys)
     isempty(cbs) && return nothing
@@ -571,7 +571,7 @@ Generate a single rootfinding callback; this happens if there is only one equati
 generate_rootfinding_callback and thus we can produce a ContinuousCallback instead of a VectorContinuousCallback.
 """
 function generate_single_rootfinding_callback(
-        eq, cb, sys::AbstractODESystem, dvs = unknowns(sys),
+        eq, cb, sys::AbstractTimeDependentSystem, dvs = unknowns(sys),
         ps = parameters(sys); kwargs...)
     if !isequal(eq.lhs, 0)
         eq = 0 ~ eq.lhs - eq.rhs
@@ -609,7 +609,7 @@ function generate_single_rootfinding_callback(
 end
 
 function generate_vector_rootfinding_callback(
-        cbs, sys::AbstractODESystem, dvs = unknowns(sys),
+        cbs, sys::AbstractTimeDependentSystem, dvs = unknowns(sys),
         ps = parameters(sys); rootfind = SciMLBase.RightRootFind,
         reinitialization = SciMLBase.CheckInit(), kwargs...)
     eqs = map(cb -> flatten_equations(cb.eqs), cbs)
@@ -683,7 +683,7 @@ end
 """
 Compile a single continuous callback affect function(s).
 """
-function compile_affect_fn(cb, sys::AbstractODESystem, dvs, ps, kwargs)
+function compile_affect_fn(cb, sys::AbstractTimeDependentSystem, dvs, ps, kwargs)
     eq_aff = affects(cb)
     eq_neg_aff = affect_negs(cb)
     affect = compile_affect(eq_aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...)
@@ -698,7 +698,7 @@ function compile_affect_fn(cb, sys::AbstractODESystem, dvs, ps, kwargs)
     (affect = affect, affect_neg = affect_neg)
 end
 
-function generate_rootfinding_callback(cbs, sys::AbstractODESystem, dvs = unknowns(sys),
+function generate_rootfinding_callback(cbs, sys::AbstractTimeDependentSystem, dvs = unknowns(sys),
         ps = parameters(sys); kwargs...)
     eqs = map(cb -> flatten_equations(cb.eqs), cbs)
     num_eqs = length.(eqs)
@@ -859,12 +859,12 @@ merge_cb(x, ::Nothing) = x
 merge_cb(x, y) = CallbackSet(x, y)
 
 function process_events(sys; callback = nothing, kwargs...)
-    if has_continuous_events(sys)
+    if has_continuous_events(sys) && !isempty(continuous_events(sys))
         contin_cb = generate_rootfinding_callback(sys; kwargs...)
     else
         contin_cb = nothing
     end
-    if has_discrete_events(sys)
+    if has_discrete_events(sys) && !isempty(discrete_events(sys))
         discrete_cb = generate_discrete_callbacks(sys; kwargs...)
     else
         discrete_cb = nothing
diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 94cf2b67e8..e497d86141 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -485,7 +485,6 @@ function DiffEqBase.ODEProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothi
         use_union = false,
         eval_expression = false,
         eval_module = @__MODULE__,
-        check_length = false,
         kwargs...)
     if !iscomplete(sys)
         error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblem`")
@@ -496,11 +495,10 @@ function DiffEqBase.ODEProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothi
         osys = ODESystem(equations(sys).x[4], get_iv(sys), unknowns(sys), parameters(sys);
             observed = observed(sys), name = nameof(sys), description = description(sys),
             systems = get_systems(sys), defaults = defaults(sys), 
-            continuous_events = continuous_events(sys),
             parameter_dependencies = parameter_dependencies(sys),
             metadata = get_metadata(sys), gui_metadata = get_gui_metadata(sys))
         osys = complete(osys)
-        return ODEProblem(osys, u0map, tspan, parammap; check_length, kwargs...)
+        return ODEProblem(osys, u0map, tspan, parammap; check_length = false, kwargs...)
     else
         _, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap;
             t = tspan === nothing ? nothing : tspan[1], use_union, tofloat = false, 
@@ -508,7 +506,7 @@ function DiffEqBase.ODEProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothi
         f = (du, u, p, t) -> (du .= 0; nothing)
         observedfun = ObservedFunctionCache(sys; eval_expression, eval_module)    
         df = ODEFunction(f; sys, observed = observedfun)
-        return ODEProblem(df, u0, tspan, p; check_length, kwargs...)
+        return ODEProblem(df, u0, tspan, p; kwargs...)
     end
 end
 

From feac63a6c123c1658b35fd04c76ac02cbffad0dc Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 5 Nov 2024 17:27:11 -0500
Subject: [PATCH 07/12] update get_continuous_events to only apply to systems
 that have them

---
 src/systems/callbacks.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl
index 9ca69bf151..9dd78af65e 100644
--- a/src/systems/callbacks.jl
+++ b/src/systems/callbacks.jl
@@ -1,7 +1,9 @@
 #################################### system operations #####################################
-get_continuous_events(sys::AbstractSystem) = SymbolicContinuousCallback[]
-get_continuous_events(sys::AbstractTimeDependentSystem) = getfield(sys, :continuous_events)
 has_continuous_events(sys::AbstractSystem) = isdefined(sys, :continuous_events)
+function get_continuous_events(sys::AbstractSystem)
+    has_continuous_events(sys) || return SymbolicContinuousCallback[]
+    getfield(sys, :discrete_events)
+end
 
 has_discrete_events(sys::AbstractSystem) = isdefined(sys, :discrete_events)
 function get_discrete_events(sys::AbstractSystem)

From e4ed1367908e3a40b8a344c4612b1517b1f948ef Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 5 Nov 2024 17:27:36 -0500
Subject: [PATCH 08/12] typo

---
 src/systems/callbacks.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl
index 9dd78af65e..0b6b69e272 100644
--- a/src/systems/callbacks.jl
+++ b/src/systems/callbacks.jl
@@ -2,7 +2,7 @@
 has_continuous_events(sys::AbstractSystem) = isdefined(sys, :continuous_events)
 function get_continuous_events(sys::AbstractSystem)
     has_continuous_events(sys) || return SymbolicContinuousCallback[]
-    getfield(sys, :discrete_events)
+    getfield(sys, :continuous_events)
 end
 
 has_discrete_events(sys::AbstractSystem) = isdefined(sys, :discrete_events)

From 32aeac69c6d083602777ee6fd331f115d8e3d7de Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Thu, 7 Nov 2024 07:00:01 -0500
Subject: [PATCH 09/12] update master

---
 test/jumpsystem.jl | 34 ++++++++++++++++++++++++----------
 1 file changed, 24 insertions(+), 10 deletions(-)

diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl
index 76d11f6205..8a3c2adccf 100644
--- a/test/jumpsystem.jl
+++ b/test/jumpsystem.jl
@@ -301,7 +301,7 @@ end
 # basic VariableRateJump test
 let
     N = 1000  # number of simulations for testing solve accuracy
-    Random.seed!(rng, 1111)
+    Random.seed!(rng, 1111) 
     @variables A(t) B(t) C(t)
     @parameters k
     vrj = VariableRateJump(k * (sin(t) + 1), [A ~ A + 1, C ~ C + 2])
@@ -427,20 +427,34 @@ end
 let
     @variables X(t) Y(t)
     @parameters k1 k2
-    rate1 = k1 * X
-    affect1! = [X ~ X - 1]
-    rate2 = k1 
-    affect2! = [Y ~ Y + 1]
+    vrj1 = VariableRateJump(k1 * X, [X ~ X - 1])
+    vrj2 = VariableRateJump(k2, [Y ~ Y + 1])
     eqs = [D(X) ~ k2, D(Y) ~ -k2/100*Y]
-    vrj1 = VariableRateJump(rate1, affect1!)
-    vrj2 = VariableRateJump(rate2, affect2!)
     @named jsys = JumpSystem([vrj1, vrj2, eqs[1], eqs[2]], t, [X, Y], [k1, k2])
     jsys = complete(jsys)
-    u0 = [X => 0.0, Y => 0.0]
-    p = [k1 => 1.0, k2 => 20.0]
+    X0 = 0.0; Y0 = 0.0
+    u0 = [X => X0, Y => Y0]
+    k1val = 1.0; k2val = 20.0
+    p = [k1 => k1val, k2 => k2val]
     tspan = (0.0, 20.0)
     oprob = ODEProblem(jsys, u0, tspan, p)
     jprob = JumpProblem(jsys, oprob; rng)
-    sol = solve(jprob, Tsit5())
+
+    times = range(0.0, 20.0, length = 100)
+    Nsims = 2000
+    X = zeros(length(times))
+    Y = similar(X)
+    for n in 1:Nsims
+        sol = solve(jprob, Tsit5())
+        X .+= Array(sol(times; idxs = X))
+        Y .+= Array(sol(times; idxs = Y))
+    end
+    X ./= Nsims; Y ./= Nsims;
+
+    Xact(t) = X0 * exp(-k1val * t) + (k2val / k1val) * (1 - exp(-k1val * t))
+    Yact(t) = Y0 * exp(-k2val/100 * t) + (k1 / (k2val/100)) * (1 - exp(-k2val/100 * t))
+
+    @test all(abs.(X .- Xact.(times)) .<= 0.05 .* X)
+    @test all(abs.(Y .- Yact.(times)) .<= 0.05 .* Y)
 
 end
\ No newline at end of file

From 8eab77a1377052dd5c49a783d2da4b6cf7110e59 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Thu, 7 Nov 2024 14:46:58 -0500
Subject: [PATCH 10/12] test passes

---
 src/systems/jumps/jumpsystem.jl |  2 +-
 test/jumpsystem.jl              | 30 ++++++++++++++----------------
 2 files changed, 15 insertions(+), 17 deletions(-)

diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index e497d86141..e745d0c627 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -288,7 +288,7 @@ function assemble_vrj(
     outputidxs = [unknowntoid[var] for var in outputvars]
     affect = eval_or_rgf(generate_affect_function(js, vrj.affect!, outputidxs);
         eval_expression, eval_module)
-    VariableRateJump(rate, affect)
+    VariableRateJump(rate, affect; save_positions = vrj.save_positions)
 end
 
 function assemble_vrj_expr(js, vrj, unknowntoid)
diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl
index 8a3c2adccf..75bac26011 100644
--- a/test/jumpsystem.jl
+++ b/test/jumpsystem.jl
@@ -427,8 +427,8 @@ end
 let
     @variables X(t) Y(t)
     @parameters k1 k2
-    vrj1 = VariableRateJump(k1 * X, [X ~ X - 1])
-    vrj2 = VariableRateJump(k2, [Y ~ Y + 1])
+    vrj1 = VariableRateJump(k1 * X, [X ~ X - 1]; save_positions = (false, false))
+    vrj2 = VariableRateJump(k1, [Y ~ Y + 1]; save_positions = (false, false))
     eqs = [D(X) ~ k2, D(Y) ~ -k2/100*Y]
     @named jsys = JumpSystem([vrj1, vrj2, eqs[1], eqs[2]], t, [X, Y], [k1, k2])
     jsys = complete(jsys)
@@ -436,25 +436,23 @@ let
     u0 = [X => X0, Y => Y0]
     k1val = 1.0; k2val = 20.0
     p = [k1 => k1val, k2 => k2val]
-    tspan = (0.0, 20.0)
+    tspan = (0.0, 10.0)
     oprob = ODEProblem(jsys, u0, tspan, p)
-    jprob = JumpProblem(jsys, oprob; rng)
+    jprob = JumpProblem(jsys, oprob; rng, save_positions = (false, false))
 
-    times = range(0.0, 20.0, length = 100)
+    times = range(0.0, tspan[2], length = 100)
     Nsims = 2000
-    X = zeros(length(times))
-    Y = similar(X)
+    Xv = zeros(length(times))
+    Yv = copy(Xv)
     for n in 1:Nsims
-        sol = solve(jprob, Tsit5())
-        X .+= Array(sol(times; idxs = X))
-        Y .+= Array(sol(times; idxs = Y))
+        sol = solve(jprob, Tsit5(); saveat = times)
+        Xv .+= sol[1,:] #sol(times; idxs = X)
+        Yv .+= sol[2,:] #sol(times; idxs = Y)
     end
-    X ./= Nsims; Y ./= Nsims;
+    Xv ./= Nsims; Yv ./= Nsims;
 
     Xact(t) = X0 * exp(-k1val * t) + (k2val / k1val) * (1 - exp(-k1val * t))
-    Yact(t) = Y0 * exp(-k2val/100 * t) + (k1 / (k2val/100)) * (1 - exp(-k2val/100 * t))
-
-    @test all(abs.(X .- Xact.(times)) .<= 0.05 .* X)
-    @test all(abs.(Y .- Yact.(times)) .<= 0.05 .* Y)
-
+    Yact(t) = Y0 * exp(-k2val/100 * t) + (k1val / (k2val/100)) * (1 - exp(-k2val/100 * t))
+    @test all(abs.(Xv .- Xact.(times)) .<= 0.05 .* Xv)
+    @test all(abs.(Yv .- Yact.(times)) .<= 0.05 .* Yv)
 end
\ No newline at end of file

From 79344744e5a4e38771f49509245356b94e4a9ecf Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Thu, 7 Nov 2024 14:47:08 -0500
Subject: [PATCH 11/12] updates

---
 test/jumpsystem.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl
index 75bac26011..7f049b0aa9 100644
--- a/test/jumpsystem.jl
+++ b/test/jumpsystem.jl
@@ -446,8 +446,10 @@ let
     Yv = copy(Xv)
     for n in 1:Nsims
         sol = solve(jprob, Tsit5(); saveat = times)
-        Xv .+= sol[1,:] #sol(times; idxs = X)
-        Yv .+= sol[2,:] #sol(times; idxs = Y)
+
+        # use direct indexing as much faster than symbolic indexing
+        Xv .+= sol[1,:] 
+        Yv .+= sol[2,:]
     end
     Xv ./= Nsims; Yv ./= Nsims;
 

From 0a223169fa74f22b033c90e88e4b64874b46aaa4 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Thu, 7 Nov 2024 17:50:13 -0500
Subject: [PATCH 12/12] more sampled

---
 test/jumpsystem.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl
index 7f049b0aa9..204f4f2fb5 100644
--- a/test/jumpsystem.jl
+++ b/test/jumpsystem.jl
@@ -441,13 +441,11 @@ let
     jprob = JumpProblem(jsys, oprob; rng, save_positions = (false, false))
 
     times = range(0.0, tspan[2], length = 100)
-    Nsims = 2000
+    Nsims = 4000
     Xv = zeros(length(times))
     Yv = copy(Xv)
     for n in 1:Nsims
         sol = solve(jprob, Tsit5(); saveat = times)
-
-        # use direct indexing as much faster than symbolic indexing
         Xv .+= sol[1,:] 
         Yv .+= sol[2,:]
     end