Solutions to Exercises

Part 1

Follow the quick start instructions to install Julia, a code editor, and the git version control system.

If you’re stuck on setup, please ask an instructor or TA for help!

Get Julia to print “Hello World”.

Output text with the println function.

println("Hello World!")

Hello World!

String concatentation in Julia uses the * symbol.

println("Hello" * " " * "World!")

Hello World!

Print everything in a vector by using broadcasting. To broadcast, append . to your function, and pass it something you can iterate over (like a vector). Or, use map or foreach.

# Both of these return a vector of `nothing`s:

# println.(["Hello", "World!"])
# map(println, ["Hello", "World!"])

# This one produces side effects only:
foreach(print, ["Hello", " World!"])

Hello World!

Perform string interpolation using $.

greeting = "Hello"
greetee = "World"

println("$greeting, $(greetee)!")

Hello, World!

Install a Julia package, such as Catlab or AlgebraicDynamics, using the Julia package manager.

From a Julia REPL, enter the package manager by typing ]. Create a new project, which will organize the packages you’re using for this workshop, by typing pkg> activate . Next, add Catlab to your project by typing pkg> add Catlab You should see a bunch of stuff happening, including precompilation of the Catlab code.

When that’s finished, hit backspace to get out of the package manager and type

using Catlab

to load Catlab into your REPL session. To make sure it worked right, try something like

g = path_graph(Graph, 7)

to produce your first graph in Catlab!

Fork and clone the ACT2023Tutorial repository, make a change, commit it, and create a pull request.

Log into github.com
Navigate to github.com/AlgebraicJulia/ACT2023Tutorials
By the top right of the page there click the “Fork” button to the left of the “Star” button.
On the next page click “Create fork”, this will create a copy of the repository on your own GitHub account
Clone your fork (replacing <user_name> with your own GitHub username)
```
git clone https://github.com/<user_name>/ACT2023Tutorials.git
```
Navigate into the solutions folder and create a new directory with your name (replacing <your_name> with your own name)
```
cd ACT2023Tutorials/solutions
mkdir <your_name>
```

Put some code in there!

echo "print('Hello, World!')" > hello_world.jl

Stage, commit, and push your changes

git add hello_world.jl
git commit -m "Adding solutions for <your_name>"
git push

Navigate to your fork on GitHub (replacing <user_name> with your GitHub username): https://github.com/<user_name>/ACT2023Tutorials
Click the “Pull requests” tab across the top of the page and click the “New pull request”. This will automatically create a new pull request on the original repository from your fork
Click “Create pull request”
Put a short title and description and click “Create pull request”

Write a function that computes the product of two finite sets, along with the projection maps out.

For convenience, we’ll use the types FinSet and FinFunction exported by Catlab, which are similar to but more general than those in the problem statement.

using Catlab

function product_projections(A::FinSet{Int}, B::FinSet{Int})
  m, n = length(A), length(B)
  AB = FinSet(m*n)
  
  # In this calculation, `÷` (short for `div` function) is truncating integer
  # division and `%` (short for `mod` function) is the modulus.
  # Notice that we're using the skeleton of FinSet with sets `1:n` but the
  # calculation below is most natural for sets `0:(n-1)`.
  πA = FinFunction(i -> (i-1) ÷ n + 1, AB, A)
  πB = FinFunction(i -> (i-1) % n + 1, AB, B)
  
  (πA, πB)
end

An example:

using Test

π₁, π₂ = product_projections(FinSet(3), FinSet(4))
@test dom(π₁) == dom(π₂) == FinSet(12)
@test collect(π₁) == [1,1,1,1, 2,2,2,2, 3,3,3,3]
@test collect(π₂) == [1,2,3,4, 1,2,3,4, 1,2,3,4]

Test Passed

Implement a graph data structure in Catlab either (a) using Catlab’s support for C-sets (@acset_type) or (b) more directly, using finite sets (FinSet) and functions (FinFunction).

Graphs as C-sets:

using Catlab

@present SchGph(FreeSchema) begin
   V::Ob
   E::Ob
   src::Hom(E,V)
   tgt::Hom(E,V)
end

# Index the `src` and `tgt` functions for fast lookups of inverse images
# using the `incident` function.
@acset_type Gph(SchGph, index=[:src,:tgt])

# An example graph
g = @acset Gph begin
  V = 3
  E = 2
  src = [1,1]
  tgt = [2,1]
