Composing open dynamical systems I: Directed composition

By Sophie Libkind and James Fairbanks on January 4, 2021

There is a discrepancy between how scientists conceive of their models and how they are implemented in a computer program. Informally, scientists represent their models as a composition of many primitive interactions. For example, an ecologist studying an ecosystem with one hundred species examines interactions between pairs of species and takes the full ecosystem to be a composite of these primitive systems. However, when the scientist sits down at a computer to encode the model, this modular structure is lost. The ecologist must simply write down a 100-variable ODE.

In AlgebraicDynamics.jl (see also the documentation), we introduce a modeling framework which allows the user to encode a complex dynamical system as the composite of primitive dynamical systems. Following the mathematics of operads and operad algebras, we explicitly represent the composition syntax itself as an algebraic object. While traditional modeling tools use a fixed syntax that is provided by the tool, the operadic approach provides a level-shift, where syntax is elevated to become flexible and programmable. Users can define new syntaxes adapted to their domains yet still interoperate between different syntaxes using the infrastructure built for operads.

So what is this process that builds complex systems out of primitive ones? First, start with rules for how to combine primitive systems. This theory of composition is an algebraic structure that specifies a composition syntax. Then, choose a pattern which follows the rules of the syntax — called a composition pattern — and primitive systems to compose. Lastly, the oapply method returns the composite system. The meat of this post is to give examples this general approach.

The category theoretic parallel is:

We are interested in two distinct styles of composition (1) composition via directed communication and (2) composition via undirected communication.[1] We call dynamical systems that compose via directed communication machines and dynamical systems that compose via undirected communcation resource sharers. A zoo of composition theories is given below with directed theories in purple and undirected theories in blue.





In the following series of blog posts, we will explore the different species of this zoo and their implementations. In this post, we begin with directed theories. We choose to demonstrate the simpler theories to make the code readable and the examples constrained. These posts also showcase several categorical features of attributed \(\mathsf{C}\)-sets such as functorial data migration, \(\mathsf{C}\)-set transformations, and (co)limits of \(\mathsf{C}\)-sets.

[1] A third style of composition allows both directed and undirected communication. The composition theory for this style of communication is shown in orange above. We call dynamical systems that compose via both directed and undirect communcation resource sharing machines. (Libkind 2020).

Composition of machines

Machines are dynamical systems which compose via directed communication. Informally, a machine consists of five component parts:

We will assume throughout this blog post that the evolution rule is an ODE with exogenous variables (also called a driven or forced ODE). Exogenous variables drive the dynamics but their values are determined elsewhere. For example, in the differential equation

\[\dot r(t) = \alpha r(t) - \beta r(t)h(t),\]

\(h\) is an exogenous variable which drives the system, while \(\alpha\) and \(\beta\) are (fixed) parameters. Following the informal defintion, we implement a machine as a struct below.

struct Machine{T}
    ninputs::Int
    nstates::Int
    noutputs::Int
    dynamics::Function
    readout::Function
end

nstates(m::Machine)  = m.nstates
ninputs(m::Machine)  = m.ninputs
noutputs(m::Machine) = m.noutputs

Following the general approach, we will first define a composition syntax. Here is our first set of rules for composing machines: for each machine (locally called the receiver), each input of that machine is wired to the output of a machine (locally called the sender). The sender transmits its output to the receiver which uses the information as input to its evolution function. In the special case where the sender and receiver are the same machine, the composition induces a feedback loop.

This composition syntax is captured by the category of single-input port graphs, \(\mathsf{Th}(\mathsf{SIPortGraph})\), which is defined by the schema below.


A single-input port graph is a \(\mathsf{C}\)-set over \(\mathsf{Th}(\mathsf{SIPortGraph})\), i.e. a functor

\[d: \mathsf{Th}(\mathsf{SIPortGraph}) \to \mathsf{Set}.\]

This composition pattern consists of

  1. a set of boxes \(d(\textrm{Box})\).

  2. Every box \(b\) has a set of in-ports \(d(\operatorname{box_{in}})^{-1}(b) \subseteq d(P_{\textrm{in}})\) and a set of out-ports \(d(\operatorname{box_{out}})^{-1}(b) \subseteq d(P_{\textrm{out}})\).

  3. Every in-port \(p \in d(P_{\textrm{in}})\) is fed information by a unique out-port \(d(\operatorname{wire})(p) \in d(P_{\textrm{out}})\).

