Merge pull request #72 from isaacsas/use_defaults

Use defaults

Project.toml

@@ -1,7 +1,7 @@
 name = "ReactionNetworkImporters"
 uuid = "b4db0fb7-de2a-5028-82bf-5021f5cfa881"
 authors = ["Samuel Isaacson <isaacsas@users.noreply.github.com>"]
-version = "0.12.0"
+version = "0.13.0"
 
 [deps]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"

README.md

@@ -31,28 +31,34 @@ prnbng = loadrxnetwork(BNGNetwork(), fname)
 Here `BNGNetwork` is a type specifying the file format that is being loaded.
 `prnbng` is a `ParsedReactionNetwork` structure with the following fields:
 - `rn`, a Catalyst `ReactionSystem`
-- `u₀`, the initial condition (as a `Vector{Float64}`)
-- `p`, the parameter vector (as a `Vector{Float64}`)
-- `paramexprs`, the parameter vector as a mix of `Numbers`, `Symbols` and
-  `Exprs`. `p` is generated by evaluation of these expressions and symbols.
-- `varstonames`, a `Dict` mapping from the internal `Symbol` of a species used
-  in the generated `ReactionSystem` to a `String` generated from the name in the
-  .net file. This is necessary as BioNetGen can generate exceptionally long
-  species names, involving characters that lead to malformed species names when
-  used with `Catalyst`.
+- `u₀`, a `Dict` mapping initial condition symbolic variables to numeric values
+  and/or symbolic expressions.
+- `p`, a `Dict` mapping parameter symbolic variables to numeric values and/or
+  symbolic expressions.
+- `varstonames`, a `Dict` mapping the internal symbolic variable of a species
+  used in the generated `ReactionSystem` to a `String` generated from the name
+  in the .net file. This is necessary as BioNetGen can generate exceptionally
+  long species names, involving characters that lead to malformed species names
+  when used with `Catalyst`.
 - `groupstosyms`, a `Dict` mapping the `String`s representing names for any
-  species groups defined in the .net file to a vector of symbolic variables
-  representing the corresponding observables in the generated `ReactionSystem`.
+  groups defined in the BioNetGen file to the corresponding symbolic variable
+  representing the `ModelingToolkit` symbolic observable associated with the
+  group.
 
 Given `prnbng`, we can construct and solve the corresponding ODE model for the
 reaction system by
 ```julia
 using OrdinaryDiffEq, Catalyst
 rn = prnbng.rn
 tf = 100000.0
-oprob = ODEProblem(rn, prnbng.u₀, (0.,tf), prnbng.p)
+oprob = ODEProblem(rn, Float64[], (0.,tf), Float64[])
 sol = solve(oprob, Tsit5(), saveat=tf/1000.)
 ```