end

Gph {V:3, E:2}

E	src	tgt
1	1	2
2	1	1

We use the name Gph to avoid a name conflict with Catlab’s exported type Graph, but Graph is defined in Catlab in exactly the same way!

The exported Graph type includes support for visualization using Graphviz, which we can use to see what the above graph looks like:

g = @acset Graph begin
  V = 3
  E = 2
  src = [1,1]
  tgt = [2,1]
end

to_graphviz(g, node_labels=true)

Using finite sets and functions:

# Catlab exports these generic functions. We import them explicitly so that
# we can add new methods to them.
import Catlab: src, tgt, vertices, edges

struct MyGraph
  # We don't need fields for `E` and `V` since these are the domain and
  # codomain of `src` and `tgt`.
  src::FinFunction
  tgt::FinFunction

  function MyGraph(src::FinFunction, tgt::FinFunction)
    dom(src) == dom(tgt) || error("Source and target must have same domain")
    codom(src) == codom(tgt) || error("Source and target must have same codomain")
    new(src, tgt)
  end
end

# A few accessor functions. Our users shouldn't need to know the internal
# layout of the fields in the struct!
src(g::MyGraph) = g.src
tgt(g::MyGraph) = g.tgt
vertices(g::MyGraph) = codom(src(g)) # == codom(tgt(g))
edges(g::MyGraph) = dom(src(g)) # == dom(tgt(g))

Now let’s test that it works!

using Test

s, t = FinFunction([1,1], 3), FinFunction([2,1], 3)
g = MyGraph(s, t)

@test vertices(g) == FinSet(3)
@test edges(g) == FinSet(2)

Test Passed

While the direct implementation is rather straightforward, the acset-based implementation will be dramatically faster in certain situations, especially for graph algorithms that involve iterating over all edges with a given source or target–all such lists of edges are pre-computed by the indexing of src and tgt in the exported Graph type. It’s also much better-suited for efficiently generalizing graph algorithms to other types of graphs.

Write an optimized implementation of Euler’s method that computes the derivative in place and only makes allocations when initializing the algorithm.

For the sake of comparison, let’s first measure the allocations of the naive implementation of Euler’s method.

function euler(du::Function, init::Vector{Float64}, dt::Float64, m::Int)
  n = length(init)
  x = zeros(Float64, (m+1, n))
  x[1,:] = init
  for i in 1:m
    x[i+1,:] = x[i,:] + dt * du(x[i,:])
  end
  x
end

du(u) = [u[3], u[4], 0, -9.8]
@time trajectory = euler(du, [0, 0, 5.0, 10.0], 0.1, 1000);

  0.000209 seconds (5.00 k allocations: 500.266 KiB)

In the optimized version, we pass a function du! that computes the derivative in-place.

function inplace_euler(du!::Function, init::Vector{Float64}, dt::Float64, m::Int)
  n = length(init)
  x = zeros(Float64, (m+1, n))
  du = zeros(Float64, n)
  x[1,:] = init
  for i in 1:m
    du!(du, j -> x[i,j])
    for j in 1:n
      x[i+1,j] = x[i,j] + dt * du[j]
    end
  end
  x
end

function du!(du, u)
  du[1] = u(3)
  du[2] = u(4)
  du[3] = 0
  du[4] = -9.8
end
@time trajectory′ = inplace_euler(du!, [0, 0, 5.0, 10.0], 0.1, 1000);

  0.000027 seconds (5 allocations: 31.609 KiB)

Notice that the number of allocations is constant in the length of the trajectory. We better also check that we get the same answer as before.

@test trajectory′ ≈ trajectory

Test Passed

Part 2

Implement the dynamical system

\[ \dot{R} = \alpha R \]

using AlgebraicDynamics.

Define the dynamical system as a resource sharer:

using AlgebraicDynamics

ṙ(r,p,t) = p.α * r
rabbit_growth = ContinuousResourceSharer{Float64}(1, ṙ)

ContinuousResourceSharer(ℝ^1 → ℝ^1) with 1 exposed port

Simulate and plot:

using LabelledArrays
using DifferentialEquations, Plots

u0 = [1.0]
params = LVector(α=0.3)
tspan = (0.0, 10.0)

prob = ODEProblem(rabbit_growth, u0, tspan, params)
sol = solve(prob, Tsit5())

plot(sol, title="Exponential growth of rabbit population",
     legend=false, xlabel="time", ylabel="population size")

