Use MTK's ODE func (#386)

* wip

* up

* use mtk

* wip

* pkg

* finally

* ubuntu-latest runner

.github/workflows/ci.yml

@@ -12,9 +12,8 @@ concurrency:
   cancel-in-progress: true
 
 env:
-  NBCONVERT_JOBS: '6'
-  LITERATE_PROC: '6'
-  JULIA_NUM_THREADS: '2'
+  NBCONVERT_JOBS: '4'
+  LITERATE_PROC: '4'
   OPENBLAS_NUM_THREADS: '1'
   ALLOWERRORS: 'false'
   NBCACHE: '.cache'
@@ -23,7 +22,7 @@ env:
 
 jobs:
   CI:
-    runs-on: self-hosted
+    runs-on: ubuntu-latest
     steps:
     - name: Checkout repository
       uses: actions/checkout@v4

Manifest.toml

@@ -595,9 +595,9 @@ version = "1.11.0"
 
 [[deps.Distributions]]
 deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "9d9e93d19c912ee6f0f3543af0d8839079dbd0d7"
+git-tree-sha1 = "4b138e4643b577ccf355377c2bc70fa975af25de"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.114"
+version = "0.25.115"
 
     [deps.Distributions.extensions]
     DistributionsChainRulesCoreExt = "ChainRulesCore"