+Note that we specify empty parameter and initial condition vectors as these are
+already stored in the generated `ReactionSystem`, `rn`. A `Dict` mapping each
+symbolic species and parameter to its initial value or symbolic expression can
+be obtained using `ModelingToolkit.defaults(rn)`. 
+
 See the [Catalyst documentation](https://catalyst.sciml.ai/dev/) for how to
 generate ODE, SDE, jump and other types of models.
 
@@ -85,7 +91,8 @@ substoich =[ 2  0  1  0  0;
 prodstoich =  [0  2  0  1  3;
                 1  0  0  1  0;
                 0  0  1  0  0]
-mn= MatrixNetwork(pars, substoich, prodstoich; species=species, params=pars) # matrix network
+mn= MatrixNetwork(pars, substoich, prodstoich; species=species, 
+                  params=pars) # a matrix network
 prn = loadrxnetwork(mn) # dense version
 
 # test the two networks are the same
@@ -101,7 +108,8 @@ incidencemat  = [-1   1   0   0   0;
                  0   0   1  -1   0;
                  0   0   0   0  -1;
                  0   0   0   0   1]
-cmn= ComplexMatrixNetwork(pars, stoichmat, incidencemat; species=species, params=pars)  # complex matrix network
+cmn= ComplexMatrixNetwork(pars, stoichmat, incidencemat; species=species, 
+                          params=pars)  # a complex matrix network
 prn = loadrxnetwork(cmn)
 
 # test the two networks are the same
@@ -110,55 +118,66 @@ prn = loadrxnetwork(cmn)
 
 The basic usages are
 ```julia
-mn = MatrixNetwork(rateexprs, substoich, prodstoich; species=Any[], params=Any[], t=nothing)
+mn = MatrixNetwork(rateexprs, substoich, prodstoich; species=Any[], 
+                   params=Any[], t=nothing)
 prn = loadrxnetwork(mn::MatrixNetwork) 
 
-cmn = ComplexMatrixNetwork(rateexprs, stoichmat, incidencemat; species=Any[], params=Any[], t=nothing)
+cmn = ComplexMatrixNetwork(rateexprs, stoichmat, incidencemat; species=Any[], 
+                           params=Any[], t=nothing)
 prn = loadrxnetwork(cmn::ComplexMatrixNetwork)
 ```
-Here `MatrixNetwork` and `ComplexMatrixNetwork` are the types, which select that we are
-constructing a substrate/product stoichiometric matrix-based or a reaction complex matrix-based 
-stoichiometric representation as input. See the [Catalyst.jl API](https://catalyst.sciml.ai/dev/api/catalyst_api/) 
-for more discussion on these matrix representations, and how Catalyst handles symbolic reaction rate expressions. 
-These two types have the following fields:
-- `rateexprs`, any valid [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl) expression 
-  for the rates, or any basic number type. This can be a hardcoded rate constant like `1.0`, a parameter
-  like `k1` above, or an general Symbolics expression involving parameters and species like
-  `k*A`. Note, the reaction `A+B --> C` with rate (i.e. `rateexprs` entry) `k*B` would have rate law
-  `k*A*B^2`.
+Here `MatrixNetwork` and `ComplexMatrixNetwork` are the types, which select that
+we are constructing a substrate/product stoichiometric matrix-based or a
+reaction complex matrix-based stoichiometric representation as input. See the
+[Catalyst.jl API](https://catalyst.sciml.ai/dev/api/catalyst_api/) for more
+discussion on these matrix representations, and how Catalyst handles symbolic
+reaction rate expressions. These two types have the following fields:
+- `rateexprs`, any valid
+  [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl) expression for
+  the rates, or any basic number type. This can be a hardcoded rate constant
+  like `1.0`, a parameter like `k1` above, or an general Symbolics expression
+  involving parameters and species like `k*A`. 
 - matrix inputs 
   - For `MatrixNetwork`
     - `substoich`, a number of species by number of reactions matrix with entry
-      `(i,j)` giving the stoichiometric coefficient of species `i` as a substrate in
-      reaction `j`. 
+      `(i,j)` giving the stoichiometric coefficient of species `i` as a
+      substrate in reaction `j`. 
     - `prodstoich`, a number of species by number of reactions matrix with entry
-      `(i,j)` giving the stoichiometric coefficient of species `i` as a product in
-      reaction `j`.
+      `(i,j)` giving the stoichiometric coefficient of species `i` as a product
+      in reaction `j`.
   - For `ComplexMatrixNetwork`  
-    - `stoichmat`, the complex stoichiometry matrix [defined here](https://catalyst.sciml.ai/dev/api/catalyst_api/#Catalyst.complexstoichmat).
-    - `incidencemat`, the complex incidence matrix [defined here](https://catalyst.sciml.ai/dev/api/catalyst_api/#Catalyst.reactioncomplexes). 
-- `species`, an optional vector of symbolic expressions representing each species in
-  the network. Can be constructed using the Symbolics.jl `@variables` macro. Each species
-  should be dependent on the same time variable (`t` in the example above).
-- `parameters`, an optional vector of symbolic parameters representing each
-  parameter in the network. Can be constructed with the [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl)
-  `@parameters` macro.
-- `t`, an optional Symbolics.jl variable representing time as the independent variable of the reaction network
-
-For both input types, `loadrxnetwork` returns a `ParsedReactionNetwork`, `prn`, with only 
-the field, `prn.rn`, filled in. `prn.rn` corresponds to the generated 
-[Catalyst.jl `ReactionSystem`](https://catalyst.sciml.ai/dev/api/catalyst_api/#Catalyst.ReactionSystem)
+    - `stoichmat`, the complex stoichiometry matrix [defined
+      here](https://catalyst.sciml.ai/dev/api/catalyst_api/#Catalyst.complexstoichmat).
+    - `incidencemat`, the complex incidence matrix [defined
+      here](https://catalyst.sciml.ai/dev/api/catalyst_api/#Catalyst.reactioncomplexes).      
+- `species`, an optional vector of symbolic variables representing each species
+  in the network. Can be constructed using the Symbolics.jl `@variables` macro.
+  Each species should be dependent on the same time variable (`t` in the example
+  above). 
+- `parameters`, a vector of symbolic variables representing each parameter in
+  the network. Can be constructed with the
+  [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl)
+  `@parameters` macro. If no parameters are used it is an optional keyword.
+- `t`, an optional Symbolics.jl variable representing time as the independent
+  variable of the reaction network. If not provided `Catalyst.DEFAULT_IV` is
+  used to determine the default time variable.
+
+For both input types, `loadrxnetwork` returns a `ParsedReactionNetwork`, `prn`,
+with only the field, `prn.rn`, filled in. `prn.rn` corresponds to the generated
+[Catalyst.jl
+`ReactionSystem`](https://catalyst.sciml.ai/dev/api/catalyst_api/#Catalyst.ReactionSystem)
 that represents the network.
 
 Dispatches are added if `substoich` and `prodstoich` both have the type
-`SparseMatrixCSC`in case of `MatrixNetwork` (or `stoichmat` and `incidencemat` both have the 
-type `SparseMatrixCSC` in case of `ComplexMatrixNetwork`), in which case they are efficiently 
-iterated through using the `SparseArrays` interface.
+`SparseMatrixCSC`in case of `MatrixNetwork` (or `stoichmat` and `incidencemat`
+both have the type `SparseMatrixCSC` in case of `ComplexMatrixNetwork`), in
+which case they are efficiently iterated through using the `SparseArrays`
+interface.
 
 If the keyword argument `species` is not set, the resulting reaction network
 will simply name the species `S1`, `S2`,..., `SN` for a system with `N` total
-species. `params` defaults to an empty vector, so that it does not
-need to be set for systems with no parameters.
+species. `params` defaults to an empty vector, so that it does not need to be
+set for systems with no parameters.
 
 <!-- ### Loading a RSSA format network file
 As the licensing is unclear we can not redistribute any example RSSA formatted networks. They can be downloaded from the model collection link listed above. Assuming you've saved both a reaction network file and corresponding initial condition file, they can be loaded as

examples/test_bcr_odes.jl

@@ -27,12 +27,13 @@ reset_timer!(to)
 # BioNetGen network
 @timeit to "bionetgen" prnbng = loadrxnetwork(BNGNetwork(),fname); 
 show(to)
-rnbng = prnbng.rn; u0 = prnbng.u₀; p = prnbng.p; 
+rnbng = prnbng.rn
 @timeit to "bODESystem" bosys = convert(ODESystem, rnbng)
 show(to)
-@timeit to "bODEProb" boprob = ODEProblem(bosys, Pair.(species(rnbng),u0), (0.,tf), Pair.(parameters(rnbng),p))
+@timeit to "bODEProb" boprob = ODEProblem(bosys, Float64[], (0.,tf), Float64[])
 show(to)
-u = copy(u0);
+u = zeros(numspecies(rnbng))
+p = zeros(length(parameters(rnbng)))
 du = similar(u);
 @timeit to "f1" boprob.f(du,u,p,0.)
 @timeit to "f2" boprob.f(du,u,p,0.)
@@ -45,15 +46,6 @@ end
 show(to)
 println()
 
-# BNG simulation results for Activated Syk
-asykgroups = prnbng.groupstoids[:Activated_Syk]
-#asyksyms = findall(x -> x ∈ asykgroups, rnbng.syms_to_ints)
-# asynbng = zeros(length(gdatdf[:time]))
-# for sym in asyksyms
-#     global asynbng
-#     asynbng += gdatdf[sym]
-# end
-
 # DiffEq solver 
 #reset_timer!(to); 
 #@timeit to "BNG_CVODE_BDF-LU-1" begin bsol = solve(boprob, CVODE_BDF(), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
@@ -67,18 +59,17 @@ asykgroups = prnbng.groupstoids[:Activated_Syk]
 @timeit to "KenCarp4-2" begin bsol = solve(boprob,KenCarp4(autodiff=false), abstol=1e-8, reltol=1e-8, saveat=1.); end; show(to)
 
 # Activated Syk from DiffEq
-basyk = sum(bsol[asykgroups,:], dims=1)
+@unpack Activated_Syk = rnbng
 
 if doplot
     plotlyjs()
     plot(gdatdf[!,:time][2:end], gdatdf[!,:Activated_Syk][2:end], xscale=:log10, label="AsykGroup", linestyle=:dot)
-#     # plot!(cdatdf[:time][2:end], asynbng[2:end], xscale=:log10, label=:AsykSum)
-    plot!(bsol.t[2:end], basyk'[2:end], label="AsykDEBio", xscale=:log10)
+    plot!(bsol.t[2:end], sol[Activated_Syk][2:end], label="Activated_Syk", xscale=:log10)
 end
 
 # test the error, note may be large in abs value though small relatively
 # #norm(gdatdf[:Activated_Syk] - asynbng, Inf)
 # #norm(asynbng - basyk', Inf)
-norm(gdatdf[!,:Activated_Syk] - basyk', Inf)
+norm(gdatdf[!,:Activated_Syk] - sol[Activated_Syk], Inf)
 
 #@assert all(abs.(gdatdf[!,:Activated_Syk] - basyk') .< 1e-6 * abs.(gdatdf[:Activated_Syk]))

src/ReactionNetworkImporters.jl

@@ -23,29 +23,25 @@ struct ParsedReactionNetwork
     "Catalyst Network"
     rn::ReactionSystem
 
-    "Initial Conditions"
+    "Dict mapping initial condition symbolic variables to values."
     u₀
 
-    "Parameters"
+    "Dict mapping parameter symbolic variables to values."
     p
 
-    "Expressions for the Parameters"
-    paramexprs
-
-    "Dict from symbolic variable in species(rn) to full string for species name"
+    "Dict mapping symbolic variable for species names to full string for species name"
     varstonames
 
     "Dict from group name (as string) to corresponding symbolic variable"
     groupstosyms
 
 end
-ParsedReactionNetwork(rn::ReactionSystem, u₀; p=nothing, paramexprs=nothing, varstonames=nothing, groupstosyms=nothing) = 
-                        ParsedReactionNetwork(rn, u₀, p, paramexprs, varstonames, groupstosyms)
+ParsedReactionNetwork(rn::ReactionSystem; u₀=nothing, p=nothing, varstonames=nothing, groupstosyms=nothing) = 
+                        ParsedReactionNetwork(rn, u₀, p, varstonames, groupstosyms)
 
 export BNGNetwork, MatrixNetwork, ParsedReactionNetwork, ComplexMatrixNetwork
 
 # parsers
-#include("parsing_routines_rssafiles.jl")
 include("parsing_routines_bngnetworkfiles.jl")
 include("parsing_routines_matrixnetworks.jl")
 

src/parsing_routines_bngnetworkfiles.jl

@@ -181,19 +181,33 @@ function parse_groups(ft::BNGNetwork, lines, idx, shortsymstosyms, idstoshortsym
     obseqs,namestosyms,idx
 end
 
-function exprs_to_nums(ptoids, pvals, u0exprs)
-    p    = zeros(Float64, length(ptoids))
-    pmod = Module()
+function exprs_to_defs(opmod, ptoids, pvals, specs, u0exprs)
+    pmap = Dict()
     for (psym,pid) in ptoids
-        p[pid] = Base.eval(pmod, :($psym = $(pvals[pid])))
+        pvar      = getproperty(opmod, psym)
+        parsedval = pvals[pid]
+        if (parsedval isa Expr) || (parsedval isa Symbol)
+            pval = Base.eval(opmod, :($parsedval))
+        else 
+            @assert parsedval isa Number
+            pval = parsedval
+        end
+        push!(pmap, pvar => pval)
     end
 
-    u0 = zeros(Float64, length(u0exprs))
+    u0map = Dict()
     for (i,u0expr) in enumerate(u0exprs)
-        u0[i] = Base.eval(pmod, :($u0expr))
+        uvar = specs[i]
+        if (u0expr isa Expr) || (u0expr isa Symbol)
+            u0val = Base.eval(opmod, :($u0expr))
+        else 
+            @assert u0expr isa Number
+            u0val = u0expr
+        end
+        push!(u0map, uvar => u0val)
     end
 
-    p,u0
+    union(pmap,u0map),pmap,u0map
 end
 
 
@@ -253,12 +267,15 @@ function loadrxnetwork(ft::BNGNetwork, rxfilename; name = gensym(:ReactionSystem
 
     close(file)
 
-    rn = ReactionSystem(rxs, t, specs, ps; name=name, observed=obseqs, kwargs...)
+    # setup default values / expressions for params and initial conditions
+    defmap,pmap,u0map = exprs_to_defs(opmod, ptoids, pvals, specs, u0exprs)
+
+    # build the model    
+    rn = ReactionSystem(rxs, t, specs, ps; name=name, observed=obseqs, defaults=defmap, kwargs...)
 
     # get numeric values for parameters and u₀
-    p,u₀ = exprs_to_nums(ptoids, pvals, u0exprs)
     sm = speciesmap(rn)
     @assert all( sm[funcsym(sym,t)] == i for (i,sym) in enumerate(idstoshortsyms) )
 
-    ParsedReactionNetwork(rn, u₀; p = p, paramexprs = pvals, varstonames = shortsymstosyms, groupstosyms = groupstosyms)
+    ParsedReactionNetwork(rn; u₀=u0map, p=pmap, varstonames=shortsymstosyms, groupstosyms=groupstosyms)
 end