Implemented in Catlab, we have:

using Catlab, Catlab.CategoricalAlgebra

@present TheorySIPortGraph(FreeSchema) begin        
    Box::Ob
    InPort::Ob
    OutPort::Ob
    
    in_port_box::Hom(InPort, Box)
    out_port_box::Hom(OutPort, Box)
    wire::Hom(InPort, OutPort)
    end

const AbstractSIPortGraph = AbstractCSetType(TheorySIPortGraph)
const SIPortGraph = CSetType(TheorySIPortGraph, index=[:in_port_box, :out_port_box, :wire])

A machine may fill a box[2] of the composition pattern if

  1. the number of in-ports of the box equals the number of inputs to the machine, and

  2. the number of out-ports of the box equals the number of outputs of the machine.

The following method checks if a machine fills a chosen box of a composition pattern.

function fills(m::Machine, d::AbstractSIPortGraph, b::Int)
    b ≤ nparts(d, :Box) || error("Trying to fill box $b, when $d has fewer than $b boxes")
    
    return (ninputs(m) == length(incident(d, b, :in_port_box))) && (noutputs(m) == length(incident(d, b, :out_port_box)))
end

Now, given (1) a composition pattern and (2) a primitive machine filling each box, we can produce a new machine — the composite of the primitive machines — using the oapply method defined below. This composite machine has no inputs and no outputs. Its dynamics are induced by the driven dynamics of the primitive machines where the wires define how to set the inputs of each primitive machine.[3] We think of information of type T flowing along the wires.

using  Catlab.CategoricalAlgebra.FinSets

function oapply(d::SIPortGraph, xs::Vector{Machine{T}}) where T  
    
    for box in parts(d, :Box)
        fills(xs[box], d, box) || error("$(xs[box]) does not fill box $box")
    end
    
    States = coproduct((FinSet∘nstates).(xs))
    Outputs = coproduct((FinSet∘noutputs).(xs))
    
    function internal_readout(u)
        readouts = zeros(length(apex(Outputs)))
        for box in parts(d, :Box)
            view(readouts, legs(Outputs)[box](:)) .= xs[box].readout(view(u, legs(States)[box](:)))
        end
        return readouts
    end
    
    function v(u,p,t)
        dotu = zero(u)
        readouts = internal_readout(u)
        for box in parts(d, :Box)
            inputs = view(readouts, subpart(d, incident(d, box, :in_port_box), :wire))
            vars = legs(States)[box](:)
            view(dotu, vars) .= xs[box].dynamics(view(u, vars), inputs, t)
        end
        return dotu
    end
    
    return Machine{T}(0, length(apex(States)), 0, v, x -> 0)
end

We check that the machines have the correct signature for the box it fills when we construct the total system, not when we run it. This validation allows modeling software to reject malformed models and gives dynamical systems modelers the advantages associated with static type checking.

[2] What does "filling a box" correspond to in the operadic setting? Consider the operad algebra \(\mathsf{CDS}\) defined in (Schultz et al 2016). A box corresponds to a type \(t\) in the operad, and a machine can fill the box if it corresponds to an element of \(\mathsf{CDS}(t)\).
[3] The method oapply(d::SIPortGraph, xs::Vector{Machine{T}}) implements the operad algebra \(\mathsf{CDS}\) defined in (Schultz et al 2016) — see also (Vagner et al 2015) — restricted to the special case of morphisms where the codomain is the box with no incoming or outgoing wires. Such morphisms are equivalent to single-input port graphs. A generalization of the algebra \(\mathsf{CDS}\) to the operad of directed wiring diagrams is implemented in full in AlgebraicDynamics.jl.

Example: the Lotka-Volterra predator-prey model

A standard Lotka-Volterra predator-prey model is the composition of two machines:

  1. Evolution of a rabbit population — this machine takes one input which represents a population of predators, \(h\), that hunt rabbits. This machine has one output which emits the rabbit population \(r\). The dynamics of this machine is the driven ODE

