AlgebraicPetri: the Petri-net view

petri_net returns a LabelledPetriNet alongside the name-indexed rate and initial-condition Dicts that AlgebraicPetri's own vectorfield needs (it indexes by species and transition name, not by position). It starts from the same Gamma(3, 1.5) lowering as the other backends — see Lowering a distribution to a dynamical system.

This backend's demo runs in its own environment: the package keeps one at test/algebraic_petri, and the demo below runs against it.

If you are reading this as a downloaded script rather than building the docs, this section needs a checked-out copy of the repository (it resolves paths relative to the package source).

using LoweredDistributions
using Distributions
using CairoMakie
using AlgebraOfGraphics
using DataFramesMeta

CairoMakie.activate!(type = "png", px_per_unit = 2)
set_theme!(theme_latexfonts(); fontsize = 14)

d = Gamma(3.0, 1.5)

petri_env = joinpath(pkgdir(LoweredDistributions), "test", "algebraic_petri")
petri_script = joinpath(pkgdir(LoweredDistributions), "docs",
    "algebraic_petri", "demo.jl")

"/home/runner/work/LoweredDistributions.jl/LoweredDistributions.jl/docs/algebraic_petri/demo.jl"

The demo script itself, printed from the file that is about to be run, so the page cannot drift from the code:

print(read(petri_script, String))

# Run by the getting-started tutorial in the isolated `test/algebraic_petri`
# environment (AlgebraicPetri 0.10 caps Catalyst at 13, so it cannot share an
# environment with this package's Catalyst 16 extension). Run it by hand with:
#
#   julia --project=test/algebraic_petri docs/algebraic_petri/demo.jl

using LoweredDistributions, Distributions, AlgebraicPetri

# The same delay the rest of the tutorial lowers: Gamma(3, 1.5), an exact
# three-compartment Erlang chain.
d = Gamma(3.0, 1.5)
built = petri_net(d)

println("places:      ", snames(built.petri_net))
println("transitions: ", tnames(built.petri_net))
println("rates:       ", sort(collect(built.rates), by = first))
println("u0:          ", sort(collect(built.u0), by = first))

# AlgebraicPetri's own mass-action vectorfield, evaluated at t = 0. All the
# mass starts in the first compartment, so it drains into the second at the
# chain's per-stage rate and nothing else has moved yet.
f! = vectorfield(built.petri_net)
du = Dict(k => 0.0 for k in keys(built.u0))
f!(du, built.u0, built.rates, 0.0)
println("du at t = 0: ", sort(collect(du), by = first))

# A trajectory through the vectorfield, so the calling tutorial can plot it
# alongside the other backends'. This environment has no ODE solver (only
# AlgebraicPetri and Distributions, kept minimal to dodge the Catalyst version
# clash — see the module header), so a plain fixed-step RK4 over the
# `vectorfield` closure is enough: the generator is linear, so RK4 at a small
# step is accurate well past plotting resolution.
species_names = snames(built.petri_net)

function deriv(u::Vector{Float64})
    udict = Dict(species_names[i] => u[i] for i in eachindex(species_names))
    dudict = Dict(n => 0.0 for n in species_names)
    f!(dudict, udict, built.rates, 0.0)
    return [dudict[n] for n in species_names]
end

function rk4_step(u::Vector{Float64}, dt::Float64)
    k1 = deriv(u)
    k2 = deriv(u .+ (dt / 2) .* k1)
    k3 = deriv(u .+ (dt / 2) .* k2)
    k4 = deriv(u .+ dt .* k3)
    return u .+ (dt / 6) .* (k1 .+ 2 .* k2 .+ 2 .* k3 .+ k4)
end

function print_trajectory(; dt::Float64 = 0.01, print_every::Int = 25)
    tspan = 0.0:dt:15.0                  # print_every * dt = 0.25 time units
    u = [built.u0[n] for n in species_names]
    for (i, t) in enumerate(tspan)
        if (i - 1) % print_every == 0
            for (n, v) in zip(species_names, u)
                println("TRAJ,", t, ",", n, ",", v)
            end
        end
        t == last(tspan) && break
        u = rk4_step(u, dt)
    end
    return nothing