Write a call to the @relation macro that constructs the shown UWD.
Construct the UWD manually using the API for UWDs or the lower-level API for acsets.

Using the @relation macro:

uwd = @relation (x,y,z) where (w,x,y,z) begin
  R(x,w)
  S(y,w)
  T(z,w)
end

to_graphviz(uwd, box_labels=:name, junction_labels=:variable)

Using the UWD API:

d = RelationDiagram(3) # Create diagram with 3 outer ports
add_junctions!(d, 4, variable=[:w,:x,:y,:z]) # Add four junctions
set_junction!(d, 1:3, 2:4, outer=true) # Set junctions of outer ports

for name in [:R, :S, :T]
  add_box!(d, 2, name=name) # Add box with two ports
end
for (i, box) in enumerate(boxes(d))
  set_junction!(d, (box, 1), i+1)
  set_junction!(d, (box, 2), 1)
end

@test d == uwd

Test Passed

The advantange of using the programmatic API is that this code is easily generalized to a function that creates a “star-shaped UWD” of any size. You might try implementing this function as a further exercise.

Implement a dynamical system modeling the trajectory of a pendulum.

We’ll use the differential equation for the motion of a simple pendulum

\[ \frac{d^2 \theta}{dt^2} + \frac{g}{\ell} \sin\theta = 0, \]

where $g$ is the acceleration due to gravity and $\ell$ is the length of the rod. It’s fun that since we’re working numerically, we don’t need to bother with the small-angle approximation $\sin(\theta) \approx \theta$ that you probably remember from physics class.

We turn this second-order ODE into a first-order system of ODEs by introducing a new variable, the angular velocity $\omega$:

\[\begin{align*} \dot\theta &= \omega \\ \dot\omega &= -\frac{g}{\ell} \sin\theta \end{align*}\]

Now we can plug this into AlgebraicDynamics:

du(u,p,t) = [u[2], -p.g/p.ℓ * sin(u[1])]

pendulum = ContinuousResourceSharer{Float64}(2, du)

u0 = [π/4, 0]
params = LVector(g=9.8, ℓ=1)
tspan = (0.0, 10.0)

prob = ODEProblem(pendulum, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, title="Simple pendulum", xlabel="time", label=["θ" "ω"])

┌ Warning: Using arrays or dicts to store parameters of different types can hurt performance.
│ Consider using tuples instead.
└ @ SciMLBase ~/.julia/packages/SciMLBase/kTUaf/src/performance_warnings.jl:32

Implement a graph traversal algorithm, such as depth-first search or breadth-first search, using the acsets API (subpart, incident).

A stack-based depth-first search algorithm:

using Catlab
using DataStructures: Stack

""" Depth-first search in graph `g` starting at vertex `v`

The function `f` is called at every vertex in the search path.
"""
function dfs(f::Function, g::ACSet, v::Int)
  seen = Set{Int}()
  next = Stack{Int}()
  push!(next, v)
  while !isempty(next)
    v = pop!(next)
    v ∈ seen && continue
    f(v)
    push!(seen, v)
    # Add all outgoing edges to the stack.
    for e in incident(g, v, :src)
      # Syntactic sugar: `g[e, :tgt] == subpart(g, e, :tgt)``
      push!(next, g[e, :tgt])
    end
  end
end

Traversal of a directed cycle, starting at the third vertex, yields the vertices in the expected order.

g = cycle_graph(Graph, 6)
dfs(println, g, 3)

Use the union-find data structure to find the connected components of a graph, by making src(e) and tgt(e) be in the same equivalence class for each edge e.

using Catlab
using DataStructures: IntDisjointSets, union!, find_root

""" Find connected components of a graph.

Returns a dictionary from the component roots to lists of vertices in each
component.
"""
function connected_components(g::HasGraph)
  sets = IntDisjointSets(nparts(g, :V))
  for e in parts(g, :E)
    union!(sets, g[e,:src], g[e,:tgt])
  end

  components = Dict{Int,Vector{Int}}()
  for v in parts(g, :V)
    # The first argument inserts an empty array
    # when v's connected component has not yet been visited.
    component = get!(() -> Int[], components, find_root(sets, v))
    push!(component, v)
  end
  components
end

A graph with three components constructed as a coproduct of basic graphs:

g = path_graph(Graph, 3) ⊕ cycle_graph(Graph, 4) ⊕ star_graph(Graph, 4)

