Tutorial Pluto notebooks are in the MDCExamples repository.
At a terminal, in the directory of your choice, type the following commands:
git clone https://github.com/Dhruva2/MDCExamples.git cd MDCExamples julia
This clones the github repository with the examples, changes location to the cloned repository, and starts a julia session.
In your new julia session, open the Pluto environment:
using Pluto; Pluto.run()
This should open a new window in your web browser, in which you can run the notebooks.
Open and run the notebooks! Note that precompiling packages for the notebooks, the first time round, can take a few minutes. For information on how to run Pluto notebooks, see the docs
3b. Pluto notebooks run annotated .jl files. You can also run the .jl files themselves without using Pluto:
Suggested order to go through notebooks
No. 1 just shows you how to play around with cost functions
No. 2 is an example on an extremely simple model. Running curves takes fractions of a second here, so it's a good playground for playing with curve settings.
No. 3 is fairly detailed. It's useful if you want to see how scientific insight can be gained from looking at Minimally Disruptive Curves. Each curve on this notebook takes 2-10 mins to generate on my laptop (2017 macbook pro).
Nos 4. and 5. show you how to use the (much quicker to evaluate) collocation cost. I've generated, but haven't really had time to properly analyse, the minimally disruptive curves.
This example is less didactic than the notebooks, but useful for somebody who just wants to get things working quickly.
We do a cursory analysis of the Lotka-Volterra, predator-prey differential equation model. This is often used as a toy example for global sensitivity analysis tools, see e.g. the following web links: a, b, c. Because it's really simple.
You can copy and paste the code at the bottom, and run it as a .jl script. It runs in seconds (after package/function pre-compilation). Just make sure you add all the package dependencies first:
] add OrdinaryDiffEq, ForwardDiff, Statistics, Plots, LinearAlgebra, LaTeXStrings ] add https://github.com/Dhruva2/MinimallyDisruptiveCurves.jl.git
Note that one of the model features we analyse is non-differentiable. Mathematically. Computationally everything works fine, as we'll explain.
The dynamics for the system are:
and are the birth and death rates of 🐁
and are the birth and death rates of 🦉
... and here is a time-course of their populations when
p = [1.5,1.0,3.0,1.0]. We call these values the nominal parameters.
Nominal system (i.e. simulated with nominal parameter values)
We will follow b in having two features of interest:
mean population of 🐁 over time
max population of 🦉 over time
The model features of the nominal system are
Our cost function on model behaviour is the mean square deviation from these nominal model features.
By the way In your ecology class, you may have heard of selection. This is a (criticised) theory of how there are two extremal strategies that individual species trade off between:
r-Strategists (e.g. rabbits): high birth and death rate. low investment in individual organisms
K-strategists (e.g. humans): low birth and death rate. high investment in individual organisms
We will use this terminology in our description of the MD curves.
By running the code at the bottom of the page twice, setting
whichdir = 1and
whichdir=2respectively, we get two MD curves (see below).
As we move along the curve, predators move towards an -strategy, and prey moves towards a -strategy.
We can also move in the opposite direction (prey to strategy, predator to ). To avoid figure clutter, run the pasted code yourself to see this!
This curve instead finds a way to modulate the frequency of predator prey interactions, while preserving model features
Note how it abruptly changes direction when slowing frequency eventually prevents model features being maintained (as the period approaches the length of simulation). It goes from decreasing the frequency to increasing it.
Gives static plots of the minimally disruptive curves. Code for the gifs was a bit more involved!
using OrdinaryDiffEq, ForwardDiff, MinimallyDisruptiveCurves, Statistics, Plots, LinearAlgebra, LaTeXStrings #which md curve to plot which_dir = 2 ## define dynamics of differential equation function f(du,u,p,t) du = p*u - p*u*u #prey du = -p*u + p*u*u #predator end u0 = [1.0;1.0] # initial populations tspan = (0.0,10.0) #time span to simulate over t = collect(range(0, stop=10., length=200)) # time points to measure p = [1.5,1.0,3.0,1.0] # initial parameter values nom_prob = ODEProblem(f,u0,tspan,p) # package as an ODE problem nom_sol = solve(nom_prob, Tsit5()) # solve ## Model features of interest are mean prey population, and max predator population (over time) function features(p) prob = remake(nom_prob; p=p) sol = solve(prob, Tsit5(); saveat = t) return [mean(sol[1,:]), maximum(sol[2,:])] end nom_features = features(p) ## loss function, we can take as l2 difference of features vs nominal features function loss(p) prob = remake(nom_prob; p=p) p_features = features(p) loss = sum(abs2, p_features - nom_features) return loss end ## gradient of loss function function lossgrad(p,g) g[:] = ForwardDiff.gradient(p) do p loss(p) end return loss(p) end ## package the loss and gradient into a DiffCost structure cost = DiffCost(loss, lossgrad) """ We evaluate the hessian once only, at p. Why? to find locally insensitive directions of parameter perturbation The small eigenvalues of the Hessian are one easy way of defining these directions """ hess0 = ForwardDiff.hessian(loss,p) ev(i) = eigen(hess0).vectors[:,i] ## Now we set up a minimally disruptive curve, with nominal parameters p and initial direction ev(1) init_dir = ev(which_dir); momentum = 1.; span = (-15.,15.) curve_prob = MDCProblem(cost, p, init_dir, momentum, span) mdc = evolve(curve_prob, Tsit5) function sol_at_p(p) prob = remake(nom_prob; p=p) sol = solve(prob, Tsit5()) end p1 = plot(mdc; pnames = [L"p_1" L"p_2" L"p_3" L"p_4"]) cost_vec = [mdc.cost(el) for el in eachcol(trajectory(mdc))] p2 = plot(distances(mdc), log.(cost_vec), ylabel = "log(cost)", xlabel = "distance", title = "cost over MD curve"); mdc_plot = plot(p1,p2, layout=(2,1), size = (800,800)) nominal_trajectory = plot(sol_at_p(mdc(0.)[:states]), label = ["prey" "predator"]) perturbed_trajectory = plot(sol_at_p(mdc(-15.)[:states]), label = ["prey" "predator"]) traj_comparison = plot(nominal_trajectory, perturbed_trajectory, layout = (2,1), xlabel = "time", ylabel = "population")
MinimallyDisruptiveCurves.jl showed how we could continuously vary the parameters of the Lotka-Volterra model, while preserving
mean prey population over time
max predator population over time
One strategy (MDC1) was to alter the parameters to increase (decrease) birth and death rates of the predators while decreasing (increasing) those of the prey. Along the axis of -strategies to -strategies, this corresponds to one species moving in one direction, with the other going in reverse.
Another strategy (MDC2) was to modulate the frequency of the oscillations. MDC2 showed us exactly how to change the parameters to continuously modulate this frequency.
By the way
The maximum predator population over time is a non differentiable function of the solution. How come this worked?
First note that ReLu units in neural networks are also technically nondifferentiable at a point, but are still used in differentiable machine learning models. Same reason!
OK the details? The maximum of a collection of elements (predator measurements over time), is differentiable except at crossing points, where two elements share a maximal value. Different values of predator(t) will never have exactly the same value, due to numerical error as well as equation dynamics. So it works out! Indeed, notice the timepoints at which maximal values are attained changes with the model parameters.