end

print_trajectory()

And its real output, from a subprocess against that isolated environment. read throws on a non-zero exit, so a broken demo fails the docs build rather than silently printing nothing.

setup = `$(Base.julia_cmd()) --project=$petri_env -e "using Pkg; Pkg.instantiate()"`
demo = `$(Base.julia_cmd()) --project=$petri_env $petri_script`

demo_output = read(demo, String)
demo_lines = split(demo_output, '\n')
traj_lines = filter(l -> startswith(l, "TRAJ,"), demo_lines)
print(join(filter(l -> !startswith(l, "TRAJ,"), demo_lines), '\n'))

    Updating registry at `~/.julia/registries/General.toml`
   Installed XML2_jll ───────────── v2.15.3+0
   Installed JSON3 ──────────────── v1.14.3
   Installed NameResolution ─────── v0.1.5
   Installed Permutations ───────── v0.4.23
   Installed GeneralizedGenerated ─ v0.3.3
   Installed LabelledArrays ─────── v1.20.2
   Installed GATlab ─────────────── v0.2.4
   Installed Measures ───────────── v0.3.3
   Installed Combinatorics ──────── v1.1.0
   Installed Catlab ─────────────── v0.17.6
   Installed CompTime ───────────── v0.1.2
   Installed MLStyle ────────────── v0.4.17
   Installed AlgebraicInterfaces ── v0.1.4
   Installed StructEquality ─────── v2.1.0
   Installed AlgebraicPetri ─────── v0.10.0
   Installed JuliaVariables ─────── v0.2.4
   Installed DataStructures ─────── v0.18.22
   Installed StructTypes ────────── v1.11.0
   Installed Compose ────────────── v0.9.7
   Installed PrettyPrint ────────── v0.2.0
   Installed LightXML ───────────── v0.9.3
   Installed PEG ────────────────── v1.0.4
┌ Warning: tarball content does not match git-tree-sha1
└ @ Pkg.Operations /opt/hostedtoolcache/julia/1.12.6/x64/share/julia/stdlib/v1.12/Pkg/src/Operations.jl:814
     Cloning [227ef7b5-1206-438b-ac65-934d6da304b8] ACSets from https://github.com/AlgebraicJulia/ACSets.jl.git
   Installed ACSets ─────────────── v0.2.29
  Installing 1 artifacts
   Installed artifact XML2        2.2 MiB
    Updating `~/work/LoweredDistributions.jl/LoweredDistributions.jl/test/algebraic_petri/Project.toml`
  [4f99eebe] + AlgebraicPetri v0.10.0
  [31c24e10] + Distributions v0.25.130
  [8539b7eb] + LoweredDistributions v0.1.0 `../..`
  [f8b46487] + TestItemRunner v1.1.5
  [1c621080] + TestItems v1.0.0
  [8dfed614] ~ Test ⇒ v1.11.0
    Updating `~/work/LoweredDistributions.jl/LoweredDistributions.jl/test/algebraic_petri/Manifest.toml`
  [227ef7b5] + ACSets v0.2.29
  [7d9f7c33] + Accessors v0.1.45
  [79e6a3ab] + Adapt v4.7.0
  [23cfdc9f] + AlgebraicInterfaces v0.1.4
  [4f99eebe] + AlgebraicPetri v0.10.0
  [66dad0bd] + AliasTables v1.1.3
  [4fba245c] + ArrayInterface v7.28.1
  [134e5e36] + Catlab v0.17.6
  [d360d2e6] + ChainRulesCore v1.26.1
  [3da002f7] + ColorTypes v0.12.1
  [5ae59095] + Colors v0.13.1
  [861a8166] + Combinatorics v1.1.0
  [38540f10] + CommonSolve v0.2.12
  [bbf7d656] + CommonSubexpressions v0.3.1
  [0fb5dd42] + CompTime v0.1.2
  [34da2185] + Compat v4.18.1
  [a81c6b42] + Compose v0.9.7
  [a33af91c] + CompositionsBase v0.1.2
  [187b0558] + ConstructionBase v1.6.0
  [a8cc5b0e] + Crayons v4.2.0
  [9a962f9c] + DataAPI v1.16.0
