Merge pull request #290 from SciML/s/newmtk2

Updates for upstream breaking changes

Project.toml

@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "6.3.1"
+version = "6.4"
 
 [deps]
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
@@ -13,7 +13,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Catlab = "0.7.3, 0.8, 0.9, 0.10"
 Latexify = "0.13.5, 0.14"
 MacroTools = "0.5.5"
-ModelingToolkit = "4.4.1"
+ModelingToolkit = "5"
 Reexport = "0.2, 1.0"
 julia = "1.3"
 

docs/src/tutorials/generated_systems.md

@@ -15,18 +15,18 @@ API method listed first:
 
 Each `Reaction` within `reactions(rn)` has a number of subfields. For `rx` a
 `Reaction` we have:
-* `rx.substrates`, a vector of `ModelingToolkit.Term`s storing each
+* `rx.substrates`, a vector of ModelingToolkit expressions storing each
   substrate variable.
-* `rx.products`, a vector of `ModelingToolkit.Term`s storing each product
+* `rx.products`, a vector of ModelingToolkit expressions storing each product
   variable.
 * `rx.substoich`, a vector storing the corresponding integer stoichiometry of
   each substrate species in `rx.substrates`.
 * `rx.prodstoich`, a vector storing the corresponding integer stoichiometry of
   each product species in `rx.products`.
-* `rx.rate`, a `Number, `ModelingToolkit.Sym` or `ModelingToolkit.Term`
+* `rx.rate`, a `Number, `ModelingToolkit.Sym` or ModelingToolkit expression
   representing the reaction rate. E.g., for a reaction like `k*X, Y --> X+Y`,
   we'd have `rate = k*X`.
-* `rx.netstoich`, a vector of pairs mapping the `ModelingToolkit.Term` for
+* `rx.netstoich`, a vector of pairs mapping the ModelingToolkit expression for
   each species that changes numbers by the reaction to how much it changes. E.g.,
   for `k, X + 2Y --> X + W`, we'd have `rx.netstoich = [Y(t) => -2, W(t) => 1]`.
 * `rx.only_use_rate`, a boolean that is `true` if the reaction was made with

src/Catalyst.jl

@@ -1,7 +1,7 @@
 module Catalyst
 
 using Reexport, ModelingToolkit
-using ModelingToolkit: Symbolic, value
+using ModelingToolkit: Symbolic, value, istree
 @reexport using ModelingToolkit
 import MacroTools
 import Base: (==), merge!, merge

src/networkapi.jl

