Merge pull request #127 from isaacsas/update_for_catalyst14

isaacsas · web-flow · commit 84d52593b388 · 2024-07-20T20:52:47.000-04:00
Update for Catalyst 14
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "ReactionNetworkImporters"
 uuid = "b4db0fb7-de2a-5028-82bf-5021f5cfa881"
 authors = ["Samuel Isaacson <isaacsas@users.noreply.github.com>"]
-version = "0.14.1"
+version = "0.15.0"
 
 [deps]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
@@ -13,20 +13,18 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 [compat]
 Aqua = "0.8"
 CSVFiles = "1"
-Catalyst = "13"
+Catalyst = "14"
 DataFrames = "1.3"
 DataStructures = "0.18.13"
-DiffEqBase = "6.117"
+DiffEqBase = "6.151"
 LinearAlgebra = "1.10"
-ModelingToolkit = "8.51"
-OrdinaryDiffEq = "6.47"
-Plots = "1.30"
+ModelingToolkit = "9.24"
+OrdinaryDiffEq = "6.8"
 SafeTestsets = "0.1"
 SparseArrays = "1.10"
 Sundials = "4.4"
-Symbolics = "5.2"
+Symbolics = "5.33"
 Test = "1"
-TimerOutputs = "0.5"
 julia = "1.10"
 
 [extras]
@@ -36,11 +34,9 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [targets]
-test = ["Aqua", "CSVFiles", "DataFrames", "DiffEqBase", "LinearAlgebra", "OrdinaryDiffEq", "Plots", "SafeTestsets", "Sundials", "Test", "TimerOutputs"]
+test = ["Aqua", "CSVFiles", "DataFrames", "DiffEqBase", "LinearAlgebra", "OrdinaryDiffEq", "SafeTestsets", "Sundials", "Test"]
diff --git a/README.md b/README.md
@@ -51,7 +51,9 @@ 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`
+  - `rn`, a Catalyst `ReactionSystem`. **Note** this system is *not* marked
+    complete by default (see the [Catalyst docs](https://catalyst.sciml.ai/) for
+    a discussion of completeness of systems).
   - `u0`, 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
@@ -71,7 +73,7 @@ reaction system by
 
 ```julia
 using OrdinaryDiffEq, Catalyst
-rn = prnbng.rn
+rn = complete(prnbng.rn)   # get the reaction network and mark it complete
 tf = 100000.0
 oprob = ODEProblem(rn, Float64[], (0.0, tf), Float64[])
 sol = solve(oprob, Tsit5(), saveat = tf / 1000.0)
@@ -169,18 +171,18 @@ reaction rate expressions. These two types have the following fields:
     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`.
           * `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`.
-    
+
       + For `ComplexMatrixNetwork`
-        
+
           * `stoichmat`, the complex stoichiometry matrix [defined
             here](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Catalyst.complexstoichmat).
           * `incidencemat`, the complex incidence matrix [defined
@@ -194,14 +196,16 @@ reaction rate expressions. These two types have the following fields:
     [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/)
     `@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
+    variable of the reaction network. If not provided `Catalyst.default_t()` 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://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Catalyst.ReactionSystem)
-that represents the network.
+that represents the network. **Note**, `prn.rn` is not marked as complete by
+default and must be manually completed by setting `rn = complete(prn.rn)` before
+creating an `ODEProblem` or such, see the Catalyst docs for details.
 
 Dispatches are added if `substoich` and `prodstoich` both have the type
 `SparseMatrixCSC`in case of `MatrixNetwork` (or `stoichmat` and `incidencemat`
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -4,6 +4,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ReactionNetworkImporters = "b4db0fb7-de2a-5028-82bf-5021f5cfa881"
 
 [compat]
-Catalyst = "13.0"
+Catalyst = "14.0"
 Documenter = "1"
-ReactionNetworkImporters = "0.14"
+ReactionNetworkImporters = "0.15"
diff --git a/docs/make.jl b/docs/make.jl
@@ -1,8 +1,10 @@
 using Documenter
 using ReactionNetworkImporters
 
-cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true)
-cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
+docpath = Base.source_dir()
+assetpath = joinpath(docpath, "src", "assets")
+cp(joinpath(docpath, "Manifest.toml"), joinpath(assetpath, "Manifest.toml"), force = true)
+cp(joinpath(docpath, "Project.toml"), joinpath(assetpath, "Project.toml"), force = true)
 
 include("pages.jl")
 