⌅ [864edb3b] + DataStructures v0.18.22
  [e2d170a0] + DataValueInterfaces v1.0.0
  [163ba53b] + DiffResults v1.1.0
  [b552c78f] + DiffRules v1.16.0
  [31c24e10] + Distributions v0.25.130
  [ffbed154] + DocStringExtensions v0.9.5
  [e2ba6199] + ExprTools v0.1.11
  [1a297f60] + FillArrays v1.17.0
⌅ [53c48c17] + FixedPointNumbers v0.8.6
  [f6369f11] + ForwardDiff v1.4.4
  [f0ffcf3b] + GATlab v0.2.4
  [46192b85] + GPUArraysCore v0.2.0
  [a0844989] + Gamma v1.1.0
  [6b9d7cbe] + GeneralizedGenerated v0.3.3
  [34004b35] + HypergeometricFunctions v0.3.30
  [3587e190] + InverseFunctions v0.1.17
  [92d709cd] + IrrationalConstants v0.2.6
  [c8e1da08] + IterTools v1.10.0
  [82899510] + IteratorInterfaceExtensions v1.0.0
  [692b3bcd] + JLLWrappers v1.8.0
⌅ [682c06a0] + JSON v0.21.4
  [0f8b85d8] + JSON3 v1.14.3
  [b14d175d] + JuliaVariables v0.2.4
  [b964fa9f] + LaTeXStrings v1.4.0
  [2ee39098] + LabelledArrays v1.20.2
  [9c8b4983] + LightXML v0.9.3
  [2ab3a3ac] + LogExpFunctions v1.0.1
  [8539b7eb] + LoweredDistributions v0.1.0 `../..`
  [d8e11817] + MLStyle v0.4.17
  [1914dd2f] + MacroTools v0.5.16
  [442fdcdd] + Measures v0.3.3
  [e1d29d7a] + Missings v1.2.0
  [77ba4419] + NaNMath v1.1.4
  [71a1bf82] + NameResolution v0.1.5
