From 6dde4d6fbc1c61d011d744eda2256d0ca420323f Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 13 Sep 2024 17:13:37 -0400 Subject: [PATCH 01/14] try dropping scalarizing --- src/network_analysis.jl | 7 +++---- src/reactionsystem.jl | 8 ++++---- test/reactionsystem_core/reactionsystem.jl | 8 +++----- 3 files changed, 10 insertions(+), 13 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 8fb8f2a208..2806a0b4b6 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -657,16 +657,15 @@ function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_o indepspecs = sts[indepidxs] depidxs = col_order[(r + 1):end] depspecs = sts[depidxs] - constants = MT.unwrap.(MT.scalarize(only( - @parameters $(CONSERVED_CONSTANT_SYMBOL)[1:nullity] [conserved = true]))) + constants = MT.unwrap(only( + @parameters $(CONSERVED_CONSTANT_SYMBOL)[1:nullity] [conserved = true])) conservedeqs = Equation[] constantdefs = Equation[] for (i, depidx) in enumerate(depidxs) scaleby = (N[i, depidx] != 1) ? N[i, depidx] : one(eltype(N)) (scaleby != 0) || error("Error, found a zero in the conservation law matrix where " - * - "one was not expected.") + * "one was not expected.") coefs = @view N[i, indepidxs] terms = sum(p -> p[1] / scaleby * p[2], zip(coefs, indepspecs)) eq = depspecs[i] ~ constants[i] - terms diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 7497544c24..fce066e3b4 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -405,11 +405,11 @@ function ReactionSystem(eqs, iv, unknowns, ps; sivs′ = if spatial_ivs === nothing Vector{typeof(iv′)}() else - value.(MT.scalarize(spatial_ivs)) + value.(spatial_ivs) end - unknowns′ = sort!(value.(MT.scalarize(unknowns)), by = !isspecies) + unknowns′ = sort!(value.(unknowns), by = !isspecies) spcs = filter(isspecies, unknowns′) - ps′ = value.(MT.scalarize(ps)) + ps′ = value.(ps) # Checks that no (by Catalyst) forbidden symbols are used. allsyms = Iterators.flatten((ps′, unknowns′)) @@ -494,7 +494,7 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; t = value(iv) ivs = Set([t]) if (spatial_ivs !== nothing) - for siv in (MT.scalarize(spatial_ivs)) + for siv in (spatial_ivs) push!(ivs, value(siv)) end end diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl index 685b373ffe..120ff3be91 100644 --- a/test/reactionsystem_core/reactionsystem.jl +++ b/test/reactionsystem_core/reactionsystem.jl @@ -868,11 +868,9 @@ end let @species (A(t))[1:20] using ModelingToolkit: value - @test isspecies(value(A)) - @test isspecies(value(A[2])) - Av = value.(ModelingToolkit.scalarize(A)) - @test isspecies(Av[2]) - @test isequal(value(Av[2]), value(A[2])) + Av = value(A) + @test isspecies(Av) + @test all(i -> isspecies(Av[i]), 1:length(Av)) end # Test mixed models are formulated correctly. From 3776534e74dc3023353ac7fe86efb60328f4f91f Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sat, 14 Sep 2024 09:34:47 -0400 Subject: [PATCH 02/14] don't assume things are Nums --- src/reactionsystem.jl | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index fce066e3b4..8c69214522 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -467,7 +467,7 @@ end # Two-argument constructor (reactions/equations and time variable). # Calls the `make_ReactionSystem_internal`, which in turn calls the four-argument constructor. function ReactionSystem(rxs::Vector, iv = Catalyst.DEFAULT_IV; kwargs...) - make_ReactionSystem_internal(rxs, iv, Vector{Num}(), Vector{Num}(); kwargs...) + make_ReactionSystem_internal(rxs, iv, [], []; kwargs...) end # One-argument constructor. Creates an emtoy `ReactionSystem` from a time independent variable only. @@ -503,6 +503,7 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; # Preallocates the `vars` set, which is used by `findvars!` us = OrderedSet{eltype(us_in)}(us_in) ps = OrderedSet{eltype(ps_in)}(ps_in) + @show eltype(ps_in) vars = OrderedSet() # Extracts the reactions and equations from the combined reactions + equations input vector. @@ -548,7 +549,24 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; # Converts the found unknowns and parameters to vectors. usv = collect(us) - psv = collect(ps) + + new_ps = OrderedSet() + for p in ps + if iscall(p) && operation(p) === getindex + par = arguments(p)[begin] + if Symbolics.shape(Symbolics.unwrap(par)) !== Symbolics.Unknown() && + all(par[i] in ps for i in eachindex(par)) + @show "here" + push!(new_ps, par) + else + push!(new_ps, p) + end + else + push!(new_ps, p) + end + end + @show new_ps + psv = collect(new_ps) # Passes the processed input into the next `ReactionSystem` call. ReactionSystem(fulleqs, t, usv, psv; spatial_ivs, continuous_events, From e51718576137e46eb9b204772dd1fdf0c301d029 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sat, 14 Sep 2024 09:35:17 -0400 Subject: [PATCH 03/14] remove debug printing --- src/reactionsystem.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 8c69214522..ea90b25763 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -556,7 +556,6 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; par = arguments(p)[begin] if Symbolics.shape(Symbolics.unwrap(par)) !== Symbolics.Unknown() && all(par[i] in ps for i in eachindex(par)) - @show "here" push!(new_ps, par) else push!(new_ps, p) From d769f5ab1e99bb79e38c7b292cb185e5b548e57b Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sat, 14 Sep 2024 09:35:41 -0400 Subject: [PATCH 04/14] more debug prints --- src/reactionsystem.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index ea90b25763..e502bed138 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -564,7 +564,6 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; push!(new_ps, p) end end - @show new_ps psv = collect(new_ps) # Passes the processed input into the next `ReactionSystem` call. From 745f449410e81ddc06b5853e59283ea8f453d97f Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sat, 14 Sep 2024 10:28:13 -0400 Subject: [PATCH 05/14] start fixing tests --- src/reactionsystem.jl | 1 - test/reactionsystem_core/reactionsystem.jl | 16 +++++++--------- 2 files changed, 7 insertions(+), 10 deletions(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index e502bed138..15afb2429f 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -503,7 +503,6 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; # Preallocates the `vars` set, which is used by `findvars!` us = OrderedSet{eltype(us_in)}(us_in) ps = OrderedSet{eltype(ps_in)}(ps_in) - @show eltype(ps_in) vars = OrderedSet() # Extracts the reactions and equations from the combined reactions + equations input vector. diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl index 120ff3be91..f282528b91 100644 --- a/test/reactionsystem_core/reactionsystem.jl +++ b/test/reactionsystem_core/reactionsystem.jl @@ -40,7 +40,7 @@ rxs = [Reaction(k[1], nothing, [A]), # 0 -> A Reaction(k[19] * t, [A], [B]), # A -> B with non constant rate. Reaction(k[20] * t * A, [B, C], [D], [2, 1], [2]), # 2A +B -> 2C with non constant rate. ] -@named rs = ReactionSystem(rxs, t, [A, B, C, D], k) +@named rs = ReactionSystem(rxs, t, [A, B, C, D], [k]) rs = complete(rs) odesys = complete(convert(ODESystem, rs)) sdesys = complete(convert(SDESystem, rs)) @@ -109,11 +109,12 @@ end # Defaults test. let - def_p = [ki => float(i) for (i, ki) in enumerate(k)] + kvals = Float64.(1:length(k)) + def_p = [k => kvals] def_u0 = [A => 0.5, B => 1.0, C => 1.5, D => 2.0] defs = merge(Dict(def_p), Dict(def_u0)) - @named rs = ReactionSystem(rxs, t, [A, B, C, D], k; defaults = defs) + @named rs = ReactionSystem(rxs, t, [A, B, C, D], [k]; defaults = defs) rs = complete(rs) odesys = complete(convert(ODESystem, rs)) sdesys = complete(convert(SDESystem, rs)) @@ -126,15 +127,11 @@ let defs u0map = [A => 5.0] - pmap = [k[1] => 5.0] + kvals[1] = 5.0 + pmap = [k => kvals] prob = ODEProblem(rs, u0map, (0, 10.0), pmap) @test prob.ps[k[1]] == 5.0 @test prob.u0[1] == 5.0 - u0 = [10.0, 11.0, 12.0, 13.0] - ps = [float(x) for x in 100:119] - prob = ODEProblem(rs, u0, (0, 10.0), ps) - @test [prob.ps[k[i]] for i in 1:20] == ps - @test prob.u0 == u0 end ### Check ODE, SDE, and Jump Functions ### @@ -144,6 +141,7 @@ end let u = rnd_u0(rs, rng) p = rnd_ps(rs, rng) + @show u,p du = oderhs(last.(u), last.(p), 0.0) G = sdenoise(last.(u), last.(p), 0.0) sdesys = complete(convert(SDESystem, rs)) From 871d3da60725c6446e04cb5cc0a3400f0b37cbb4 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sat, 14 Sep 2024 10:33:42 -0400 Subject: [PATCH 06/14] remove shows --- test/reactionsystem_core/reactionsystem.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl index f282528b91..7553982658 100644 --- a/test/reactionsystem_core/reactionsystem.jl +++ b/test/reactionsystem_core/reactionsystem.jl @@ -141,7 +141,6 @@ end let u = rnd_u0(rs, rng) p = rnd_ps(rs, rng) - @show u,p du = oderhs(last.(u), last.(p), 0.0) G = sdenoise(last.(u), last.(p), 0.0) sdesys = complete(convert(SDESystem, rs)) From ad58d367afc32570a8be004f446d1fcf9a098c23 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 20 Sep 2024 10:50:30 -0400 Subject: [PATCH 07/14] more test fixes --- test/reactionsystem_core/reactionsystem.jl | 10 ++++++---- test/test_functions.jl | 8 ++++---- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl index 7553982658..97d32eb99f 100644 --- a/test/reactionsystem_core/reactionsystem.jl +++ b/test/reactionsystem_core/reactionsystem.jl @@ -46,11 +46,12 @@ odesys = complete(convert(ODESystem, rs)) sdesys = complete(convert(SDESystem, rs)) # Hard coded ODE rhs. -function oderhs(u, k, t) +function oderhs(u, kv, t) A = u[1] B = u[2] C = u[3] D = u[4] + k = kv[1] du = zeros(eltype(u), 4) du[1] = k[1] - k[3] * A + k[4] * C + 2 * k[5] * C - k[6] * A * B + k[7] * B^2 / 2 - k[9] * A * B - k[10] * A^2 - k[11] * A^2 / 2 - k[12] * A * B^3 * C^4 / 144 - @@ -68,11 +69,12 @@ function oderhs(u, k, t) end # SDE noise coefs. -function sdenoise(u, k, t) +function sdenoise(u, kv, t) A = u[1] B = u[2] C = u[3] D = u[4] + k = kv[1] G = zeros(eltype(u), length(k), length(u)) z = zero(eltype(u)) @@ -178,7 +180,7 @@ let Reaction(k[19] * t, [D], [E]), # D -> E with non constant rate. Reaction(k[20] * t * A, [D, E], [F], [2, 1], [2]), # 2D + E -> 2F with non constant rate. ] - @named rs = ReactionSystem(rxs, t, [A, B, C, D, E, F], k) + @named rs = ReactionSystem(rxs, t, [A, B, C, D, E, F], [k]) rs = complete(rs) js = complete(convert(JumpSystem, rs)) @@ -190,7 +192,7 @@ let @test all(map(i -> typeof(equations(js)[i]) <: JumpProcesses.VariableRateJump, vidxs)) p = rand(rng, length(k)) - pmap = parameters(js) .=> p + pmap = [k => p] u0 = rand(rng, 2:10, 6) u0map = unknowns(js) .=> u0 ttt = rand(rng) diff --git a/test/test_functions.jl b/test/test_functions.jl index ca35b3c5c8..5eb3d2c7d5 100644 --- a/test/test_functions.jl +++ b/test/test_functions.jl @@ -7,22 +7,22 @@ using Random, Test # Generates a random initial condition (in the form of a map). Each value is a Float64. function rnd_u0(sys, rng; factor = 1.0, min = 0.0) - return [u => min + factor * rand(rng) for u in unknowns(sys)] + return [u => (min .+ factor .* rand(rng, size(u)...)) for u in unknowns(sys)] end # Generates a random initial condition (in the form of a map). Each value is a Int64. function rnd_u0_Int64(sys, rng; n = 5, min = 0) - return [u => min + rand(rng, 1:n) for u in unknowns(sys)] + return [u => (min .+ rand(rng, 1:n, size(u)...)) for u in unknowns(sys)] end # Generates a random parameter set (in the form of a map). Each value is a Float64. function rnd_ps(sys, rng; factor = 1.0, min = 0.0) - return [p => min + factor * rand(rng) for p in parameters(sys)] + return [p => ( min .+ factor .* rand(rng, size(p)...)) for p in parameters(sys)] end # Generates a random parameter set (in the form of a map). Each value is a Float64. function rnd_ps_Int64(sys, rng; n = 5, min = 0) - return [p => min + rand(rng, 1:n) for p in parameters(sys)] + return [p => (min .+ rand(rng, 1:n, size(p)...)) for p in parameters(sys)] end # Used to convert a generated initial condition/parameter set to a vector that can be used for normal From f7e68f75d6447c38502820247cc18e4f572c0406 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 20 Sep 2024 12:14:39 -0400 Subject: [PATCH 08/14] disable broken tests --- src/dsl.jl | 2 ++ src/reactionsystem.jl | 5 +++-- test/dsl/dsl_options.jl | 9 +++------ 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/dsl.jl b/src/dsl.jl index bb3fbd6df8..4954e367e4 100644 --- a/src/dsl.jl +++ b/src/dsl.jl @@ -382,6 +382,8 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) push!(rxexprs.args, equation) end + #println(observed_vars) + # Output code corresponding to the reaction system. quote $ivexpr diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 15afb2429f..0b63d35034 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -486,8 +486,9 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; observed = [], kwargs...) # Filters away any potential observables from `states` and `spcs`. - obs_vars = [obs_eq.lhs for obs_eq in observed] - us_in = filter(u -> !any(isequal(u, obs_var) for obs_var in obs_vars), us_in) + obs_vars = Set(obs_eq.lhs for obs_eq in observed) + any(in(obs_vars), us_in) && error("Found an observable in the list of unknowns. This is not allowed.") + # us_in = filter(u -> !any(isequal(u, obs_var) for obs_var in obs_vars), us_in) # Creates a combined iv vector (iv and sivs). This is used later in the function (so that # independent variables can be excluded when encountered quantities are added to `us` and `ps`). diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl index 04247d0778..8174598fd4 100644 --- a/test/dsl/dsl_options.jl +++ b/test/dsl/dsl_options.jl @@ -604,7 +604,7 @@ let k, 0 --> X1 + X2 end @test isequal(observed(rn1)[1].rhs, observed(rn2)[1].rhs) - @test isequal(observed(rn1)[1].lhs.metadata, observed(rn2)[1].lhs.metadata) + @test_broken isequal(observed(rn1)[1].lhs.metadata, observed(rn2)[1].lhs.metadata) @test isequal(unknowns(rn1), unknowns(rn2)) # Case with metadata. @@ -618,7 +618,7 @@ let k, 0 --> X1 + X2 end @test isequal(observed(rn3)[1].rhs, observed(rn4)[1].rhs) - @test isequal(observed(rn3)[1].lhs.metadata, observed(rn4)[1].lhs.metadata) + @test_broken isequal(observed(rn3)[1].lhs.metadata, observed(rn4)[1].lhs.metadata) @test isequal(unknowns(rn3), unknowns(rn4)) end @@ -822,10 +822,7 @@ let @variables X(t) @equations 2X ~ $c - X end) - - u0 = [rn.X => 0.0] - ps = [] - oprob = ODEProblem(rn, u0, (0.0, 100.0), ps; structural_simplify=true) + oprob = ODEProblem(rn, [], (0.0, 100.0); structural_simplify=true) sol = solve(oprob, Tsit5(); abstol=1e-9, reltol=1e-9) @test sol[rn.X][end] ≈ 2.0 end From bac10911119a6c39cd3b8f0f5a695b9ec55f7753 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 20 Sep 2024 12:50:00 -0400 Subject: [PATCH 09/14] remove scalarizing of parameters in macro --- src/dsl.jl | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/src/dsl.jl b/src/dsl.jl index 4954e367e4..db703383d1 100644 --- a/src/dsl.jl +++ b/src/dsl.jl @@ -369,10 +369,10 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) sexprs = get_sexpr(species_extracted, options; iv_symbols = ivs) vexprs = get_sexpr(vars_extracted, options, :variables; iv_symbols = ivs) pexprs = get_pexpr(parameters_extracted, options) - ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps") - vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars") - sps, spssym = scalarize_macro(!isempty(species), sexprs, "specs") - comps, compssym = scalarize_macro(!isempty(compound_species), compound_expr, "comps") + ps, pssym = assign_expr_to_var(!isempty(parameters), pexprs, "ps") + vars, varssym = assign_expr_to_var(!isempty(variables), vexprs, "vars", true) + sps, spssym = assign_expr_to_var(!isempty(species), sexprs, "specs", true) + comps, compssym = assign_expr_to_var(!isempty(compound_species), compound_expr, "comps", true) rxexprs = :(CatalystEqType[]) for reaction in reactions push!(rxexprs.args, get_rxexprs(reaction)) @@ -382,8 +382,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) push!(rxexprs.args, equation) end - #println(observed_vars) - # Output code corresponding to the reaction system. quote $ivexpr @@ -593,14 +591,21 @@ function get_rxexprs(rxstruct) end # takes a ModelingToolkit declaration macro like @parameters and returns an expression -# that calls the macro and then scalarizes all the symbols created into a vector of Nums -function scalarize_macro(nonempty, ex, name) +# that calls the macro and saves it in a variable named namesym. +# also scalarizes if desired +function assign_expr_to_var(nonempty, ex, name, scalarize = false) namesym = gensym(name) if nonempty - symvec = gensym() - ex = quote - $symvec = $ex - $namesym = reduce(vcat, Symbolics.scalarize($symvec)) + if scalarize + symvec = gensym() + ex = quote + $symvec = $ex + $namesym = reduce(vcat, Symbolics.scalarize($symvec)) + end + else + ex = quote + $namesym = $ex + end end else ex = :($namesym = Num[]) From 6467291766abd7ce3e5126695563322406c1bf26 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 20 Sep 2024 12:51:52 -0400 Subject: [PATCH 10/14] make scalarizing a keyword --- src/dsl.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/dsl.jl b/src/dsl.jl index db703383d1..ae73c7b71e 100644 --- a/src/dsl.jl +++ b/src/dsl.jl @@ -370,9 +370,11 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) vexprs = get_sexpr(vars_extracted, options, :variables; iv_symbols = ivs) pexprs = get_pexpr(parameters_extracted, options) ps, pssym = assign_expr_to_var(!isempty(parameters), pexprs, "ps") - vars, varssym = assign_expr_to_var(!isempty(variables), vexprs, "vars", true) - sps, spssym = assign_expr_to_var(!isempty(species), sexprs, "specs", true) - comps, compssym = assign_expr_to_var(!isempty(compound_species), compound_expr, "comps", true) + vars, varssym = assign_expr_to_var(!isempty(variables), vexprs, "vars"; + scalarize = true) + sps, spssym = assign_expr_to_var(!isempty(species), sexprs, "specs"; scalarize = true) + comps, compssym = assign_expr_to_var(!isempty(compound_species), compound_expr, + "comps"; scalarize = true) rxexprs = :(CatalystEqType[]) for reaction in reactions push!(rxexprs.args, get_rxexprs(reaction)) @@ -591,9 +593,9 @@ function get_rxexprs(rxstruct) end # takes a ModelingToolkit declaration macro like @parameters and returns an expression -# that calls the macro and saves it in a variable named namesym. -# also scalarizes if desired -function assign_expr_to_var(nonempty, ex, name, scalarize = false) +# that calls the macro and saves it in a variable given by namesym based on name. +# scalarizes if desired +function assign_expr_to_var(nonempty, ex, name; scalarize = false) namesym = gensym(name) if nonempty if scalarize From 2c904fba212df2e5bafb99f2c7eb8c2ece35b56c Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 20 Sep 2024 12:56:13 -0400 Subject: [PATCH 11/14] format --- src/dsl.jl | 6 +++--- src/network_analysis.jl | 3 ++- src/reactionsystem.jl | 3 ++- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/dsl.jl b/src/dsl.jl index ae73c7b71e..91a2fa3851 100644 --- a/src/dsl.jl +++ b/src/dsl.jl @@ -370,10 +370,10 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) vexprs = get_sexpr(vars_extracted, options, :variables; iv_symbols = ivs) pexprs = get_pexpr(parameters_extracted, options) ps, pssym = assign_expr_to_var(!isempty(parameters), pexprs, "ps") - vars, varssym = assign_expr_to_var(!isempty(variables), vexprs, "vars"; + vars, varssym = assign_expr_to_var(!isempty(variables), vexprs, "vars"; scalarize = true) - sps, spssym = assign_expr_to_var(!isempty(species), sexprs, "specs"; scalarize = true) - comps, compssym = assign_expr_to_var(!isempty(compound_species), compound_expr, + sps, spssym = assign_expr_to_var(!isempty(species), sexprs, "specs"; scalarize = true) + comps, compssym = assign_expr_to_var(!isempty(compound_species), compound_expr, "comps"; scalarize = true) rxexprs = :(CatalystEqType[]) for reaction in reactions diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 2806a0b4b6..3e42e064c1 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -665,7 +665,8 @@ function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_o for (i, depidx) in enumerate(depidxs) scaleby = (N[i, depidx] != 1) ? N[i, depidx] : one(eltype(N)) (scaleby != 0) || error("Error, found a zero in the conservation law matrix where " - * "one was not expected.") + * + "one was not expected.") coefs = @view N[i, indepidxs] terms = sum(p -> p[1] / scaleby * p[2], zip(coefs, indepspecs)) eq = depspecs[i] ~ constants[i] - terms diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 0b63d35034..3c6307ec37 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -487,7 +487,8 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; # Filters away any potential observables from `states` and `spcs`. obs_vars = Set(obs_eq.lhs for obs_eq in observed) - any(in(obs_vars), us_in) && error("Found an observable in the list of unknowns. This is not allowed.") + any(in(obs_vars), us_in) && + error("Found an observable in the list of unknowns. This is not allowed.") # us_in = filter(u -> !any(isequal(u, obs_var) for obs_var in obs_vars), us_in) # Creates a combined iv vector (iv and sivs). This is used later in the function (so that From 55d8e4810b83934f62d9488f511d2aa56c992c05 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 20 Sep 2024 12:57:42 -0400 Subject: [PATCH 12/14] remove dropped code --- src/reactionsystem.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 3c6307ec37..fc9fb69f4e 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -485,11 +485,10 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; spatial_ivs = nothing, continuous_events = [], discrete_events = [], observed = [], kwargs...) - # Filters away any potential observables from `states` and `spcs`. + # Error if any observables have been declared a species or variable obs_vars = Set(obs_eq.lhs for obs_eq in observed) any(in(obs_vars), us_in) && error("Found an observable in the list of unknowns. This is not allowed.") - # us_in = filter(u -> !any(isequal(u, obs_var) for obs_var in obs_vars), us_in) # Creates a combined iv vector (iv and sivs). This is used later in the function (so that # independent variables can be excluded when encountered quantities are added to `us` and `ps`). From 2c1a46f4e01dba4ecacf5f7e247c530c588b8aee Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 20 Sep 2024 14:12:55 -0400 Subject: [PATCH 13/14] fix another test --- test/simulation_and_solving/simulate_SDEs.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/simulation_and_solving/simulate_SDEs.jl b/test/simulation_and_solving/simulate_SDEs.jl index 278c8ac915..16f662e7a4 100644 --- a/test/simulation_and_solving/simulate_SDEs.jl +++ b/test/simulation_and_solving/simulate_SDEs.jl @@ -204,9 +204,9 @@ let (k1, k2), X1 ↔ X2, ([noise_scaling=η[1]],[noise_scaling=η[2]]) end @unpack k1, k2, η = noise_scaling_network_2 - sprob_2_1 = SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η[1] => 2.0, η[2] => 2.0]) - sprob_2_2 = SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η[1] => 2.0, η[2] => 0.2]) - sprob_2_3 = SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η[1] => 0.2, η[2] => 0.2]) + sprob_2_1 = SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η => [2.0, 2.0]]) + sprob_2_2 = SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η => [2.0, 0.2]]) + sprob_2_3 = SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η => [0.2, 0.2]]) for repeat in 1:5 sol_2_1 = solve(sprob_2_1, ImplicitEM(); saveat = 1.0, seed = rand(rng, 1:100)) sol_2_2 = solve(sprob_2_2, ImplicitEM(); saveat = 1.0, seed = rand(rng, 1:100)) From 9e2240ffc0a82da4c7457bac0fcdf0ac53d4683a Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 20 Sep 2024 16:22:28 -0400 Subject: [PATCH 14/14] fix doc issue --- docs/src/model_creation/dsl_advanced.md | 2 +- .../programmatic_generative_linear_pathway.md | 2 +- .../smoluchowski_coagulation_equation.md | 2 +- src/reactionsystem.jl | 4 ++-- test/dsl/dsl_advanced_model_construction.jl | 17 +++++++++++++++-- test/network_analysis/conservation_laws.jl | 2 +- test/reactionsystem_core/reactionsystem.jl | 2 +- test/upstream/mtk_problem_inputs.jl | 2 +- 8 files changed, 23 insertions(+), 10 deletions(-) diff --git a/docs/src/model_creation/dsl_advanced.md b/docs/src/model_creation/dsl_advanced.md index 0a443e02a2..a04a632aa9 100644 --- a/docs/src/model_creation/dsl_advanced.md +++ b/docs/src/model_creation/dsl_advanced.md @@ -256,7 +256,7 @@ Sometimes, one wishes to declare a large number of similar parameters or species using Catalyst # hide two_state_model = @reaction_network begin @parameters k[1:2] - @species X(t)[1:2] + @species (X(t))[1:2] (k[1],k[2]), X[1] <--> X[2] end ``` diff --git a/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md b/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md index e6a12b876e..eb5180c468 100644 --- a/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md +++ b/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md @@ -84,7 +84,7 @@ t = default_t() @parameters τ function generate_lp(n) # Creates a vector `X` with n+1 species. - @species X(t)[1:n+1] + @species (X(t))[1:n+1] @species Xend(t) # Generate diff --git a/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md b/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md index 6b0ab3de76..7ed5a96bff 100644 --- a/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md +++ b/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md @@ -98,7 +98,7 @@ for n = 1:nr [1, 1], [1])) end end -@named rs = ReactionSystem(rx, t, collect(X), collect(k)) +@named rs = ReactionSystem(rx, t, collect(X), [k]) rs = complete(rs) ``` We now convert the [`ReactionSystem`](@ref) into a `ModelingToolkit.JumpSystem`, and solve it using Gillespie's direct method. For details on other possible solvers (SSAs), see the [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/types/jump_types/) documentation diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index fc9fb69f4e..4e28741cd7 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -502,8 +502,8 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; # Initialises the new unknowns and parameter vectors. # Preallocates the `vars` set, which is used by `findvars!` - us = OrderedSet{eltype(us_in)}(us_in) - ps = OrderedSet{eltype(ps_in)}(ps_in) + us = OrderedSet{Any}(us_in) + ps = OrderedSet{Any}(ps_in) vars = OrderedSet() # Extracts the reactions and equations from the combined reactions + equations input vector. diff --git a/test/dsl/dsl_advanced_model_construction.jl b/test/dsl/dsl_advanced_model_construction.jl index bac27b7590..092f7bd0be 100644 --- a/test/dsl/dsl_advanced_model_construction.jl +++ b/test/dsl/dsl_advanced_model_construction.jl @@ -340,10 +340,23 @@ let k[1]*a+k[2], X[1] + V[1]*X[2] --> V[2]*W*Y + B*C end - @parameters k[1:3] a B + @parameters k[1:2] a B @variables (V(t))[1:2] W(t) @species (X(t))[1:2] Y(t) C(t) rx = Reaction(k[1]*a+k[2], [X[1], X[2]], [Y, C], [1, V[1]], [V[2] * W, B]) @named arrtest = ReactionSystem([rx], t) - arrtest == rn + @test arrtest == rn + + rn = @reaction_network twostate begin + @parameters k[1:2] + @species (X(t))[1:2] + (k[1],k[2]), X[1] <--> X[2] + end + + @parameters k[1:2] + @species (X(t))[1:2] + rx1 = Reaction(k[1], [X[1]], [X[2]]) + rx2 = Reaction(k[2], [X[2]], [X[1]]) + @named twostate = ReactionSystem([rx1, rx2], t) + @test twostate == rn end diff --git a/test/network_analysis/conservation_laws.jl b/test/network_analysis/conservation_laws.jl index e962252ae0..fe79920a02 100644 --- a/test/network_analysis/conservation_laws.jl +++ b/test/network_analysis/conservation_laws.jl @@ -343,7 +343,7 @@ end let # Prepares the model. t = default_t() - @species X(t)[1:2] + @species (X(t))[1:2] @parameters k[1:2] rxs = [ Reaction(k[1], [X[1]], [X[2]]), diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl index 97d32eb99f..3acdf008e3 100644 --- a/test/reactionsystem_core/reactionsystem.jl +++ b/test/reactionsystem_core/reactionsystem.jl @@ -391,7 +391,7 @@ end # Needs fix for https://github.com/JuliaSymbolics/Symbolics.jl/issues/842. let @parameters a - @species A(t) B(t) C(t)[1:2] + @species A(t) B(t) (C(t))[1:2] rx1 = Reaction(a, [A, C[1]], [C[2], B], [1, 2], [2, 3]) io = IOBuffer() show(io, rx1) diff --git a/test/upstream/mtk_problem_inputs.jl b/test/upstream/mtk_problem_inputs.jl index ffe2037519..c973596f03 100644 --- a/test/upstream/mtk_problem_inputs.jl +++ b/test/upstream/mtk_problem_inputs.jl @@ -180,7 +180,7 @@ end begin # Declares the model (with vector species/parameters, with/without default values, and observables). t = default_t() - @species X(t)[1:2] Y(t)[1:2] = [10.0, 20.0] XY(t)[1:2] + @species (X(t))[1:2] (Y(t))[1:2] = [10.0, 20.0] (XY(t))[1:2] @parameters p[1:2] d[1:2] = [0.2, 0.5] rxs = [ Reaction(p[1], [], [X[1]]),