@@ -284,7 +284,7 @@ id of the parameter within the system.
 """
 function addparam!(network::ReactionSystem, p::Symbolic; disablechecks=false)
     # we don't check subsystems since we will add it to the top-level system...
-    if p isa Term && !(p.op isa Sym)
+    if istree(p) && !(operation(p) isa Sym)
         error("If the passed in parameter is an expression, it must correspond to an underlying Variable.")
     end
     curidx = disablechecks ? nothing : findfirst(S -> isequal(S, p), network.ps)

test/latexify.jl

@@ -38,11 +38,11 @@ end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 v5 K5 n5 k1 k2 k3 k4 k5 k6 d1 d2 d3 d4 d
 @test latexify(r) == replace(
 raw"\begin{align}
 \require{mhchem}
-\ce{ \varnothing &->[\frac{v1 K1^{n1}}{K1^{n1} + \left( \mathrm{X2}\left( t \right) \right)^{n1}} \frac{v1 \left( \mathrm{X4}\left( t \right) \right)^{n1}}{K1^{n1} + \left( \mathrm{X4}\left( t \right) \right)^{n1}}] X1}\\
-\ce{ \varnothing &->[\frac{v2 \left( \mathrm{X5}\left( t \right) \right)^{n2}}{K2^{n2} + \left( \mathrm{X5}\left( t \right) \right)^{n2}}] X2}\\
-\ce{ \varnothing &->[\frac{v3 \left( \mathrm{X3}\left( t \right) \right)^{n3}}{K3^{n3} + \left( \mathrm{X3}\left( t \right) \right)^{n3}}] X3}\\
-\ce{ \varnothing &->[\frac{v4 K4^{n4}}{K4^{n4} + \left( \mathrm{X1}\left( t \right) \right)^{n4}}] X4}\\
-\ce{ \varnothing &->[\frac{v5 \left( \mathrm{X2}\left( t \right) \right)^{n5}}{K5^{n5} + \left( \mathrm{X2}\left( t \right) \right)^{n5}}] X5}\\
+\ce{ \varnothing &->[K1^{n1} v1^{2} \left( \mathrm{X4}\left( t \right) \right)^{n1} \mathrm{inv}\left( K1^{n1} + \left( \mathrm{X2}\left( t \right) \right)^{n1} \right) \mathrm{inv}\left( K1^{n1} + \left( \mathrm{X4}\left( t \right) \right)^{n1} \right)] X1}\\
+\ce{ \varnothing &->[v2 \left( \mathrm{X5}\left( t \right) \right)^{n2} \mathrm{inv}\left( K2^{n2} + \left( \mathrm{X5}\left( t \right) \right)^{n2} \right)] X2}\\
+\ce{ \varnothing &->[v3 \left( \mathrm{X3}\left( t \right) \right)^{n3} \mathrm{inv}\left( K3^{n3} + \left( \mathrm{X3}\left( t \right) \right)^{n3} \right)] X3}\\
+\ce{ \varnothing &->[v4 K4^{n4} \mathrm{inv}\left( K4^{n4} + \left( \mathrm{X1}\left( t \right) \right)^{n4} \right)] X4}\\
+\ce{ \varnothing &->[v5 \left( \mathrm{X2}\left( t \right) \right)^{n5} \mathrm{inv}\left( K5^{n5} + \left( \mathrm{X2}\left( t \right) \right)^{n5} \right)] X5}\\
 \ce{ X2 &<=>[k1][k2] X1 + 2 X4}\\
 \ce{ X4 &<=>[k3][k4] X3}\\
 \ce{ 3 X5 + X1 &<=>[k5][k6] X2}\\
@@ -54,15 +54,14 @@ raw"\begin{align}
 \end{align}
 ", "\r\n"=>"\n")
 
-
 # Latexify.@generate_test latexify(r, mathjax=false)
 @test latexify(r, mathjax = false) == replace(
 raw"\begin{align}
-\ce{ \varnothing &->[$\frac{v1 K1^{n1}}{K1^{n1} + \left( \mathrm{X2}\left( t \right) \right)^{n1}} \frac{v1 \left( \mathrm{X4}\left( t \right) \right)^{n1}}{K1^{n1} + \left( \mathrm{X4}\left( t \right) \right)^{n1}}$] X1}\\
-\ce{ \varnothing &->[$\frac{v2 \left( \mathrm{X5}\left( t \right) \right)^{n2}}{K2^{n2} + \left( \mathrm{X5}\left( t \right) \right)^{n2}}$] X2}\\
-\ce{ \varnothing &->[$\frac{v3 \left( \mathrm{X3}\left( t \right) \right)^{n3}}{K3^{n3} + \left( \mathrm{X3}\left( t \right) \right)^{n3}}$] X3}\\
-\ce{ \varnothing &->[$\frac{v4 K4^{n4}}{K4^{n4} + \left( \mathrm{X1}\left( t \right) \right)^{n4}}$] X4}\\
-\ce{ \varnothing &->[$\frac{v5 \left( \mathrm{X2}\left( t \right) \right)^{n5}}{K5^{n5} + \left( \mathrm{X2}\left( t \right) \right)^{n5}}$] X5}\\
+\ce{ \varnothing &->[$K1^{n1} v1^{2} \left( \mathrm{X4}\left( t \right) \right)^{n1} \mathrm{inv}\left( K1^{n1} + \left( \mathrm{X2}\left( t \right) \right)^{n1} \right) \mathrm{inv}\left( K1^{n1} + \left( \mathrm{X4}\left( t \right) \right)^{n1} \right)$] X1}\\
+\ce{ \varnothing &->[$v2 \left( \mathrm{X5}\left( t \right) \right)^{n2} \mathrm{inv}\left( K2^{n2} + \left( \mathrm{X5}\left( t \right) \right)^{n2} \right)$] X2}\\
+\ce{ \varnothing &->[$v3 \left( \mathrm{X3}\left( t \right) \right)^{n3} \mathrm{inv}\left( K3^{n3} + \left( \mathrm{X3}\left( t \right) \right)^{n3} \right)$] X3}\\
+\ce{ \varnothing &->[$v4 K4^{n4} \mathrm{inv}\left( K4^{n4} + \left( \mathrm{X1}\left( t \right) \right)^{n4} \right)$] X4}\\
+\ce{ \varnothing &->[$v5 \left( \mathrm{X2}\left( t \right) \right)^{n5} \mathrm{inv}\left( K5^{n5} + \left( \mathrm{X2}\left( t \right) \right)^{n5} \right)$] X5}\\
 \ce{ X2 &<=>[$k1$][$k2$] X1 + 2 X4}\\
 \ce{ X4 &<=>[$k3$][$k4$] X3}\\
 \ce{ 3 X5 + X1 &<=>[$k5$][$k6$] X2}\\
@@ -85,7 +84,7 @@ end p_a k n d_a p_b d_b r_a r_b
 @test latexify(r) == replace(
 raw"\begin{align}
 \require{mhchem}
-\ce{ \varnothing &<=>[\frac{p_{a} \left( B\left( t \right) \right)^{n}}{k^{n} + \left( B\left( t \right) \right)^{n}}][d_{a}] A}\\
+\ce{ \varnothing &<=>[p_{a} \left( B\left( t \right) \right)^{n} \mathrm{inv}\left( k^{n} + \left( B\left( t \right) \right)^{n} \right)][d_{a}] A}\\
 \ce{ \varnothing &<=>[p_{b}][d_{b}] B}\\
 \ce{ 3 B &<=>[r_{a}][r_{b}] A}
 \end{align}
@@ -94,8 +93,9 @@ raw"\begin{align}
 # Latexify.@generate_test latexify(r, mathjax=false)
 @test latexify(r, mathjax = false) == replace(
 raw"\begin{align}
-\ce{ \varnothing &<=>[$\frac{p_{a} \left( B\left( t \right) \right)^{n}}{k^{n} + \left( B\left( t \right) \right)^{n}}$][$d_{a}$] A}\\
+\ce{ \varnothing &<=>[$p_{a} \left( B\left( t \right) \right)^{n} \mathrm{inv}\left( k^{n} + \left( B\left( t \right) \right)^{n} \right)$][$d_{a}$] A}\\
 \ce{ \varnothing &<=>[$p_{b}$][$d_{b}$] B}\\
 \ce{ 3 B &<=>[$r_{a}$][$r_{b}$] A}
 \end{align}
 ", "\r\n"=>"\n")
+

test/make_model.jl

@@ -1,15 +1,12 @@
 ### Fetch required packages and reaction networks ###
 using DiffEqBase, Catalyst, Random, Test, UnPack
-using ModelingToolkit: operation, Sym, Term
+using ModelingToolkit: operation, Sym, istree
 include("test_networks.jl")
 
 
 ### Debugg functions ###
 
-function opname(x::Term)
-    nameof(operation(x))
-end
-opname(s::Sym) = nameof(s)
+opname(x) = istree(x) ? nameof(operation(x)) : nameof(x)
 
 alleq(xs,ys) = all(isequal(x,y) for (x, y) in zip(xs, ys))
 
@@ -178,9 +175,9 @@ for networks in identical_networks_1
         u0 = factor*rand(length(networks[1].states))
         p = factor*rand(length(networks[1].ps))
         t = rand()
-        @test all(abs.(f1(u0,p,t) .- f2(u0,p,t)) .< 100*eps())
-        @test all(abs.(f1.jac(u0,p,t) .- f2.jac(u0,p,t)) .< 100*eps())
-        @test all(abs.(g1(u0,p,t) .- g2(u0,p,t)) .< 100*eps())
+        @test all(abs.(f1(u0,p,t) .≈ f2(u0,p,t)))
+        @test all(abs.(f1.jac(u0,p,t) .≈ f2.jac(u0,p,t)))
+        @test all(abs.(g1(u0,p,t) .≈ g2(u0,p,t)))
     end
 end
 
@@ -233,9 +230,9 @@ for networks in identical_networks_2
         u0 = factor*rand(length(networks[1].states))
         p = factor*rand(length(networks[1].ps))
         t = rand()
-        @test all(abs.(f1(u0,p,t) .- f2(u0,p,t)) .< 100*eps())
-        @test all(abs.(f1.jac(u0,p,t) .- f2.jac(u0,p,t)) .< 100*eps())
-        @test all(abs.(g1(u0,p,t) .- g2(u0,p,t)) .< 100*eps())
+        @test all(f1(u0,p,t) .≈ f2(u0,p,t))
+        @test all(f1.jac(u0,p,t) .≈ f2.jac(u0,p,t))
+        @test all(g1(u0,p,t) .≈ g2(u0,p,t))
     end
 end
 
@@ -349,9 +346,9 @@ for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3]
     κ2 = factor*rand(); κ3 = factor*rand(); κ6 = factor*rand();
     τ = rand()
     p1 = [τ, κ2, κ3, τ, τ, κ6]; p2 = [κ2, κ3, κ6];
-    @test all(abs.(f1(u0,p1,τ) .- f2(u0,p2,τ)) .< 100*eps())
-    @test all(abs.(f1.jac(u0,p1,τ) .- f2.jac(u0,p2,τ)) .< 100*eps())
-    @test all(abs.(g1(u0,p1,τ) .- g2(u0,p2,τ)) .< 100*eps())
+    @test all(f1(u0,p1,τ) .≈ f2(u0,p2,τ))
+    @test all(f1.jac(u0,p1,τ) .≈ f2.jac(u0,p2,τ))
+    @test all(g1(u0,p1,τ) .≈ g2(u0,p2,τ))
 end
 
 

test/solve_SDEs.jl

@@ -95,15 +95,14 @@ for (i,networks) in enumerate(identical_networks)
         (i==2) && (u0[1] += 1000.)
         (i==3) ? (p[2:2:6] .*= 1000.; u0 .+= 1000) : (p[1] += 500.)
         prob1 = SDEProblem(networks[1],u0,(0.,100.),p)
-        sol1 = solve(prob1,ImplicitEM(),saveat=0.01,maxiters=1e7)
         prob2 = SDEProblem(networks[2][1],networks[2][2],u0,(0.,100.),p,noise_rate_prototype=networks[2][3])
-        sol2 = solve(prob2,ImplicitEM(),saveat=0.01,maxiters=1e7)
-        for i = 1:length(u0)
-            vals1 = getindex.(sol1.u[1000:end],i);
-            vals2 = getindex.(sol1.u[1000:end],i);
-            @test 0.8 < mean(vals1)/mean(vals2) < 1.25
-            @test 0.8 < std(vals1)/std(vals2) < 1.25
-        end
+        du1 = similar(u0); du2 = similar(u0);
+        prob1.f.f(du1,u0,p,0.0); prob2.f.f(du2,u0,p,0.0)
+        @test all(isapprox.(du1,du2))
+        g1 = zeros(numspecies(networks[1]),numreactions(networks[1]))
+        g2 = copy(g1)
+        prob1.f.g(g1,u0,p,0.0); prob2.f.g(g2,u0,p,0.0)
+        @test all(isapprox.(g1,g2))
     end
 end