From c0cac4f7e4b77b6376b6530ba098884cb2096090 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Mon, 4 May 2026 15:57:35 -0400 Subject: [PATCH 1/2] JumpProblemLibrary v2.0: drop stored DiscreteProblem, add build_jump_problem MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit SciMLBase v3 + Catalyst v16 broke the v1.x design of JumpProblemLibrary in two independent ways: 1. The `DiscreteProblem(::ReactionSystem, u0, tspan, p; eval_module=...)` constructor was removed. Lowering via `complete(jump_model(rs))` and then calling `DiscreteProblem(jsys, u0, tspan, p)` does not work either: SciMLBase v3 explicitly throws "A system with jumps cannot be used to construct a `SciMLBase.DiscreteProblem`. Consider a `JumpProblem` instead." (lib/ModelingToolkitBase/src/problems/compatibility.jl:125). 2. The new idiom for going from a `ReactionSystem` to a stochastic-jump simulation is `JumpProblem(rs, u0, tspan, p; aggregator=...)` directly — no separate `DiscreteProblem` step. This commit redesigns the library around the new idiom: - `JumpProblemNetwork.discrete_prob` is now `nothing` for every problem. The struct shape is preserved so consumers can still pattern-match on the field, but a v2.0 major bump signals the semantic change. - Adds `build_jump_problem(jpn, aggregator; kwargs...)` which constructs a `JumpProcesses.JumpProblem` from a `JumpProblemNetwork` with the supplied aggregator. Preferred consumer entry point in v2.0+. - Adds `JumpProcesses` as a direct dep (was transitive via Catalyst). - Bumps `Catalyst` compat from `15` to `16` and `DiffEqBase` to `6, 7` to allow the OrdinaryDiffEq v7 stack. Migration for downstream consumers (e.g. SciMLBenchmarks Jumps benchmark): ```julia # v1.x rn = jpn.network; prob = jpn.discrete_prob JumpProblem(rn, prob, Direct()) # v2.0 build_jump_problem(jpn, Direct()) # or directly: JumpProblem(jpn.network, jpn.u0, (0.0, jpn.tstop), jpn.rates; aggregator = Direct()) ``` `Aqua.test_all` passes locally on Julia 1.11 with the v7 stack (SciMLBase 3.7.1 / DiffEqBase 7.1.0 / Catalyst 16.1.1 / JumpProcesses 9.28.0). Co-Authored-By: Claude Opus 4.7 (1M context) Co-Authored-By: Chris Rackauckas --- lib/JumpProblemLibrary/Project.toml | 8 ++- .../src/JumpProblemLibrary.jl | 68 +++++++++++++++---- 2 files changed, 59 insertions(+), 17 deletions(-) diff --git a/lib/JumpProblemLibrary/Project.toml b/lib/JumpProblemLibrary/Project.toml index d2657d5..44d0d66 100644 --- a/lib/JumpProblemLibrary/Project.toml +++ b/lib/JumpProblemLibrary/Project.toml @@ -1,16 +1,18 @@ name = "JumpProblemLibrary" uuid = "faf0f6d7-8cee-47cb-b27c-1eb80cef534e" -version = "1.2.0" +version = "2.0.0" [deps] Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" +JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" [compat] Aqua = "0.8" -Catalyst = "15" -DiffEqBase = "6" +Catalyst = "16" +DiffEqBase = "6, 7" +JumpProcesses = "9" RuntimeGeneratedFunctions = "0.5" julia = "1.10" diff --git a/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl b/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl index 88a8ee9..315cd4f 100644 --- a/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl +++ b/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl @@ -1,6 +1,6 @@ module JumpProblemLibrary -using DiffEqBase, Catalyst +using DiffEqBase, Catalyst, JumpProcesses import RuntimeGeneratedFunctions RuntimeGeneratedFunctions.init(@__MODULE__) @@ -12,22 +12,62 @@ export prob_jump_dnarepressor, prob_jump_constproduct, prob_jump_nonlinrxs, # examples used in published benchmarks / comparisons prob_jump_multistate, prob_jump_twentygenes, prob_jump_dnadimer_repressor, # examples approximating diffusion by continuous time random walks - prob_jump_diffnetwork + prob_jump_diffnetwork, + # Builder helper for v2.0+ consumers + build_jump_problem """ -General structure to hold JumpProblem info. Needed since -the JumpProblem constructor requires the algorithm, so we -don't create the JumpProblem here. + JumpProblemNetwork + +General structure to hold the inputs needed to construct a `JumpProblem` from +a Catalyst `ReactionSystem`. Stores raw inputs (`network`, `rates`, `u0`, +`tstop`); use [`build_jump_problem`](@ref) to materialize the actual +`JumpProblem` with the aggregator of your choice. + +!!! note "v2.0 migration" + In v1.x the `discrete_prob` field held a pre-constructed + `DiscreteProblem` built directly from the `ReactionSystem`. SciMLBase v3 + no longer supports that constructor and prohibits `DiscreteProblem` from + jump-bearing systems entirely, so this field is now `nothing` for every + problem in the library. Consumer code that previously did + + ```julia + rn = jpn.network; prob = jpn.discrete_prob + JumpProblem(rn, prob, Direct()) + ``` + + should switch to + + ```julia + build_jump_problem(jpn, Direct()) + # or directly: + JumpProblem(jpn.network, jpn.u0, (0.0, jpn.tstop), jpn.rates; aggregator = Direct()) + ``` """ struct JumpProblemNetwork - network::Any # Catalyst network + network::Any # Catalyst network (ReactionSystem) rates::Any # vector of rate constants or nothing tstop::Any # time to end simulation u0::Any # initial values - discrete_prob::Any # initialized discrete problem + discrete_prob::Any # always `nothing` in v2.0+; see migration note above prob_data::Any # additional problem data, stored as a Dict end +""" + build_jump_problem(jpn::JumpProblemNetwork, aggregator; kwargs...) + +Construct a `JumpProcesses.JumpProblem` from a [`JumpProblemNetwork`](@ref) +using the supplied stochastic-simulation aggregator (e.g. `Direct()`, +`SortingDirect()`, `RSSA()`, `RSSACR()`, `DirectCR()`, ...). All keyword +arguments are forwarded to `JumpProblem`. +""" +function build_jump_problem(jpn::JumpProblemNetwork, aggregator; + save_positions = (true, true), kwargs...) + p = jpn.rates === nothing ? DiffEqBase.NullParameters() : jpn.rates + return JumpProblem(jpn.network, jpn.u0, (0.0, jpn.tstop), p; + aggregator, save_positions, kwargs...) +end + dna_rs = @reaction_network begin k1, DNA --> mRNA + DNA k2, mRNA --> mRNA + P @@ -46,7 +86,7 @@ rates = [ ] tf = 1000.0 u0 = [:DNA => 1, :mRNA => 0, :P => 0, :DNAR => 0] -prob = DiscreteProblem(dna_rs, u0, (0.0, tf), rates, eval_module = @__MODULE__) +prob = nothing Nsims = 8000 expected_avg = 5.92655375e+2 prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg) @@ -62,7 +102,7 @@ end rates = [:k1 => 1000.0, :k2 => 10.0] tf = 1.0 u0 = [:A => 0] -prob = DiscreteProblem(bd_rs, u0, (0.0, tf), rates, eval_module = @__MODULE__) +prob = nothing Nsims = 16000 expected_avg = t -> rates[1] / rates[2] .* (1.0 - exp.(-rates[2] * t)) prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean_at_t" => expected_avg) @@ -81,7 +121,7 @@ end rates = [:k1 => 1.0, :k2 => 2.0, :k3 => 0.5, :k4 => 0.75, :k5 => 0.25] tf = 0.01 u0 = [:A => 200, :B => 100, :C => 150] -prob = DiscreteProblem(nonlin_rs, u0, (0.0, tf), rates, eval_module = @__MODULE__) +prob = nothing Nsims = 32000 expected_avg = 84.876015624999994 prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg) @@ -107,7 +147,7 @@ u0 = [ :SP2 => 50.0, ] # Hill equations force use of floats! tf = 4000.0 -prob = DiscreteProblem(oscil_rs, u0, (0.0, tf), eval_module = @__MODULE__) +prob = nothing """ Oscillatory system, uses a mixture of jump types. """ @@ -166,7 +206,7 @@ u0 = [ :S6 => 0, :S7 => 0, :S8 => 0, :S9 => 0, ] tf = 100.0 -prob = DiscreteProblem(rs, u0, (0.0, tf), rates, eval_module = @__MODULE__) +prob = nothing """ Multistate model from Gupta and Mendes, "An Overview of Network-Based and -Free Approaches for Stochastic Simulation of Biochemical Systems", @@ -222,7 +262,7 @@ for i in 1:(2 * N) u0[findfirst(isequal(G[i]), unknowns(rs))] = (G[i] => 1) end tf = 2000.0 -prob = DiscreteProblem(rs, u0, (0.0, tf), eval_module = @__MODULE__) +prob = nothing """ Twenty-gene model from McCollum et al, @@ -248,7 +288,7 @@ rnpar = [ varlabels = ["G", "M", "P", "P2", "P2G"] u0 = [:G => 1000, :M => 0, :P => 0, :P2 => 0, :P2G => 0] tf = 4000.0 -prob = DiscreteProblem(rn, u0, (0.0, tf), rnpar, eval_module = @__MODULE__) +prob = nothing """ Negative feedback autoregulatory gene expression model. Dimer is the repressor. Taken from Marchetti, Priami and Thanh, From f9768c3009f3abe442368845329d385c619d862a Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Fri, 8 May 2026 10:04:54 -0400 Subject: [PATCH 2/2] JumpProblemLibrary: drop build_jump_problem helper, lean on Catalyst API MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The helper just forwarded to `JumpProblem(rs, u0, tspan, p; aggregator=...)`, the same one-liner consumers would write anyway. With the legacy `discrete_prob` field also gone, `JumpProblemNetwork` becomes a pure data container — and JumpProcesses + DiffEqBase fall out of the dependency list. Migration for consumers in this v2.0: using JumpProcesses JumpProblem(jpn.network, jpn.u0, (0.0, jpn.tstop), jpn.rates; aggregator = Direct()) Co-Authored-By: Chris Rackauckas --- lib/JumpProblemLibrary/Project.toml | 4 - .../src/JumpProblemLibrary.jl | 78 +++++-------------- 2 files changed, 18 insertions(+), 64 deletions(-) diff --git a/lib/JumpProblemLibrary/Project.toml b/lib/JumpProblemLibrary/Project.toml index 44d0d66..9c9e90a 100644 --- a/lib/JumpProblemLibrary/Project.toml +++ b/lib/JumpProblemLibrary/Project.toml @@ -4,15 +4,11 @@ version = "2.0.0" [deps] Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" -DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" -JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" [compat] Aqua = "0.8" Catalyst = "16" -DiffEqBase = "6, 7" -JumpProcesses = "9" RuntimeGeneratedFunctions = "0.5" julia = "1.10" diff --git a/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl b/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl index 315cd4f..f5700b5 100644 --- a/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl +++ b/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl @@ -1,6 +1,6 @@ module JumpProblemLibrary -using DiffEqBase, Catalyst, JumpProcesses +using Catalyst import RuntimeGeneratedFunctions RuntimeGeneratedFunctions.init(@__MODULE__) @@ -12,62 +12,28 @@ export prob_jump_dnarepressor, prob_jump_constproduct, prob_jump_nonlinrxs, # examples used in published benchmarks / comparisons prob_jump_multistate, prob_jump_twentygenes, prob_jump_dnadimer_repressor, # examples approximating diffusion by continuous time random walks - prob_jump_diffnetwork, - # Builder helper for v2.0+ consumers - build_jump_problem + prob_jump_diffnetwork """ JumpProblemNetwork -General structure to hold the inputs needed to construct a `JumpProblem` from -a Catalyst `ReactionSystem`. Stores raw inputs (`network`, `rates`, `u0`, -`tstop`); use [`build_jump_problem`](@ref) to materialize the actual -`JumpProblem` with the aggregator of your choice. +Container for the inputs needed to construct a `JumpProblem` from a +Catalyst `ReactionSystem`. Consumers materialize the actual `JumpProblem` +via the Catalyst API: -!!! note "v2.0 migration" - In v1.x the `discrete_prob` field held a pre-constructed - `DiscreteProblem` built directly from the `ReactionSystem`. SciMLBase v3 - no longer supports that constructor and prohibits `DiscreteProblem` from - jump-bearing systems entirely, so this field is now `nothing` for every - problem in the library. Consumer code that previously did - - ```julia - rn = jpn.network; prob = jpn.discrete_prob - JumpProblem(rn, prob, Direct()) - ``` - - should switch to - - ```julia - build_jump_problem(jpn, Direct()) - # or directly: - JumpProblem(jpn.network, jpn.u0, (0.0, jpn.tstop), jpn.rates; aggregator = Direct()) - ``` +```julia +using JumpProcesses +JumpProblem(jpn.network, jpn.u0, (0.0, jpn.tstop), jpn.rates; aggregator = Direct()) +``` """ struct JumpProblemNetwork - network::Any # Catalyst network (ReactionSystem) + network::Any # Catalyst ReactionSystem rates::Any # vector of rate constants or nothing tstop::Any # time to end simulation u0::Any # initial values - discrete_prob::Any # always `nothing` in v2.0+; see migration note above prob_data::Any # additional problem data, stored as a Dict end -""" - build_jump_problem(jpn::JumpProblemNetwork, aggregator; kwargs...) - -Construct a `JumpProcesses.JumpProblem` from a [`JumpProblemNetwork`](@ref) -using the supplied stochastic-simulation aggregator (e.g. `Direct()`, -`SortingDirect()`, `RSSA()`, `RSSACR()`, `DirectCR()`, ...). All keyword -arguments are forwarded to `JumpProblem`. -""" -function build_jump_problem(jpn::JumpProblemNetwork, aggregator; - save_positions = (true, true), kwargs...) - p = jpn.rates === nothing ? DiffEqBase.NullParameters() : jpn.rates - return JumpProblem(jpn.network, jpn.u0, (0.0, jpn.tstop), p; - aggregator, save_positions, kwargs...) -end - dna_rs = @reaction_network begin k1, DNA --> mRNA + DNA k2, mRNA --> mRNA + P @@ -86,14 +52,13 @@ rates = [ ] tf = 1000.0 u0 = [:DNA => 1, :mRNA => 0, :P => 0, :DNAR => 0] -prob = nothing Nsims = 8000 expected_avg = 5.92655375e+2 prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg) """ DNA negative feedback autoregulatory model. Protein acts as repressor. """ -prob_jump_dnarepressor = JumpProblemNetwork(dna_rs, rates, tf, u0, prob, prob_data) +prob_jump_dnarepressor = JumpProblemNetwork(dna_rs, rates, tf, u0, prob_data) bd_rs = @reaction_network begin k1, 0 --> A @@ -102,14 +67,13 @@ end rates = [:k1 => 1000.0, :k2 => 10.0] tf = 1.0 u0 = [:A => 0] -prob = nothing Nsims = 16000 expected_avg = t -> rates[1] / rates[2] .* (1.0 - exp.(-rates[2] * t)) prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean_at_t" => expected_avg) """ Simple birth-death process with constant production and degradation. """ -prob_jump_constproduct = JumpProblemNetwork(bd_rs, rates, tf, u0, prob, prob_data) +prob_jump_constproduct = JumpProblemNetwork(bd_rs, rates, tf, u0, prob_data) nonlin_rs = @reaction_network begin k1, 2A --> B @@ -121,14 +85,13 @@ end rates = [:k1 => 1.0, :k2 => 2.0, :k3 => 0.5, :k4 => 0.75, :k5 => 0.25] tf = 0.01 u0 = [:A => 200, :B => 100, :C => 150] -prob = nothing Nsims = 32000 expected_avg = 84.876015624999994 prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg) """ Example with a mix of nonlinear reactions, including third order """ -prob_jump_nonlinrxs = JumpProblemNetwork(nonlin_rs, rates, tf, u0, prob, prob_data) +prob_jump_nonlinrxs = JumpProblemNetwork(nonlin_rs, rates, tf, u0, prob_data) oscil_rs = @reaction_network begin 0.01, (X, Y, Z) --> 0 @@ -147,11 +110,10 @@ u0 = [ :SP2 => 50.0, ] # Hill equations force use of floats! tf = 4000.0 -prob = nothing """ Oscillatory system, uses a mixture of jump types. """ -prob_jump_osc_mixed_jumptypes = JumpProblemNetwork(oscil_rs, nothing, tf, u0, prob, nothing) +prob_jump_osc_mixed_jumptypes = JumpProblemNetwork(oscil_rs, nothing, tf, u0, nothing) specs_sym_to_name = Dict( :S1 => "R(a,l)", @@ -206,7 +168,6 @@ u0 = [ :S6 => 0, :S7 => 0, :S8 => 0, :S9 => 0, ] tf = 100.0 -prob = nothing """ Multistate model from Gupta and Mendes, "An Overview of Network-Based and -Free Approaches for Stochastic Simulation of Biochemical Systems", @@ -214,7 +175,7 @@ Computation 2018, 6, 9; doi:10.3390/computation6010009 Translated from supplementary data file: Models/Multi-state/fixed_multistate.xml """ prob_jump_multistate = JumpProblemNetwork( - rs, rates, tf, u0, prob, + rs, rates, tf, u0, Dict( "specs_to_sym_name" => specs_sym_to_name, "rates_sym_to_idx" => rsi, "params" => params @@ -262,14 +223,13 @@ for i in 1:(2 * N) u0[findfirst(isequal(G[i]), unknowns(rs))] = (G[i] => 1) end tf = 2000.0 -prob = nothing """ Twenty-gene model from McCollum et al, "The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior" Comp. Bio. and Chem., 30, pg. 39-49 (2006). """ -prob_jump_twentygenes = JumpProblemNetwork(rs, nothing, tf, u0, prob, nothing) +prob_jump_twentygenes = JumpProblemNetwork(rs, nothing, tf, u0, nothing) rn = @reaction_network begin c1, G --> G + M @@ -288,7 +248,6 @@ rnpar = [ varlabels = ["G", "M", "P", "P2", "P2G"] u0 = [:G => 1000, :M => 0, :P => 0, :P2 => 0, :P2G => 0] tf = 4000.0 -prob = nothing """ Negative feedback autoregulatory gene expression model. Dimer is the repressor. Taken from Marchetti, Priami and Thanh, @@ -296,7 +255,7 @@ Taken from Marchetti, Priami and Thanh, Springer (2017). """ prob_jump_dnadimer_repressor = JumpProblemNetwork( - rn, rnpar, tf, u0, prob, + rn, rnpar, tf, u0, Dict("specs_names" => varlabels) ) @@ -326,8 +285,7 @@ network given the number of lattice sites. u0 is a similar function that returns the initial condition vector. """ prob_jump_diffnetwork = JumpProblemNetwork( - getDiffNetwork, params, tf, getDiffu0, nothing, - nothing + getDiffNetwork, params, tf, getDiffu0, nothing ) end # module