Merge pull request #1129 from anandijain/aj/rxnsys_defaults

adding defaults to ReactionSystem constructors

Author: ChrisRackauckas

src/systems/reaction/reactionsystem.jl

@@ -149,44 +149,38 @@ struct ReactionSystem <: AbstractSystem
     name::Symbol
     """systems: The internal systems"""
     systems::Vector
+    """
+    defaults: The default values to use when initial conditions and/or
+    parameters are not supplied in `ODEProblem`.
+    """
+    defaults::Dict
 
-    function ReactionSystem(eqs, iv, states, ps, observed, name, systems)
+    function ReactionSystem(eqs, iv, states, ps, observed, name, systems, defaults)
         iv′ = value(iv)
         states′ = value.(states)
         ps′ = value.(ps)
         check_variables(states′, iv′)
         check_parameters(ps′, iv′)
-        new(collect(eqs), iv′, states′, ps′, observed, name, systems)
+        new(collect(eqs), iv′, states′, ps′, observed, name, systems, defaults)
     end
 end
 
 function ReactionSystem(eqs, iv, species, params;
                         observed = [],
                         systems = [],
-                        name = gensym(:ReactionSystem))
+                        name = gensym(:ReactionSystem),
+                        default_u0=Dict(),
+                        default_p=Dict(),
+                        defaults=_merge(Dict(default_u0), Dict(default_p)))
 
     #isempty(species) && error("ReactionSystems require at least one species.")
-    ReactionSystem(eqs, iv, species, params, observed, name, systems)
+    ReactionSystem(eqs, iv, species, params, observed, name, systems, defaults)
 end
 
 function ReactionSystem(iv; kwargs...)
     ReactionSystem(Reaction[], iv, [], []; kwargs...)
 end
 
-function equations(sys::ModelingToolkit.ReactionSystem)
-    eqs = get_eqs(sys)
-    systems = get_systems(sys)
-    if isempty(systems)
-        return eqs
-    else
-        eqs = [eqs;
-               reduce(vcat,
-                      namespace_equations.(get_systems(sys));
-                      init=[])]
-        return eqs
-    end
-end
-
 """
     oderatelaw(rx; combinatoric_ratelaw=true)
 
@@ -419,7 +413,7 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem;
                       name=nameof(rs), combinatoric_ratelaws=true, include_zero_odes=true, kwargs...)
     eqs     = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws, include_zero_odes=include_zero_odes)
     systems = map(sys -> (sys isa ODESystem) ? sys : convert(ODESystem, sys), get_systems(rs))
-    ODESystem(eqs, get_iv(rs), get_states(rs), get_ps(rs); name=name, systems=systems, kwargs...)
+    ODESystem(eqs, get_iv(rs), get_states(rs), get_ps(rs); name=name, systems=systems, defaults=get_defaults(rs), kwargs...)
 end
 
 """
@@ -439,7 +433,7 @@ function Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem;
                       name=nameof(rs), combinatoric_ratelaws=true, include_zero_odes=true, kwargs...)
     eqs     = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws, as_odes=false, include_zero_odes=include_zero_odes)
     systems = convert.(NonlinearSystem, get_systems(rs))
-    NonlinearSystem(eqs, get_states(rs), get_ps(rs); name=name, systems=systems, kwargs...)
+    NonlinearSystem(eqs, get_states(rs), get_ps(rs); name=name, systems=systems, defaults=get_defaults(rs), kwargs...)
 end
 
 """
@@ -487,6 +481,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
               (noise_scaling===nothing) ? get_ps(rs) : union(get_ps(rs), toparam(noise_scaling));
               name=name, 
               systems=systems,
+              defaults=get_defaults(rs), 
               kwargs...)
 end
 
@@ -507,7 +502,7 @@ function Base.convert(::Type{<:JumpSystem},rs::ReactionSystem;
                       name=nameof(rs), combinatoric_ratelaws=true, kwargs...)
     eqs     = assemble_jumps(rs; combinatoric_ratelaws=combinatoric_ratelaws)
     systems = convert.(JumpSystem, get_systems(rs))
-    JumpSystem(eqs, get_iv(rs), get_states(rs), get_ps(rs); name=name, systems=systems, kwargs...)
+    JumpSystem(eqs, get_iv(rs), get_states(rs), get_ps(rs); name=name, systems=systems, defaults=get_defaults(rs), kwargs...)
 end
 
 

test/reactionsystem.jl

@@ -34,6 +34,35 @@ show(io, rs)
 str = String(take!(io))
 @test count(isequal('\n'), str) < 30
 
+# defaults test
+def_p = [ki => float(i) for (i, ki) in enumerate(k)]
+def_u0 = [A => 0.5, B => 1., C=> 1.5, D => 2.0]
+defs = merge(Dict(def_p), Dict(def_u0))
+
+rs = ReactionSystem(rxs,t,[A,B,C,D],k; defaults=defs)
+odesys = convert(ODESystem,rs)
+sdesys = convert(SDESystem,rs)
+js = convert(JumpSystem,rs)
+nlsys = convert(NonlinearSystem,rs)
+
+@test ModelingToolkit.get_defaults(rs) ==
+       ModelingToolkit.get_defaults(odesys) ==
+       ModelingToolkit.get_defaults(sdesys) ==
+       ModelingToolkit.get_defaults(js) ==
+       ModelingToolkit.get_defaults(nlsys) ==
+       defs
+
+u0map = [A=>5.] # was 0.5
+pmap = [k[1]=>5.] # was 1.
+prob = ODEProblem(rs, u0map, (0,10.), pmap)
+@test prob.p[1] == 5.
+@test prob.u0[1] == 5.
+u0 = [10., 11., 12., 13.]
+ps = [float(x) for x in 100:119]
+prob = ODEProblem(rs, u0, (0,10.), ps)
+@test prob.p == ps
+@test prob.u0 == u0
+
 # hard coded ODE rhs
 function oderhs(u,k,t)
        A = u[1]; B = u[2]; C = u[3]; D = u[4];