⌅ [bac558e1] + OrderedCollections v1.8.2
  [90014a1f] + PDMats v0.11.41
  [12d937ae] + PEG v1.0.4
  [69de0a69] + Parsers v2.8.6
  [2ae35dd2] + Permutations v0.4.23
  [d236fae5] + PreallocationTools v1.4.0
  [aea7be01] + PrecompileTools v1.3.4
  [21216c6a] + Preferences v1.5.2
  [8162dcfd] + PrettyPrint v0.2.0
  [08abe8d2] + PrettyTables v3.4.4
  [43287f4e] + PtrArrays v1.4.0
  [1fd47b50] + QuadGK v2.11.3
  [3cdcf5f2] + RecipesBase v1.3.4
  [731186ca] + RecursiveArrayTools v4.3.5
  [189a3867] + Reexport v1.2.2
  [ae029012] + Requires v1.3.1
  [79098fc4] + Rmath v0.9.0
  [f2b01f46] + Roots v3.0.6
  [7e49a35a] + RuntimeGeneratedFunctions v0.5.24
  [431bcebd] + SciMLPublic v1.2.4
  [a2af1166] + SortingAlgorithms v1.2.3
  [276daf66] + SpecialFunctions v2.8.1
  [90137ffa] + StaticArrays v1.9.18
  [1e83bf80] + StaticArraysCore v1.4.4
  [10745b16] + Statistics v1.11.1
  [82ae8749] + StatsAPI v1.8.0
  [2913bbd2] + StatsBase v0.34.12
  [4c63d2b9] + StatsFuns v2.2.1
  [892a3eda] + StringManipulation v0.4.7
  [6ec83bb0] + StructEquality v2.1.0
  [856f2bd8] + StructTypes v1.11.0
  [2efcf032] + SymbolicIndexingInterface v0.3.51
  [3783bdb8] + TableTraits v1.0.1
  [bd369af6] + Tables v1.13.0
  [f8b46487] + TestItemRunner v1.1.5
  [1c621080] + TestItems v1.0.0
  [94ce4f54] + Libiconv_jll v1.18.0+0
  [efe28fd5] + OpenSpecFun_jll v0.5.6+0
  [f50d1b31] + Rmath_jll v0.5.2+0
  [02c8fc9c] + XML2_jll v2.15.3+0
  [0dad84c5] + ArgTools v1.1.2
  [56f22d72] + Artifacts v1.11.0
  [2a0f44e3] + Base64 v1.11.0
  [ade2ca70] + Dates v1.11.0
  [f43a241f] + Downloads v1.7.0
  [7b1f6079] + FileWatching v1.11.0
  [b77e0a4c] + InteractiveUtils v1.11.0
  [ac6e5ff7] + JuliaSyntaxHighlighting v1.12.0
  [b27032c2] + LibCURL v0.6.4
  [76f85450] + LibGit2 v1.11.0
  [8f399da3] + Libdl v1.11.0
  [37e2e46d] + LinearAlgebra v1.12.0
  [56ddb016] + Logging v1.11.0
  [d6f4376e] + Markdown v1.11.0
  [a63ad114] + Mmap v1.11.0
  [ca575930] + NetworkOptions v1.3.0
  [44cfe95a] + Pkg v1.12.1
  [de0858da] + Printf v1.11.0
  [3fa0cd96] + REPL v1.11.0
  [9a3f8284] + Random v1.11.0
  [ea8e919c] + SHA v0.7.0
  [9e88b42a] + Serialization v1.11.0
  [6462fe0b] + Sockets v1.11.0
  [2f01184e] + SparseArrays v1.12.0
  [f489334b] + StyledStrings v1.11.0
  [4607b0f0] + SuiteSparse
  [fa267f1f] + TOML v1.0.3
  [a4e569a6] + Tar v1.10.0
  [8dfed614] ~ Test ⇒ v1.11.0
  [cf7118a7] + UUIDs v1.11.0
  [4ec0a83e] + Unicode v1.11.0
  [e66e0078] + CompilerSupportLibraries_jll v1.3.0+1
  [deac9b47] + LibCURL_jll v8.15.0+0
  [e37daf67] + LibGit2_jll v1.9.0+0
  [29816b5a] + LibSSH2_jll v1.11.3+1
  [14a3606d] + MozillaCACerts_jll v2025.11.4
  [4536629a] + OpenBLAS_jll v0.3.29+0
  [05823500] + OpenLibm_jll v0.8.7+0
  [458c3c95] + OpenSSL_jll v3.5.4+0
  [bea87d4a] + SuiteSparse_jll v7.8.3+2
  [83775a58] + Zlib_jll v1.3.1+2
  [8e850b90] + libblastrampoline_jll v5.15.0+0
  [8e850ede] + nghttp2_jll v1.64.0+1
  [3f19e933] + p7zip_jll v17.7.0+0
        Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated -m`
Precompiling packages...
    644.9 ms  ✓ PEG
    661.8 ms  ✓ PrettyPrint
    802.2 ms  ✓ AlgebraicInterfaces
    859.3 ms  ✓ Measures
    579.5 ms  ✓ StructEquality
   1816.5 ms  ✓ Combinatorics
   1221.2 ms  ✓ XML2_jll
   2841.0 ms  ✓ StructTypes
   2759.8 ms  ✓ DataStructures
    895.8 ms  ✓ NameResolution
    890.9 ms  ✓ Permutations
   1010.9 ms  ✓ LightXML
   1117.2 ms  ✓ SortingAlgorithms
   5247.5 ms  ✓ LabelledArrays
   3509.0 ms  ✓ QuadGK
   3845.2 ms  ✓ Compose
   5075.6 ms  ✓ StatsBase
  11125.4 ms  ✓ JSON3
    865.4 ms  ✓ PDMats → StatsBaseExt
   6291.7 ms  ✓ Distributions
   3434.8 ms  ✓ Distributions → DistributionsTestExt
   3440.5 ms  ✓ LoweredDistributions
   4771.0 ms  ✓ Distributions → DistributionsChainRulesCoreExt
  25747.9 ms  ✓ TestItemRunner
  29983.4 ms  ✓ MLStyle
   1939.8 ms  ✓ JuliaVariables
   2337.7 ms  ✓ CompTime
   4297.8 ms  ✓ GeneralizedGenerated
  11227.7 ms  ✓ ACSets
  32790.6 ms  ✓ GATlab
  43959.9 ms  ✓ Catlab
  13820.3 ms  ✓ AlgebraicPetri
   2878.9 ms  ✓ LoweredDistributions → LoweredDistributionsAlgebraicPetriExt
  33 dependencies successfully precompiled in 124 seconds. 141 already precompiled.
places:      [:state1, :state2, :state3, :state4]
transitions: [:state1_state2, :state2_state3, :state3_state4]
rates:       [:state1_state2 => 0.6666666666666666, :state2_state3 => 0.6666666666666666, :state3_state4 => 0.6666666666666666]
u0:          [:state1 => 1.0, :state2 => 0.0, :state3 => 0.0, :state4 => 0.0]
du at t = 0: [:state1 => -0.6666666666666666, :state2 => 0.6666666666666666, :state3 => 0.0, :state4 => 0.0]

The transitions are the chain's interior hops plus its exit, the rates are all the same per-stage rate 1/1.5, and du at time zero shows the mass draining out of the first compartment into the second — the same generator the ODE and jump backends read, in Petri-net clothing.

The vectorfield reproduces the same trajectory

The demo script also integrates vectorfield(built.petri_net) forward (a plain fixed-step RK4, since the isolated environment carries no ODE solver of its own — see the demo source above) and prints the result as TRAJ rows, parsed back here into a DataFrame for plotting.

traj_df = reduce(vcat,
    let parts = split(line, ',')
        DataFrame(t = parse(Float64, parts[2]), place = parts[3],
            value = parse(Float64, parts[4]))
    end
    for line in traj_lines)

places_by_order = sort(unique(traj_df.place);
    by = s -> parse(Int, replace(s, "state" => "")))
absorbed = places_by_order[end]
traj_df.label = [p == absorbed ? "absorbed" : p for p in traj_df.place]

draw(
    data(traj_df) *
    mapping(:t, :value, color = :label) *
    visual(Lines, linewidth = 2);
    axis = (xlabel = "Time", ylabel = "Occupancy probability")
)

The absorbed place is the Petri net's own view of the CDF: plotting it against the exact Gamma(3, 1.5) CDF is the same faithfulness check the ODE page makes, reached through the mass-action vectorfield instead of the linear Kolmogorov ODE.

absorbed_df = @rsubset(traj_df, :place == absorbed)
cdf_df = vcat(
    DataFrame(t = absorbed_df.t, value = absorbed_df.value,
        kind = "Absorbed mass (Petri vectorfield)"),
    DataFrame(t = absorbed_df.t, value = cdf.(d, absorbed_df.t),
        kind = "Exact Gamma CDF")
)
draw(
    data(cdf_df) *
    mapping(:t, :value, color = :kind, linestyle = :kind) *
    visual(Lines, linewidth = 2);
    axis = (xlabel = "Time", ylabel = "Probability")
)