values(connected_components(g))

ValueIterator for a Dict{Int64, Vector{Int64}} with 3 entries. Values:
  [4, 5, 6, 7]
  [8, 9, 10, 11]
  [1, 2, 3]

Write a macro that creates a graph using a Graphviz-like syntax.

In Catlab, we typically use the pattern-matching features of MLStyle.jl to simplify the parsing of Julia expressions. For the sake of illustration, we deconstruct the expressions “by hand” below.

macro make_graph(expr)
  make_graph(expr)
end

function make_graph(block::Expr)
  g = Graph()
  vnames = Dict{Symbol,Int}()
  block.head == :block || error("Input to `make_graph` must be a block")
  for arg in block.args
    if arg isa LineNumberNode
      continue
    elseif arg isa Symbol
      haskey(vnames, arg) && error("Vertex $arg already defined")
      vnames[arg] = add_vertex!(g)
    elseif arg isa Expr && arg.head == :(->)
      sname, rhs = arg.args
      tname = rhs.args[2]
      add_edge!(g, vnames[sname], vnames[tname])
    else
      error("Cannot parse statement $arg")
    end
  end
  g
end

Let’s try it out on the given example:

g = @make_graph begin
  a
  b
  c
  a -> b
  b -> c
end

to_graphviz(g)

For comparison, Catlab ships with a similar macro @graph and a host of related but more complicated macros for building categories, diagrams in a category, and so on.

Part 3

Here is the schema for a DDS:

using Catlab
@present SchDDS(FreeSchema) begin
  X::Ob
  ϕ::Hom(X,X)
end

@acset_type DDS(SchDDS)

DDS

We can instantiate an empty DDS and then mutate this object to provide values for $X$ and $\phi$.

using Catlab

myDDS = DDS()
add_parts!(myDDS, :X, 5, ϕ=[3,3,4,1,2])
myDDS

DDS {X:5}

X	ϕ
1	3
2	3
3	4
4	1
5	2

We could alternatively make use of the @acset macro. This lets us specify an instance declaratively.

using Catlab

myDDS = @acset DDS begin
  X = 5
  ϕ = [3,3,4,1,2]
end
myDDS

DDS {X:5}

X	ϕ
1	3
2	3
3	4
4	1
5	2

Suppose that you have a doubly nested undirected wiring diagram of resource sharers. By this, I mean you have an undirected wiring diagram, where each box contains an undirected wiring diagram, and each box of this contains a resource sharer.

Then you can put together these resource sharers in two different ways. The first is by composing each of the inner undirected wiring diagrams, so that the outer wiring diagram now has a resource sharer in each box, and then composing the outer undirected wiring diagram. In other words, call oapply on each inner wiring diagram, and then call oapply on the other diagram.

The second way is to first collapse the nesting, by composing the two layers of wiring diagrams into one big wiring diagram, and then compose all of the resource sharers at one. In other words, call ocompose on the UWD of UWDs, and then call oapply on the remaining UWD with a resource sharer in each box.

Both methods of composition should ultimately result in the same composed resource sharer.

This is analogous to the law for a monad algebra:

Verify this rule by testing it out in code; construct a nested undirected wiring diagram, and then try out both ways of composing it.

This part of the AlgebraicDynamics code tests this out (more or less), with the caveat that they simplify things by only putting a sub-UWD in one box of the overall UWD, and then use a method of “ocompose” that just collapses a single nested box at a time, in accordance with the treatment of operads in which you compose two operations at a time, rather than $n+1$.

Any Petri net induces a dynamical system where the state space is $\mathbb{R}^S$ ($S$ is the set of species), and the vector field is given by mass action semantics (see Chapter 2: The rate equation). In AlgebraicPetri, we have a function called vectorfield which produces a function that computes the vector field for the mass action equations. Use this function (or for an extra challenge, reimplement this function) to turn an arbitrary Petri net into a resource sharer that exposes all of the populations of the species. Then pick a Petri net, turn it into a resource sharer, and compose it with another resource sharer of your choice.

Now think about the following question. You can compose Petri nets together by gluing their species together. If you compose Petri nets and then take the dynamical system given by mass action on the composed Petri net, is that the same dynamical system as turning each Petri net into a resource sharer and then composing the resource sharers? Prove or give a counterexample.

using AlgebraicPetri, AlgebraicDynamics.UWDDynam