\[\dot r = \alpha r - \beta r h.\]
  1. Evoluation of a fox population — this machine takes one input which represents a population of prey, \(e\), that are eaten by foxes. This machine has one output which emits the fox population \(f\). The dynamics of this machine is the driven ODE

\[\dot f =\gamma fe - \delta f .\]

Since foxes hunt rabbit, these machines compose by setting the fox population to be the input for rabbit evolution. Likewise, we set the rabbit population to be the input for fox evolution. We depict this composition — both the composition pattern and the primitive machines — as


and implement this composition in Julia as

α, β, γ, δ = 0.3, 0.015, 0.015, 0.7

dotr(x, p, t) = [α*x[1] - β*x[1]*p[1]]
dotf(x, p, t) = [γ*x[1]*p[1] - δ*x[1]]

# Define the primitive systems
rabbit = Machine{Float64}(1,1,1, dotr, x -> x)
fox    = Machine{Float64}(1,1,1, dotf, x -> x)

# Define the composition pattern
rabbitfox_pattern = SIPortGraph()
add_parts!(rabbitfox_pattern, :Box,     2)
add_parts!(rabbitfox_pattern, :OutPort, 2, out_port_box=[1,2])
add_parts!(rabbitfox_pattern, :InPort,  2, in_port_box=[1,2], wire=[2,1])

# Compose
rabbitfox_sys = oapply(rabbitfox_pattern, [rabbit, fox])

The machine rabbitfox_sys now represents the complete Lotka-Volterra model. We can approximate a trajectory of the ecosystem using the solver in DifferentialEquations.jl and plot the result.

Composition patterns from graphs

The innovation of composing dynamical systems as machines alleviates some of the bookkeeping endemic to implementing a complicated ODE. However, the mappings of ports can still be quite intricate in a composition of machines with many variables, boxes, and wires. In this section, we use directed graphs to present composition patterns more succinctly.[4]

Recall the schema for graphs introduced in an earlier blog post:

We can transform an instance of Graph into an instance of SIPortGraph via pullback functorial data migration (see Appendix) induced by the functor \(\mathsf{Th(SIPortGraph)} \to \mathsf{Th(Graph)}\) defined on objects by

and on morphisms by

Single-input port graphs induced by ordinary graphs are a restricted class of single-input port graphs. Generically, a single-input port graph allows out-port splitting, i.e. two distinct in-ports may be wired to a single out-port. This syntactic feature implements the copying of output data and is not available in the more restrictive composition syntax defined by the theory of graphs.

[4] This strategy is inspired by the composition of dynamical systems defined in (DeVille and Lerman 2015).

Example: Three species ecosystem

Consider a ocean ecosystem containing three species — little fish, big fish, and sharks — with two predation interactions — sharks eat big fish and big fish eat little fish.

In order to save ourselves the bookkeeping of ports, wires, and boxes, we can model the composition pattern as the graph ocean_graph and then migrate the data of ocean_graph to an instance of SIPortGraph.

using Catlab.Graphs.BasicGraphs: TheoryGraph, Graph

ocean_graph = Graph(3)
add_parts!(ocean_graph, :E, 4, src=[1,2,2,3], tgt=[2,1,3,2])
using Catlab.Theories: id
E = TheoryGraph[:E]

# Define the composition pattern via data migration
ocean_pattern = SIPortGraph()
migrate!(ocean_pattern, ocean_graph, 
    Dict(:InPort => :E, :OutPort => :E, :Box => :V),
    Dict(:wire => id(E), :in_port_box => :tgt, :out_port_box => :src))

show_html(ocean_pattern)
CSet with elements Box = 1:3, InPort = 1:4, OutPort = 1:4
InPort in_port_box wire
1 2 1
2 1 2
3 3 3
4 2 4
OutPort out_port_box
1 1
2 2
3 2
4 3

We can now construct a model of the total ocean ecosystem by applying the primitive machines corresponding to little fish, big fish, and shark evolution into the composition pattern. Again, using a standard ODE solver we can compute and graph a trajectory of this ecosystem.

α, β, γ, δ, β′, γ′, δ′ = 0.3, 0.015, 0.015, 0.7, 0.017, 0.017, 0.35

