{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# The Unique Features and Performance of DifferentialEquations.jl\n",
"\n",
" \n",
"### JuliaCon 2017\n",
" \n",
"### Chris Rackauckas \n",
"### University of California, Irvine"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Quick About Me \n",
"\n",
"- PhD Candidate in Mathematics at UC Irvine\n",
"- Research in:\n",
" - Numerical methods for stochastic ODEs/PDEs\n",
" - Multiscale and spatial biological simulations \n",
" - Effects of stochasticity in gene regulatory networks\n",
"- I created DiffEq because it was necessary:\n",
" - No software was using the latest methods\n",
" - I wanted the latest and fastest methods to benchmark against\n",
" - I am too lazy to change my code around: everything needs one interface\n",
" - I wanted to be able to use these methods to build large biological models\n",
" - I wanted parallel computing to be part of the design\n",
" - I wanted tools to be able to match models to data"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# DifferentialEquations.jl is about modeling\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Scientific laws, theories, and experiments describe how things change\n",
" \n",
"### Differential equations give a mechanistic model\n",
"\n",
" \n",
"\n",
"### Differential equations map \"change\" to simulations\n",
"\n",
" \n",
"\n",
"### Differential equations bridge from models to predictions\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Differential Equations vs Machine Learning\n",
"\n",
"- The models are different\n",
" - Machine learning uses general non-mechanistic models\n",
" - Differential equations utilize knowledge about the structure of the problem\n",
"- The reliance on data is different\n",
" - Machine learning needs to learn from data\n",
" - Machine learning models stack simple models together\n",
" - DE models are complex and nonlinear \n",
" - DE can model without data, or can calibrate to data\n",
" - DE have fewer parameters, and can be learned with less data\n",
"- Both are used to create predictions\n",
" - Ordinary and partial differential equations are the standard tools of physics\n",
" - Stochastic differential equations are a standard tool in stock market analysis"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## There are many types of differential equations \n",
"\n",
"- Ordinary differential equations evolve over one variable (let's call this time)\n",
"- Systems of differential equations evolve multiple quantities at the same time\n",
"- Stochastic differential equations add randomness (referred to as \"noise\")\n",
"- Delay differential equations allow for delayed interactions\n",
"- Gillespie simulations evolve discrete quantities via a stochastic model\n",
"- Partial differential equations evolve over multiple variables (space and time)\n",
"\n",
"DifferentialEquations.jl incorporates all of these under a single interface"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Problem 1: Ordinary Differential Equations\n",
"\n",
"http://docs.juliadiffeq.org/latest/tutorials/ode_example.html\n",
"\n",
"$$\\frac{du(t)}{dt} = f(u(t)p,t)$$\n",
"\n",
"where $u(t_0) = u_0$ is the known initial value. We want to know $u(t)$ for $t\\in[t_0,t_f]$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Exponential Growth Model \n",
"\n",
"- My money has been temporarily stolen by the banking cartel (it's in the bank)\n",
"- My money grows at a rate of 98% per year due to interest (really good bank)\n",
"- I start with $1: How much will I have in 5 years?\n",
"\n",
"$$ \\frac{du}{dt} = 0.98u $$\n",
"\n",
"$$u(0)=1$$\n",
"\n",
"Solve for $u(t)$ for $t\\in[0,5]$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"using DifferentialEquations\n",
"f(u,p,t) = p*u\n",
"u0=1.0\n",
"p = 0.98\n",
"tspan = (0.0,5.0)\n",
"prob = ODEProblem(f,u0,tspan,p)\n",
"sol = solve(prob);"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots; gr()\n",
"# Use the plot recipe, along with some plotting attributes\n",
"plot(sol,linewidth=5,title=\"My Bank Account\",\n",
" xaxis=\"Time (t)\",yaxis=\"Money!\",legend=false)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Controlling the Solvers\n",
"\n",
"See http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"11\n",
"15\n",
"[0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0]\n",
"[0.0,0.0634774,0.221501,0.43294,0.5,0.696266,0.923443,1.20696,1.53056,1.89898,2.30449,2.74501,3.21551,3.71248,4.23208,4.77109,5.0]\n",
"[0.0,5.0]\n"
]
}
],
"source": [
"sol = solve(prob)\n",
"println(length(sol.t))\n",
"sol = solve(prob,abstol=1e-6,reltol=1e-6)\n",
"println(length(sol.t))\n",
"sol = solve(prob,saveat=0.5,reltol=1e-6)\n",
"println(sol.t)\n",
"sol = solve(prob,tstops=[0.5],reltol=1e-4)\n",
"println(sol.t)\n",
"sol = solve(prob,reltol=1e-6,save_everystep=false)\n",
"println(sol.t)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Analyzing the Solution\n",
"\n",
"http://docs.juliadiffeq.org/latest/basics/solution.html"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.106052911995544\n",
"2.956279956560924\n",
"[2.0,2.30727,3.17657,4.63952,7.01861,11.2071,18.9435,34.2493,66.3366,137.407,273.572]\n"
]
}
],
"source": [
"sol = solve(prob,Tsit5())\n",
"println(sol.t[5]) # 5th timepoint\n",
"println(sol[5]) # 5th timepoint value\n",
"println([t+2u for (t,u) in zip(sol.t,sol.u)])"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"DifferentialEquations.jl comes with interpolations on the solution type"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/html": [
"1.5542610480525965"
],
"text/plain": [
"1.5542610480525965"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol(0.45) # The value of the solution at t=0.45"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Choosing Solvers\n",
"\n",
"http://docs.juliadiffeq.org/latest/solvers/ode_solve.html\n",
"\n",
"DifferentialEquations.jl is composed of many (30+?) packages which together provide its functionality under one interface. For Ordinary Differential Equations, we have:\n",
"\n",
"- OrdinaryDiffEq.jl: native Julia library. Fully featured, and in many cases faster than the classic C/Fortran libraries.\n",
"- Sundials.jl: provides Adams and BDF methods. A standard library\n",
"- LSODA.jl: provides the stability-detecting LSODA algorithm\n",
"- ODEInterface.jl: provides classic Fortran algorithms like `dopri` and `radau`\n",
"- ODE.jl: backwards compatibility\n",
"\n",
"DifferentialEquations.jl also provides automatic detection of the appropriate algorithm"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"# Automatically uses CVODE_BDF from Sundials\n",
"sol = solve(prob,alg_hints=[:stiff],reltol=1e-8,abstol=1e-8) \n",
"# Use Tsit5 from OrdinaryDiffEq\n",
"sol = solve(prob,Tsit5()); \n",
"# Use CVODE_BDF from Sundials\n",
"sol = solve(prob,CVODE_BDF()); \n",
"using ODEInterfaceDiffEq\n",
"# Use radau from ODEInterfaceDiffEq\n",
"sol = solve(prob,radau()); "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## The `*DiffEq` Libraries\n",
"\n",
"The `*DiffEq` libraries (OrdinaryDiffEq.jl, StochasticDiffEq.jl, DelayDiffEq.jl) are the native Julia solvers. They have a lot of unique features:\n",
"\n",
"- They routinely benchmark as the fastest implementations\n",
"- They have cheap/free high order interpolations (will show in a second!)\n",
"- They have expressive event handling\n",
"- They can handle generic `Number` and `AbstractArray` types\n",
"\n",
"The wrapped libraries, while classics (Sundials.jl, ODEInterfaceDiffEq.jl) are restricted to `Vector{Float64}`.\n",
"\n",
"http://docs.juliadiffeq.org/latest/basics/compatibility_chart.html"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Method Recommendations \n",
"\n",
"- `BS3()` for fast low accuracy non-stiff\n",
"- `Tsit5()` for standard non-stiff. This is the first algorithm to try in most cases.\n",
"- `Vern7()` for high accuracy non-stiff\n",
"- `Rosenbrock23()` for stiff equations with Julia-defined types, events, etc.\n",
"- `radau()` for stiff equations at higher accuracy on `Vector{Float64}`\n",
"- `CVODE_BDF()` for large stiff equations (discretizations of PDEs) on `Vector{Float64}`"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Translations from MATLAB/Python/R\n",
"\n",
"- `ode23` --> `BS3()`\n",
"- `ode45`/`dopri5` --> `DP5()`, though in most cases `Tsit5()` is more efficient\n",
"- `ode23s` --> `Rosenbrock23()`\n",
"- `ode113` --> `CVODE_Adams()`, though in many cases `Vern7()` or `Vern8()` are more efficient\n",
"- `dop853` --> `DP8()`, though in most cases `Vern7()` or `Vern8()` are more efficient\n",
"- `ode15s`/`vode` --> `CVODE_BDF()`, though in many cases `radau()` is more efficient\n",
"- `ode23t` --> `Trapezoid()`\n",
"- `lsoda` --> `lsoda()` (requires `Pkg.add(\"LSODA\"); using LSODA`)\n",
"- `ode15i` --> `IDA()`"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Example 2: Systems of Differential Equations\n",
"\n",
"$$\n",
"\\begin{align}\n",
"\\frac{dx}{dt} &= σ(y-x) \\\\\n",
"\\frac{dy}{dt} &= x(ρ-z) - y \\\\\n",
"\\frac{dz}{dt} &= xy - βz \\\\\n",
"\\end{align}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"function lorenz(du,u,p,t)\n",
" du[1] = 10.0(u[2]-u[1])\n",
" du[2] = u[1]*(28.0-u[3]) - u[2]\n",
" du[3] = u[1]*u[2] - (8/3)*u[3]\n",
"end\n",
"u0 = [1.0;0.0;0.0]\n",
"tspan = (0.0,25.0)\n",
"prob = ODEProblem(lorenz,u0,tspan)\n",
"sol = solve(prob,Tsit5());"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## The Plot Recipe Interpolates and Does Phase Plots"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sol,vars=(1,2,3))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"plot(sol,vars=(1,2,3),denseplot=false)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## ParameterizedFunctions\n",
"\n",
"$$\n",
"\\begin{align}\n",
"\\frac{dx}{dt} &= σ(y-x) \\\\\n",
"\\frac{dy}{dt} &= x(ρ-z) - y \\\\\n",
"\\frac{dz}{dt} &= xy - βz \\\\\n",
"\\end{align}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"(::LorenzExample) (generic function with 15 methods)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g = @ode_def LorenzExample begin\n",
" dx = σ*(y-x)\n",
" dy = x*(ρ-z) - y\n",
" dz = x*y - β*z\n",
"end σ ρ β\n",
"p = [10.0,28.0,8/3]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## DifferentialEquations.jl respects the user's input type\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using StaticArrays, DifferentialEquations\n",
"A = big.(@SMatrix [ 1.0 0.0 0.0 -5.0\n",
" 4.0 -2.0 4.0 -3.0\n",
" -4.0 0.0 0.0 1.0\n",
" 5.0 -2.0 2.0 3.0])\n",
"u0 = big.(@SMatrix rand(4,2))\n",
"tspan = (0.0,1.0); f_SA(u,p,t) = A*u\n",
"prob = ODEProblem(f_SA,u0,tspan); \n",
"sol = solve(prob,Tsit5()); plot(sol)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.0 s\n",
"1.1813604129914028 N\n"
]
},
{
"data": {
"text/html": [
"1.0869040537266441 N"
],
"text/plain": [
"1.0869040537266441 N"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Unitful\n",
"f_units(y,p,t) = 0.5*y / 3.0u\"s\"\n",
"u0 = 1.0u\"N\"\n",
"prob = ODEProblem(f_units,u0,(0.0u\"s\",1.0u\"s\"))\n",
"sol =solve(prob,Tsit5());\n",
"println(sol.t[end])\n",
"println(sol[end])\n",
"sol(0.5u\"s\")"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"14.778112197857341\n",
"14.7781121978613\n",
"22.167168296786013\n",
"22.16716829679195\n"
]
}
],
"source": [
"using DifferentialEquations, SymEngine\n",
"y0 = symbols(:y0)\n",
"u0 = y0\n",
"f_sym(u,p,t) = 2u\n",
"prob = ODEProblem(f_sym,u0,(0.0,1.0))\n",
"sol = solve(prob,RK4(),dt=1/10);\n",
"\n",
"sol = solve(prob,RK4(),dt=1/1000)\n",
"end_solution = lambdify(sol[end])\n",
"\n",
"println(end_solution(2)) # 14.778112197857341\n",
"println(2exp(2)) # 14.7781121978613\n",
"println(end_solution(3)) # 22.167168296786013\n",
"println(3exp(2)) # 22.16716829679195"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## DifferentialEquations.jl is unique when speed matters\n",
"\n",
"DifferentialEquations.jl is a research lab which includes many new methods"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"using DifferentialEquations, Plots, ODEInterfaceDiffEq, ODE\n",
"# 2D Linear ODE\n",
"function test_f(du,u,p,t)\n",
" for i in 1:length(u)\n",
" du[i] = 1.01*u[i]\n",
" end\n",
"end\n",
"function test_f(::Type{Val{:analytic}},u₀,p,t)\n",
" u₀*exp(1.01*t)\n",
"end\n",
"tspan = (0.0,10.0)\n",
"prob = ODEProblem(test_f,rand(100,100),tspan)\n",
"\n",
"abstols = 1./10.^(3:13)\n",
"reltols = 1./10.^(0:10);\n",
"setups = [Dict(:alg=>DP5())\n",
" Dict(:alg=>ode45())\n",
" Dict(:alg=>dopri5())\n",
" Dict(:alg=>Tsit5())]\n",
"names = [\"DifferentialEquations\";\"ODE\";\"ODEInterface\";\"DifferentialEquations Tsit5\"]\n",
"wp = ode_workprecision_set(prob,abstols,reltols,setups;names=names,save_everystep=false);"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(wp)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Callbacks and Event Handling\n",
"\n",
"You can directly apply dynamic behavior to an ODE through its callbacks and event handling"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"function b_condition(u,t,integrator) # Event when event_f(t,u) == 0\n",
" u[1]\n",
"end\n",
"\n",
"function b_affect!(integrator)\n",
" integrator.u[2] = -integrator.u[2]\n",
"end\n",
"\n",
"cb = ContinuousCallback(b_condition,b_affect!)\n",
"\n",
"bb = @ode_def_bare BallBounce begin\n",
" dy = v\n",
" dv = -g\n",
"end g=9.81\n",
"\n",
"u0 = [50.0,0.0]\n",
"tspan = (0.0,15.0)\n",
"prob = ODEProblem(bb,u0,tspan)\n",
"sol = solve(prob,Tsit5(),callback=cb);"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sol,plotdensity=1000)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Pretty much anything can be changed with the integrator\n",
"\n",
"Let's model a growing cell population\n",
"\n",
"- Each cell has an internal protein `X`\n",
"- When the concentration of `X` hits 1, the cell divides\n",
"- When there are more cells, there is less food and so `X` is created less rapidly"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"using DifferentialEquations, Plots\n",
"function growth(du,u,p,t)\n",
" for i in 1:length(u)\n",
" du[i] = 0.3u[i]/length(u)\n",
" end\n",
"end\n",
"\n",
"function condition(u,t,integrator) # Event when event_f(t,u) == 0\n",
" 1-maximum(u)\n",
"end\n",
"\n",
"function affect!(integrator)\n",
" u = integrator.u\n",
" resize!(integrator,length(u)+1)\n",
" maxidx = findmax(u)[2]\n",
" Θ = rand()\n",
" u[maxidx] = Θ\n",
" u[end] = 1-Θ\n",
" length(u) == 10 && terminate!(integrator)\n",
" nothing\n",
"end\n",
"\n",
"callback = ContinuousCallback(condition,affect!)\n",
"u0 = [0.2]\n",
"tspan = (0.0,100.0)\n",
"prob = ODEProblem(growth,u0,tspan)\n",
"sol = solve(prob,callback=callback);"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sol.t,map((x)->length(x),sol[:]),lw=3,\n",
" ylabel=\"Number of Cells\",xlabel=\"Time\")"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ts = linspace(0,sol.t[end],1000)\n",
"plot(ts,map((x)->x[1],sol.(ts)),lw=3,\n",
" ylabel=\"Amount of X in Cell 1\",xlabel=\"Time\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## This flexibility takes DiffEq from equation solver to simulation engine\n",
"\n",
"- Abstract types can be used to encode complex models\n",
"- Integrators and event handling gives dynamic interactions\n",
"- Speed, control, adaptivity, plotting, etc. then all comes for free!\n",
"\n",
"This leads to more sustainable scientific software development:\n",
"\n",
"- Many scientific projects develop their own internal solvers for their models\n",
" - This requires a lot of extra time!\n",
" - How well tested is it? How well optimized is it? How maintainable is it?\n",
" - Is it using all of the latest PI-stepsize adaptivity? Verner's coefficients? etc.?\n",
"- DifferentialEquations.jl is thus ideal as a core engine in complex models\n",
" - Developer-friendly dependency management, MIT licensed\n",
" - Active developer channels (Gitter, or anywhere...)\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Example: Multiscale Biological Models\n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"using MultiScaleArrays\n",
"using OrdinaryDiffEq, DiffEqBase, StochasticDiffEq\n",
"\n",
"# Define a hierarchy\n",
"\n",
"immutable Cell{B} <: AbstractMultiScaleArrayLeaf{B}\n",
" x::Vector{B}\n",
"end\n",
"immutable Population{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}\n",
" x::Vector{T}\n",
" y::Vector{B}\n",
" end_idxs::Vector{Int}\n",
"end\n",
"immutable Tissue{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}\n",
" x::Vector{T}\n",
" y::Vector{B}\n",
" end_idxs::Vector{Int}\n",
"end\n",
"immutable Embryo{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}\n",
" x::Vector{T}\n",
" y::Vector{B}\n",
" end_idxs::Vector{Int}\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"cell1 = Cell([1.0;2.0;3.0])\n",
"cell2 = Cell([4.0;5])\n",
"population = construct(Population,deepcopy([cell1,cell2])) # Make a Population from cells\n",
"cell3 = Cell([3.0;2.0;5.0])\n",
"cell4 = Cell([4.0;6])\n",
"population2 = construct(Population,deepcopy([cell3,cell4]))\n",
"tissue1 = construct(Tissue,deepcopy([population,population2])) # Make a Tissue from Populations\n",
"tissue2 = construct(Tissue,deepcopy([population2,population]))\n",
"embryo = construct(Embryo,deepcopy([tissue1,tissue2])); # Make an embryo from Tissues"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Embryo{Tissue{Population{Cell{Float64},Float64},Float64},Float64}(Tissue{Population{Cell{Float64},Float64},Float64}[Tissue{Population{Cell{Float64},Float64},Float64}(Population{Cell{Float64},Float64}[Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.290911,0.581821,0.872732]),Cell{Float64}([1.16364,1.45455]),Cell{Float64}([0.0,0.0,0.0])],Float64[],[3,5,8]),Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.600023,0.400015,1.00004]),Cell{Float64}([0.80003,1.20005])],Float64[],[3,5])],Float64[],[8,13]),Tissue{Population{Cell{Float64},Float64},Float64}(Population{Cell{Float64},Float64}[Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.600023,0.400015,1.00004]),Cell{Float64}([0.80003,1.20005])],Float64[],[3,5]),Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.250004,0.500008,0.750013]),Cell{Float64}([1.00002,1.25002])],Float64[],[3,5])],Float64[],[5,10])],Float64[],[13,23])"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function cell_ode(dcell,cell,p,t)\n",
" m = mean(cell)\n",
" for i in eachindex(cell)\n",
" dcell[i] = -m*cell[i]\n",
" end\n",
"end\n",
"function embryo_ode(dembryo,embryo,p,t)\n",
" for (cell,y,z) in LevelIterIdx(embryo,2)\n",
" cell_ode(@view dembryo[y:z],cell,p,t)\n",
" end\n",
"end\n",
"\n",
"# Have a cell divide at t=0.5\n",
"tstop = [0.5]\n",
"function grow_condition(u,t,integrator)\n",
" t ∈ tstop\n",
"end\n",
"function grow_affect!(integrator)\n",
" add_daughter!(integrator,integrator.u[1,1,1],1,1)\n",
"end\n",
"growing_cb = DiscreteCallback(grow_condition,grow_affect!)\n",
"\n",
"prob = ODEProblem(embryo_ode,embryo,(0.0,1.0));\n",
"sol = solve(prob,Tsit5(),callback=growing_cb,tstops=tstop)\n",
"\n",
"sol[end]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# That's the simulation of deterministic models\n",
"\n",
"## Let's add stochasticity"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Stochastic Differential Equations\n",
"\n",
"$$ du = f(t,u)dt + \\sum_i g_i(t,u)dW^i $$\n",
"\n",
"- The $f$ term is the deterministic or \"drift\" term\n",
"- The $g$ term is the noise or \"diffusion term\n",
"\n",
"Well-known model: the Black-Scholes Model \n",
"\n",
"- My stock is growing at an average rate of 101%\n",
"- However, the value of my money has fluctuatations of 87% of its current value\n",
"\n",
"$$ du = 1.01u + 0.87dW_t $$"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using DifferentialEquations, Plots; gr()\n",
"α=1.01; β=0.87; u₀=1.0\n",
"bs_f(u,p,t) = α*u\n",
"bs_g(u,p,t) = β*u\n",
"bs_f(::Type{Val{:analytic}},u₀,p,t,W) = u₀*exp((α-(β^2)/2)*t+β*W)\n",
"prob = SDEProblem(bs_f,bs_g,u₀,(0.0,1.0))\n",
"sol = solve(prob,EM(),dt=1/2^(4))\n",
"plot(sol,plot_analytic=true)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Fast adaptive methods are unique to DifferentialEquations.jl\n"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol =solve(prob,SRIW1())\n",
"plot(sol,plot_analytic=true)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## They tend to be orders of magnitude faster than standard methods\n",
"\n",
"Many models cannot be solved without these new methods. Even those that can are faster with them.\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Gillespie Simulations\n",
"\n",
"Gillespie simulations are discrete models which change stochastically over time.\n",
"\n",
"- A protein binds to DNA at a rate X\n",
"- A protein is created at a rate Y\n",
"\n",
"These types of models are related to differential equations\n",
"\n",
"- When the number of objects gets high, it converges for an SDE on the concentrations\n",
"- When the number gets really high, it converges to a deterministic ODE"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2} with problem DiffEqBase.DiscreteProblem{uType,tType,isinplace,F,C} and aggregator DiffEqJump.Direct\n",
"Number of constant rate jumps: 2\n",
"Number of variable rate jumps: 0\n"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"prob = DiscreteProblem([100,1,0],(0.0,100.0))\n",
"\n",
"rate = (u,p,t) -> (0.1/100.0)*u[1]*u[2]\n",
"jump_affect! = function (integrator)\n",
" integrator.u[1] -= 1\n",
" integrator.u[2] += 1\n",
"end\n",
"jump = ConstantRateJump(rate,jump_affect!)\n",
"\n",
"rate = (u,p,t) -> 0.01u[2]\n",
"jump_affect! = function (integrator)\n",
" integrator.u[2] -= 1\n",
" integrator.u[3] += 1\n",
"end\n",
"jump2 = ConstantRateJump(rate,jump_affect!)\n",
"jump_prob = JumpProblem(prob,Direct(),jump,jump2)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol = solve(jump_prob,Discrete(apply_map=false));\n",
"plot(sol)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Mix jumps and differential equations"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"diff_f = function (du,u,p,t)\n",
" du[4] = u[2]*u[3]/1000 - u[1]*u[2]/100000\n",
"end\n",
"diff_g = function (du,u,p,t)\n",
" du[4] = 0.1u[4]\n",
"end\n",
"\n",
"prob = SDEProblem(diff_f,diff_g,[100.0,1.0,0.0,1.0],(0.0,100.0))\n",
"jump_prob = JumpProblem(prob,Direct(),jump,jump2)\n",
"sol = solve(jump_prob,SRIW1(),abstol=1e-1,reltol=1e-1); plot(sol)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"
"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## These features together allow for expressive scientific modeling\n",
"\n",
"Let's go back to the multiscale embyro model:\n",
"\n",
"- We can have cell proteins change according to a stochastic differential equation\n",
"- A rate of division can be dependent on a specific \"growth factor\"\n",
"- \"Jumps\" using this rate can cause cell division, growing our multiscale model\n",
"- We get high performance algorithms and optimized implementations the whole way down.\n",
"\n",
"<3 Julia"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Addons\n",
"\n",
"Up until now, I have shown how to solve many different equations and encode a wide variety of models to simulate using the various solvers in DifferentialEquations.jl. But because all of this uses a single interface, an addon suite was able to be developed.\n",
"\n",
"- Tools for automatic parallelization of Monte Carlo simulations\n",
"- Uncertainty quantification\n",
"- Parameter estimation and data fitting"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Monte Carlo Simulations\n",
"\n",
"Many times you need to solve the same equation multiple times. This is especially true for stochastic simulations."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"pf = function (du,u,p,t)\n",
" du[1] = p[1] * u[1] - p[2] * u[1]*u[2]\n",
" du[2] = -3 * u[2] + u[1]*u[2]\n",
"end\n",
"pf = ParameterizedFunction(pf_func,)\n",
"pg = function (du,u,p,t)\n",
" du[1] = p[3]*u[1]\n",
" du[2] = p[4]*u[2]\n",
"end\n",
"p = [1.5,1.0,0.1,0.1]\n",
"prob = SDEProblem(pf,pg,[1.0,1.0],(0.0,10.0),p)\n",
"\n",
"prob_func = function (prob,i)\n",
" problem_new_parameters(prob,[1.5,1.0,0.3rand(),0.3rand()])\n",
" prob\n",
"end\n",
"monte_prob = MonteCarloProblem(prob,prob_func=prob_func)\n",
"sim = solve(monte_prob,SRIW1(),num_monte=10);"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(sim,linealpha=0.6,color=:blue,vars=(0,1),title=\"The Timeseries'\")\n",
"plot!(sim,linealpha=0.6,color=:red,vars=(0,2))"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"summ = MonteCarloSummary(sim,0:0.1:10)\n",
"gr() # Note that plotly does not support ribbon plots\n",
"plot(summ,fillalpha=0.4)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Uncertainty Quantification\n",
"\n",
"Estimate the uncertainty in the solution due to numerical error"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"fitz = @ode_def FitzhughNagumo begin\n",
" dV = c*(V - V^3/3 + R)\n",
" dR = -(1/c)*(V - a - b*R)\n",
"end a b c\n",
"p = [0.2,0.2,3.0]\n",
"u0 = [-1.0;1.0]\n",
"tspan = (0.0,20.0)\n",
"prob = ODEProblem(fitz,u0,tspan,p)\n",
"cb = AdaptiveProbIntsUncertainty(5)\n",
"monte_prob = MonteCarloProblem(prob);"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sim = solve(monte_prob,Tsit5(),num_monte=100,\n",
" callback=cb,abstol=1e-3,reltol=1e-1);\n",
"using Plots; gr(); plot(sim,vars=(0,1),linealpha=0.4)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sim = solve(monte_prob,Tsit5(),num_monte=100,\n",
" callback=cb,abstol=1e-6,reltol=1e-3)\n",
"plot(sim,vars=(0,1),linealpha=0.4)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Case Study: Parameter Estimation\n",
"\n",
"Instead of simply solving differential equation models, one can ask the inverse problem: given that I see data like `y`, what needs to be true about the model such that we can match the data?\n",
"\n",
"The act of estimating the parameters for a mechanistic differential equation model is called parameter estimation.\n",
"\n",
"- Parameter estimation required repeated solving of differential equations\n",
" - Solving equations must be fast\n",
"- It requires access to optimization routines (global optimization is necessary)\n",
"- Machine learning tools (regularization, cross validation, etc.) are a must\n",
"\n",
"Thus DifferentialEquations.jl is well-suited for handling these kinds of problems."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Problem Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"using DifferentialEquations, Plots; gr()\n",
"lotka = @ode_def_nohes LotkaVolterraTest begin\n",
" dx = a*x - b*x*y\n",
" dy = -c*y + d*x*y\n",
"end a b c d\n",
"\n",
"p = [1.5,1.0,3.0,1.0]\n",
"u0 = [1.0;1.0]; tspan = (0.0,10.0); \n",
"prob = ODEProblem(lotka,u0,tspan,p)\n",
"sol = solve(prob,Vern8(),abstol=1e-14,reltol=1e-14); \n",
"t = collect(linspace(0,10,200))\n",
"randomized = [(sol(t[i]) + .01randn(2)) for i in 1:length(t)]\n",
"using RecursiveArrayTools; \n",
"data = vecvec_to_mat(randomized); "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"
![](\"http://idahoptv.org/sciencetrek/topics/food_chain/images/web2.png\",width=800)
![](\"https://user-images.githubusercontent.com/1814174/27303063-525e9e3e-54ee-11e7-9ba4-4fd307923769.png\")
![](\"https://user-images.githubusercontent.com/1814174/27211626-79fe1b9a-520f-11e7-87f1-1cb33da91609.PNG\")
![](\"https://user-images.githubusercontent.com/1814174/27208182-0d0ddc70-51f9-11e7-8750-9c94dce654a3.png\")
![](\"https://user-images.githubusercontent.com/1814174/27303007-1867b652-54ee-11e7-9380-61cf6ab0ceed.png\")