@@ -25,8 +27,7 @@ makedocs(sitename = "ReactionNetworkImporters.jl",
                                   canonical = "https://docs.sciml.ai/ReactionNetworkImporters/stable/"),
          modules = [ReactionNetworkImporters],
          checkdocs = :exports, warnonly = [:missing_docs],
-         clean = true, doctest = false, linkcheck = true,
-         pages = pages)
+         clean = true, doctest = false, pages = pages)
 
 deploydocs(repo = "github.com/SciML/ReactionNetworkImporters.jl.git";
            push_preview = true)
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -35,7 +35,7 @@ Pkg.add("ReactionNetworkImporters")
 
   - See the [SciML Style Guide](https://github.com/SciML/SciMLStyle) for common coding practices and other style decisions.
   - There are a few community forums:
-    
+
       + The #diffeq-bridged and #sciml-bridged channels in the
         [Julia Slack](https://julialang.org/slack/)
       + The #diffeq-bridged and #sciml-bridged channels in the
@@ -63,7 +63,9 @@ 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`
+  - `rn`, a Catalyst `ReactionSystem`. **Note** this system is *not* marked
+    complete by default (see the [Catalyst docs](https://catalyst.sciml.ai/) for
+    a discussion of completeness of systems).
   - `u0`, 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
@@ -83,7 +85,7 @@ reaction system by
 
 ```julia
 using OrdinaryDiffEq, Catalyst
-rn = prnbng.rn
+rn = complete(prnbng.rn)   # get the reaction network and mark it complete
 tf = 100000.0
 oprob = ODEProblem(rn, Float64[], (0.0, tf), Float64[])
 sol = solve(oprob, Tsit5(), saveat = tf / 1000.0)
@@ -181,18 +183,18 @@ reaction rate expressions. These two types have the following fields:
     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`.
           * `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`.
-    
+
       + For `ComplexMatrixNetwork`
-        
+
           * `stoichmat`, the complex stoichiometry matrix [defined
             here](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Catalyst.complexstoichmat).
           * `incidencemat`, the complex incidence matrix [defined
@@ -206,14 +208,16 @@ reaction rate expressions. These two types have the following fields:
     [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/)
     `@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
+    variable of the reaction network. If not provided, `Catalyst.default_t()` 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://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Catalyst.ReactionSystem)
-that represents the network.
+that represents the network. **Note**, `prn.rn` is not marked as complete by
+default and must be manually completed by setting `rn = complete(prn.rn)` before
+creating an `ODEProblem` or such, see the Catalyst docs for details.
 
 Dispatches are added if `substoich` and `prodstoich` both have the type
 `SparseMatrixCSC`in case of `MatrixNetwork` (or `stoichmat` and `incidencemat`
diff --git a/src/parsing_routines_bngnetworkfiles.jl b/src/parsing_routines_bngnetworkfiles.jl
@@ -221,7 +221,7 @@ function loadrxnetwork(ft::BNGNetwork, rxfilename; name = gensym(:ReactionSystem
     file = open(rxfilename, "r")
     lines = readlines(file)
     idx = 1
-    t = Catalyst.DEFAULT_IV
+    t = Catalyst.default_t()
 
     print("Parsing parameters...")
     ptoids, pvals, idx = parse_params(ft, lines, idx)
@@ -250,7 +250,7 @@ function loadrxnetwork(ft::BNGNetwork, rxfilename; name = gensym(:ReactionSystem
     print("Creating species and parameters for evaluating expressions...")
     opmod = Module()
     Base.eval(opmod, :(using Catalyst))
-    Base.eval(opmod, :(t = Catalyst.DEFAULT_IV))
+    Base.eval(opmod, :(t = Catalyst.default_t()))
     for p in ps
         psym = nameof(p)
         Base.eval(opmod, :($(psym) = $p))
@@ -275,7 +275,7 @@ function loadrxnetwork(ft::BNGNetwork, rxfilename; name = gensym(:ReactionSystem
     # setup default values / expressions for params and initial conditions
     defmap, pmap, u0map = exprs_to_defs(opmod, ptoids, pvals, specs, u0exprs)
 
-    # build the model    
+    # build the model
     rn = ReactionSystem(rxs, t, specs, ps; name = name, observed = obseqs,
                         defaults = defmap, kwargs...)
 
diff --git a/src/parsing_routines_matrixnetworks.jl b/src/parsing_routines_matrixnetworks.jl
@@ -42,7 +42,7 @@ function loadrxnetwork(mn::MatrixNetwork{S, T, U, V, W, X};
     numspecs = sz[1]
     numrxs = sz[2]
 
-    t = (mn.t === nothing) ? Catalyst.DEFAULT_IV : mn.t
+    t = (mn.t === nothing) ? Catalyst.default_t() : mn.t
     species = isempty(mn.species) ? [funcsym(:S, t, i) for i in 1:numspecs] : mn.species
 
     # create the reactions
@@ -89,7 +89,7 @@ function loadrxnetwork(mn::MatrixNetwork{S, T, U, V, W, X};
     numspecs = sz[1]
     numrxs = sz[2]
 
-    t = (mn.t === nothing) ? Catalyst.DEFAULT_IV : mn.t
+    t = (mn.t === nothing) ? Catalyst.default_t() : mn.t
     species = isempty(mn.species) ? [funcsym(:S, t, i) for i in 1:numspecs] : mn.species
 
     # create the reactions
@@ -182,7 +182,7 @@ function loadrxnetwork(cmn::ComplexMatrixNetwork{S, T, U, V, W, X};
     @assert all(∈([-1, 0, 1]), cmn.incidencemat)
     numrxs = size(cmn.incidencemat, 2)
 
-    t = (cmn.t === nothing) ? Catalyst.DEFAULT_IV : cmn.t
+    t = (cmn.t === nothing) ? Catalyst.default_t() : cmn.t
     species = isempty(cmn.species) ? [funcsym(:S, t, i) for i in 1:numspecs] : cmn.species
 
     rxs = Vector{Reaction}(undef, numrxs)
@@ -227,7 +227,7 @@ function loadrxnetwork(cmn::ComplexMatrixNetwork{S, T, U, V, W, X};
     @assert all(∈([-1, 0, 1]), cmn.incidencemat)
     numrxs = size(cmn.incidencemat, 2)
 
-    t = (cmn.t === nothing) ? Catalyst.DEFAULT_IV : cmn.t
+    t = (cmn.t === nothing) ? Catalyst.default_t() : cmn.t
     species = isempty(cmn.species) ? [funcsym(:S, t, i) for i in 1:numspecs] : cmn.species
 
     rxs = Vector{Reaction}(undef, numrxs)
diff --git a/test/test_higherorder_odes.jl b/test/test_higherorder_odes.jl
@@ -1,5 +1,6 @@
-using DiffEqBase, Catalyst, Plots, OrdinaryDiffEq, DataFrames, CSVFiles, LinearAlgebra
+using DiffEqBase, Catalyst, OrdinaryDiffEq, DataFrames, CSVFiles, LinearAlgebra
 using ReactionNetworkImporters
+# using Plots
 
 # parameters
 doplot = false
@@ -18,7 +19,7 @@ println("done")
 
 # load the BNG reaction network in DiffEqBio
 prnbng = loadrxnetwork(BNGNetwork(), fname)
-rn = prnbng.rn
+rn = complete(prnbng.rn)
 boprob = ODEProblem(rn, Float64[], (0.0, tf), Float64[])
 
 # Test that u0 == u₀ (used when the u0 indexing was introduced).
@@ -28,16 +29,16 @@ boprob = ODEProblem(rn, Float64[], (0.0, tf), Float64[])
 @unpack A = rn
 Asol = gdatdf[!, :A]
 
-# note solvers run _much_ faster the second time 
+# note solvers run _much_ faster the second time
 bsol = solve(boprob, Tsit5(), abstol = 1e-12, reltol = 1e-12, saveat = tf / nsteps);
 
-if doplot
-    plotlyjs()
-    p1 = plot(gdatdf[!, :time], gdatdf[!, :A], label = :A_BNG)
-    plot!(p1, bsol.t, bsol[A], label = :A_DE)
-    display(p1)
-    println("Err = ", norm(gdatdf[!, :A] - bsol[Aid, :], Inf))
-end
+# if doplot
+#     plotlyjs()
+#     p1 = plot(gdatdf[!, :time], gdatdf[!, :A], label = :A_BNG)
+#     plot!(p1, bsol.t, bsol[A], label = :A_DE)
+#     display(p1)
+#     println("Err = ", norm(gdatdf[!, :A] - bsol[Aid, :], Inf))
+# end
 
 @test all(bsol.t .== gdatdf[!, :time])
 @test all(abs.(gdatdf[!, :A] - bsol[A]) .< 1e-6 .* abs.(gdatdf[!, :A]))
diff --git a/test/test_nullrxs_odes.jl b/test/test_nullrxs_odes.jl
@@ -19,16 +19,16 @@ println("done")
 
 # load the BNG reaction network in DiffEqBio
 prnbng = loadrxnetwork(BNGNetwork(), fname)
-rn = prnbng.rn;
-u0 = Float64[];
-p = Float64[];
+rn = complete(prnbng.rn)
+u0 = Float64[]
+p = Float64[]
 boprob = ODEProblem(rn, u0, (0.0, tf), p)
 
 # Test that u0 == u₀ (used when the u0 indexing was introduced).
 @test isequal(prnbng.u0, prnbng.u₀)
 
-# note solvers run _much_ faster the second time 
-bsol = solve(boprob, Tsit5(), abstol = 1e-12, reltol = 1e-12, saveat = tf / nsteps);
+# note solvers run _much_ faster the second time
+bsol = solve(boprob, Tsit5(), abstol = 1e-12, reltol = 1e-12, saveat = tf / nsteps)
 
 if doplot
     plotlyjs()
diff --git a/test/test_repressilator_odes.jl b/test/test_repressilator_odes.jl
@@ -1,6 +1,8 @@
-using DiffEqBase, Catalyst, Plots, OrdinaryDiffEq, DataFrames, CSVFiles, LinearAlgebra
+using DiffEqBase, Catalyst, OrdinaryDiffEq, DataFrames, CSVFiles, LinearAlgebra
 using ReactionNetworkImporters
 
+# using Plots
+
 # parameters
 doplot = false
 networkname = "Repressilator"
@@ -18,9 +20,9 @@ println("done")
 
 # load the BNG reaction network in DiffEqBio
 prnbng = loadrxnetwork(BNGNetwork(), fname)
-rn = prnbng.rn;
-u0 = Float64[];
-p = Float64[];
+rn = complete(prnbng.rn)
+u0 = Float64[]
+p = Float64[]
 boprob = ODEProblem(rn, u0, (0.0, tf), p)
 
 # Test that u0 == u₀ (used when the u0 indexing was introduced).
@@ -29,16 +31,16 @@ boprob = ODEProblem(rn, u0, (0.0, tf), p)
 # BNG simulation data testing
 bngsol = gdatdf[!, :pTetR]
 
-# note solvers run _much_ faster the second time 
+# note solvers run _much_ faster the second time
 bsol = solve(boprob, Tsit5(), abstol = 1e-12, reltol = 1e-12, saveat = tf / nsteps);
 @unpack pTetR = rn
 
-if doplot
-    plotlyjs()
-    p1 = plot(gdatdf[!, :time], gdatdf[!, :pTetR], label = :BNGPTetR)
-    plot!(p1, bsol.t, bsol[pTetR], label = :DEBPTetR)
-    display(p1)
-end
+# if doplot
+#     plotlyjs()
+#     p1 = plot(gdatdf[!, :time], gdatdf[!, :pTetR], label = :BNGPTetR)
+#     plot!(p1, bsol.t, bsol[pTetR], label = :DEBPTetR)
+#     display(p1)
+# end
 
 @test all(bsol.t .== gdatdf[!, :time])
 @test all(abs.(gdatdf[!, :pTetR] - bsol[pTetR]) .< 1e-6 .* abs.(gdatdf[!, :pTetR]))