dotfish(f, p, t) = [α*f[1] - β*p[1]*f[1]]
dotFISH(F, p, t) = [γ*p[1]*F[1] - δ*F[1] - β′*p[2]*F[1]]
dotsharks(s, p, t) = [-δ′*s[1] + γ′*s[1]*p[1]]

# Define the primitive systems
fish   = Machine{Float64}(1,1,1, dotfish,   f->f)
FISH   = Machine{Float64}(2,1,2, dotFISH,   F->[F[1], F[1]])
sharks = Machine{Float64}(1,1,1, dotsharks, s->s)

# Compose
ocean_sys = oapply(ocean_pattern, [fish, FISH, sharks])

The bigger picture

In this post, we discussed two theories of directed composition — \(\mathsf{Th(SIPortGraph)}\) and \(\mathsf{Th(Graph)}\) — which fit into the larger zoo of directed composition syntaxes.



The functors between these theories define contravariant transformations between composition patterns in different syntaxes. The previous section exemplified this process applied to the functor from \(\mathsf{Th(SIPortGraph)}\) to \(\mathsf{Th(Graph)}\).

In the following subsections, we briefly describe the remaining theories of directed composition.

Port graphs and circular port graphs

Port graphs extend single-input port graphs by allowing in-ports to accept any number of incoming wires. In application, this means we can merge information (\(+\)) and generate trivial information (\(0\)). The schema for \(\mathsf{Th}(\mathsf{PortGraph})\) is below.



Concretely, every port graph consists of a set of boxes. Each box has set of in-ports and a set of out-ports. A port graph also has a set of wires which connect out-ports to in-ports.

Circular port graphs[5] similarly consist of a set of boxes and a set of wires which connect ports, but they do not distinguish between in-ports and out-ports — a machine may both receive and emit information through any port. The functor from \(\mathsf{Th(PortGraph)}\) to \(\mathsf{Th(CPG)}\) turns circular port graphs into port graphs by duplicating every port with one copy interpretted as an in-port and the other as an out-port. Circular port graphs are particularly useful for modeling systems where the composition pattern is given by a network, such as a mesh of points in a reaction-diffusion model or a network of cities in a transportation grid. Such networks often arise in physics modeling.

[5] A note on the nomenclature. We draw the boxes of port graphs as squares with ports on the top edge indicating the in-ports and the ports on the bottom edge indicating the out-ports. Since there is no such distinction between ports of a circular port graph, we draw the boxes of a circular port graph as circles. While port graphs are established terminology in the computer science literature, circular port graphs are not.

Open composition syntaxes

The theories discussed so far generate closed composition patterns because they produce composite machine that have no inputs and no outputs. The composite machine is closed in the sense that it cannot interact with other machines.

In AlgebraicDynamics.jl, we implement two open theories for directed composition (1) directed wiring diagrams and (2) open circular port graphs, along with corresponding oapply methods.

Up next

In the next blog post of this series, we will investigate undirected theories for composition of dynamical systems as well as the strength of open composition syntaxes.

Appendix: Pullback Functorial Data Migration

f Let \(F: \mathsf{C} \to \mathsf{D}\) be a functor between categories \(\mathsf{C}\) and \(\mathsf{D}\). Then \(F\) contravariantly induces a map from \(\mathsf{D}\)-sets to \(\mathsf{C}\)-sets by precomposition, which maps a \(\mathsf{D}\)-set

\[X: \mathsf{D} \to \mathsf{Set},\]

to the \(\mathsf{C}\)-set

\[F \operatorname{⨟} X : \mathsf{C} \to \mathsf{Set}.\]

We call this transformation pulling back \(X\) along \(F\). We saw an example of migrating the data of a half-edge graph to the data of a symmetric graph and vice versa[6] in a previous post.

[6] In that particular example, we were able to migrate the data in both directions because the functor \(F\) was invertible.

References

DeVille and Lerman, 2016. Dynamics on networks of manifolds. Symmetry, Integrability and Geometry: Methods and Applications. arXiv:1208.1513

Libkind, 2020. An algebra of resource sharing machines. arXiv:2007.14442.

Schultz, Spivak, and Vasilakopoulou, 2016. Dynamical systems and sheaves. arXiv:1609.08086.

Vagner, Spivak, and Lerman, 2015. Algebras of open dynamical systems on the operad of wiring diagrams. arXiv:1408.1598