@@ -850,9 +850,9 @@ version = "2.13.3+1"
 
 [[deps.FriBidi_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "1ed150b39aebcc805c26b93a8d0122c940f64ce2"
+git-tree-sha1 = "846f7026a9decf3679419122b49f8a1fdb48d2d5"
 uuid = "559328eb-81f9-559d-9380-de523a88c83c"
-version = "1.0.14+0"
+version = "1.0.16+0"
 
 [[deps.FunctionWrappers]]
 git-tree-sha1 = "d62485945ce5ae9c0c48f124a84998d755bae00e"
@@ -1209,9 +1209,9 @@ version = "1.11.0"
 
 [[deps.Libffi_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
+git-tree-sha1 = "27ecae93dd25ee0909666e6835051dd684cc035e"
 uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
-version = "3.2.2+1"
+version = "3.2.2+2"
 
 [[deps.Libgcrypt_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll"]
@@ -2652,9 +2652,9 @@ version = "1.3.243+0"
 
 [[deps.Wayland_jll]]
 deps = ["Artifacts", "EpollShim_jll", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
-git-tree-sha1 = "7558e29847e99bc3f04d6569e82d0f5c54460703"
+git-tree-sha1 = "85c7811eddec9e7f22615371c3cc81a504c508ee"
 uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
-version = "1.21.0+1"
+version = "1.21.0+2"
 
 [[deps.Wayland_protocols_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -2722,10 +2722,10 @@ uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
 version = "1.0.11+1"
 
 [[deps.Xorg_libXcursor_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
-git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
+git-tree-sha1 = "807c226eaf3651e7b2c468f687ac788291f9a89b"
 uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
-version = "1.2.0+4"
+version = "1.2.3+0"
 
 [[deps.Xorg_libXdmcp_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -2740,16 +2740,16 @@ uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
 version = "1.3.6+1"
 
 [[deps.Xorg_libXfixes_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
-git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll"]
+git-tree-sha1 = "6fcc21d5aea1a0b7cce6cab3e62246abd1949b86"
 uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
-version = "5.0.3+4"
+version = "6.0.0+0"
 
 [[deps.Xorg_libXi_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
-git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
+git-tree-sha1 = "984b313b049c89739075b8e2a94407076de17449"
 uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
-version = "1.7.10+4"
+version = "1.8.2+0"
 
 [[deps.Xorg_libXinerama_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXext_jll"]

docs/figs/ch4-01.jl

@@ -4,55 +4,46 @@
 Steady states and phase plots in an asymmetric network.
 ===#
 using DifferentialEquations
-using SimpleUnPack              ## @unpack
+using ModelingToolkit
 using Plots
 Plots.default(linewidth=2)
 
-# Convenience functions
-hil(x, k) = x / (x + k)
-hil(x, k, n) = hil(x^n, k^n)
+# The model for figure 4.1, 4.2, and 4.3.
+@independent_variables t
+@variables A(t) B(t)
+@parameters k1 k2 k3 k4 k5 n
+D = Differential(t)
 
-# The model and derivatives for figure 4.1, 4.2, and 4.3.
-
-function ∂A41(u, p, t)
-    a, b = u
-    @unpack k1, k2, k3, k4, k5, n = p
-    return da = k1 * hil(1, b, n) - (k3 + k5) * a
-end
-
-function ∂B41(u, p, t)
-    a, b = u
-    @unpack k1, k2, k3, k4, k5, n = p
-    return db = k2 + k5 * a - k4 * b
-end
-
-function model41!(D, u, p, t)
-    D[1] = ∂A41(u, p, t)
-    D[2] = ∂B41(u, p, t)
-    return nothing
-end
+eqs = [
+    D(A) ~ k1 / (1 + B^n) - (k3 + k5) * A,
+    D(B) ~ k2 + k5 * A - k4 * B
+]
+@mtkbuild osys = ODESystem(eqs, t)
 
 # ## Fig 4.2 A
+tend = 1.5
+ps1 = Dict(k1 => 20, k2 => 5, k3=> 5, k4 => 5, k5 => 2, n => 4)
+prob = ODEProblem(osys, [A => 0.0, B => 0.0], tend, ps1)
 
-ps1 = (k1=20., k2=5., k3=5., k4=5., k5=2., n=4.)
-u0s = (
+u0s = [
     [0.0, 0.0],
     [0.5, 0.6],
     [0.17, 1.1],
     [0.25, 1.9],
     [1.85, 1.7],
-)
-
-tend = 1.5
+]
 
 #---
-sols = [solve(ODEProblem(model41!, u0, tend, ps1)) for u0 in u0s];
+sols = map(u0s) do u0
+    solve(remake(prob, u0=u0))
+end
 
 #---
-plot(sols[1], xlabel="Time", ylabel="Concentration", title="Fig. 4.2 A (Time series)", labels=["[A]" "[B]"])
+plot(sols[1], xlabel="Time", ylabel="Concentration", title="Fig. 4.2 A (Time series)")
 
 # ## Fig. 4.2 B (Phase plot)
-plot( sols[1], idxs=(1, 2),
+plot(
+    sols[1], idxs=(A, B),
     xlabel="[A]", ylabel="[B]",
     aspect_ratio=:equal, legend=nothing,
     title="Fig. 4.2 B (Phase plot)",
@@ -74,23 +65,15 @@ plot!(fig, xlabel="Time", ylabel="Concentration")
 fig = plot(title="Fig. 4.3 B (Phase plot)")
 
 for sol in sols
-    plot!(fig, sol, idxs=(1, 2), legend=nothing)
+    plot!(fig, sol, idxs=(A, B), legend=nothing)
 end
 
-plot!(fig, xlabel="[A]", ylabel="[B]", xlims=(0., 2.), ylims=(0., 2.), aspect_ratio=:equal)
+plot!(fig, xlabel="[A]", ylabel="[B]", xlims=(0.0, 2.0), ylims=(0.0, 2.0), aspect_ratio=:equal)
 
 # Let's sketch vector fields in phase plots.
-∂A = function (x, y)
-    ∂A41((x, y), ps1, 0.0)
-end
-
-∂B = function (x, y)
-    ∂B41((x, y), ps1, 0.0)
-end
 
 ∂F44 = function (x, y; scale=20)
-    da = ∂A(x, y)
-    db = ∂B(x, y)
+    da, db = prob.f([x, y], prob.p, nothing)
     s = sqrt(hypot(da, db)) * scale
     return (da / s, db / s)
 end
@@ -106,7 +89,7 @@ fig = plot(title="Fig. 4.4 A (Phase plot with vector field)")
 quiver!(fig, xx, yy, quiver=∂F44, line=(:lightgrey), arrow=(:closed), aspect_ratio=:equal)
 
 for sol in sols
-    plot!(fig, sol, idxs=(1, 2), linealpha=0.7, legend=nothing)
+    plot!(fig, sol, idxs=(A, B), linealpha=0.7, legend=nothing)
 end
 
 plot!(fig, size=(600, 600), xlims=(rxy[1], rxy[end]), ylims=(rxy[1], rxy[end]), xlabel="[A]", ylabel="[B]")
@@ -117,25 +100,26 @@ fig = plot(title="Fig. 4.5 A (Phase plot with nullclines)")
 
 ## Phase plots
 for sol in sols
-    plot!(fig, sol, idxs=(1, 2), linealpha=0.7, lab=nothing)
+    plot!(fig, sol, idxs=(A, B), linealpha=0.7, lab=nothing)
 end
 
 ## nullclines
-contour!(fig, 0:0.01:2, 0:0.01:2, ∂A, levels=[0], cbar=false, line=(:black))
+∂A44 = (x, y) -> prob.f([x, y], prob.p, nothing)[1]
+∂B44 = (x, y) -> prob.f([x, y], prob.p, nothing)[2]
+contour!(fig, 0:0.01:2, 0:0.01:2, ∂A44, levels=[0], cbar=false, line=(:black))
 plot!(fig, Float64[], Float64[], line=(:black), label="A nullcline")  ## Adding a fake line for the legend of A nullcline
-contour!(fig, 0:0.01:2, 0:0.01:2, ∂B, levels=[0], cbar=false, line=(:black, :dash))
+contour!(fig, 0:0.01:2, 0:0.01:2, ∂B44, levels=[0], cbar=false, line=(:black, :dash))
 plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline") ## Adding a fake line for the legend of B nullcline
-plot!(fig, xlim=(0., 2.), ylim=(0., 2.), legend=:bottomright, size=(600, 600), xlabel="[A]", ylabel="[B]", aspect_ratio=:equal)
+plot!(fig, xlim=(0.0, 2.0), ylim=(0.0, 2.0), legend=:bottomright, size=(600, 600), xlabel="[A]", ylabel="[B]", aspect_ratio=:equal)
 
 # ## Figure 4.5 B
-
 # Vector field
 fig = plot(title="Fig. 4.5 B (Vector field with nullclines)")
 quiver!(fig, xx, yy, quiver=∂F44, line=(:lightgrey), arrow=(:closed), aspect_ratio=:equal)
 
 # Nullclines
-contour!(fig, 0:0.01:2, 0:0.01:2, ∂A, levels=[0], cbar=false, line=(:black))
+contour!(fig, 0:0.01:2, 0:0.01:2, ∂A44, levels=[0], cbar=false, line=(:black))
 plot!(fig, Float64[], Float64[], line=(:black), label="A nullcline")  ## Adding a fake line for the legend of A nullcline
-contour!(fig, 0:0.01:2, 0:0.01:2, ∂B, levels=[0], cbar=false, line=(:black, :dash))
+contour!(fig, 0:0.01:2, 0:0.01:2, ∂B44, levels=[0], cbar=false, line=(:black, :dash))
 plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline") ## Adding a fake line for the legend of B nullcline
-plot!(fig, xlim=(0., 2.), ylim=(0., 2.), legend=:bottomright, size=(600, 600), xlabel="[A]", ylabel="[B]")
+plot!(fig, xlim=(0.0, 2.0), ylim=(0.0, 2.0), legend=:bottomright, size=(600, 600), xlabel="[A]", ylabel="[B]")