function ContinuousResourceSharer{T}(pn::AbstractPetriNet) where {T}
  vf! = vectorfield(pn)
  vf(u, p, t) = vf!(zero(u), u, p, t)
  ContinuousResourceSharer{T}(ns(pn), vf)
end

sir_petri = PetriNet(LabelledPetriNet([:S,:I,:R], :inf=>((:S,:I)=>(:I,:I)), :rec=>(:I=>:R)))

sir_sharer = ContinuousResourceSharer{Float64}(sir_petri)

birth = ContinuousResourceSharer{Float64}(1, (u, p, t) -> sqrt(u[1]))

uwd = @relation (s,i,r) where (s,i,r) begin
  SIR(s,i,r)
  Birth(s)
end

sys = oapply(uwd, [sir_sharer, birth])

u0 = [1.0, 0.1, 0]
params = [0.5, 0.2]
tspan = (0.0, 10.0)

prob = ODEProblem(sys, u0, tspan, params)
sol = solve(prob, Tsit5())
plot(sol, title="SIR with square root growth", xlabel="time", label=["s" "i" "r"])

[ Info: Precompiling AlgebraicPetriOrdinaryDiffEqExt [b17ec253-c1e2-5e90-9bb1-27e3085d53bf]
WARNING: using UWDDynam.trajectory in module Main conflicts with an existing identifier.

Suppose that there are populations of rocks, populations of papers, and populations of scissors. Predation happens to each based on the classic pattern; i.e. scissors eat papers, papers eat rocks, and rocks eat scissors.

Write down a model for a pairwise interaction between two implements based on Lotka-Volterra dynamics, and then compose three copies of that model to produce a three-population system.

First, let’s define a resource sharing machine representing a predation process. This resource sharer has two state variables, one for the predator population and the other for the prey population. It is parameterized by a predation rate (i.e. how efficiently our predator inanimate objects eat the prey inanimate objects) and a growth rate (i.e. how efficiently our predators reproduce given how much they eat).

using AlgebraicDynamics.UWDDynam
using Catlab
using Plots

# Return a continuous resource sharer describing the interaction
# of a predator and prey based on predation rate and growth rate parameters.
function binary_predation_system(growth_rate, predation_rate)
    dynamics(u,_,_) = [
        # The predators can reproduce if there are at least 2 of them.
        u[1] >= 2 ? growth_rate*u[1]*u[2] : 0,
        # The predators cannot eat if there aren't at least 1 of them.
        u[1] > 1 ? -predation_rate*u[1]*u[2] : 0
    ]
    return ContinuousResourceSharer{Float64}(2, dynamics)
end

binary_predation_system (generic function with 1 method)

Now, we can use Catlab to define the composition pattern of our rock-paper-scissors system.

d = @relation (rocks, papers, scissors) begin
    crush(rocks, scissors)
    cover(papers, rocks)
    cut(scissors, papers)
end

to_graphviz(d, box_labels = :name, junction_labels = :variable, edge_attrs=Dict(:len => ".75"))

Now we can pick some rate parameters and fill our boxes with binary predation systems.

crush_rate = .2
cover_rate = .1
cut_rate = .3

rock_reproduction_rate = 0.1
paper_reproduction_rate = 0.2
scissor_reproduction_rate = 0.1

subsystems = Dict(
    :crush => binary_predation_system(rock_reproduction_rate, crush_rate),
    :cover => binary_predation_system(paper_reproduction_rate, cover_rate),
    :cut => binary_predation_system(scissor_reproduction_rate, cut_rate)
)

rps_system = oapply(d, subsystems)

ContinuousResourceSharer(ℝ^3 → ℝ^3) with 3 exposed ports

To see the results of our game, we can construct and solve an ODE problem for the system and plot the solution.

using OrdinaryDiffEq

u0 = [125.0, 100.0, 90.0]
tspan = (0.0,5.0)

prob = ODEProblem(rps_system, u0, tspan)
sol = solve(prob, Tsit5())

plot(sol, rps_system)

┌ Warning: To maintain consistency with solution indexing, keyword argument vars will be removed in a future version. Please use keyword argument idxs instead.
│   caller = ip:0x0
└ @ Core :-1

As we can see, this choice of parameters has led us to a steady state position where rock has won!

Figure out a good mathematical abstraction for, and implement nested acsets in Julia.

For some thoughts on what this might mean, see:

If I knew, I’d tell you.