Skip to content

Instantly share code, notes, and snippets.

@sdwfrost

sdwfrost/pomp_integrator.ipynb

Created November 29, 2017 13:19
Show Gist options
  • Save sdwfrost/94d4a61b3aef277fc3814fee8dea06a5 to your computer and use it in GitHub Desktop.
Save sdwfrost/94d4a61b3aef277fc3814fee8dea06a5 to your computer and use it in GitHub Desktop.
Download ZIP
Discrete SIR model using the DiffEq integrator interface
Display the source blob
Display the rendered blob
Raw
pomp_integrator.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"deletable": true,
"editable": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: Method definition non_autodiff_setup(Any, Any) in module OrdinaryDiffEq at /Applications/JuliaPro-0.6.0.1.app/Contents/Resources/pkgs-0.6.0.1/v0.6/OrdinaryDiffEq/src/misc_utils.jl:49 overwritten in module DiffEqCallbacks at /Applications/JuliaPro-0.6.0.1.app/Contents/Resources/pkgs-0.6.0.1/v0.6/DiffEqCallbacks/src/manifold.jl:32.\n"
]
}
],
"source": [
"using DifferentialEquations\n",
"using Distributions"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false,
"deletable": true,
"editable": true
},
"outputs": [
{
"data": {
"text/plain": [
"sir (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function sir(t,u,du,parms)\n",
" (S,I,R,Y)=u\n",
" (β,γ,α,ι,N,δt)=parms\n",
" λ=β*((convert(Float64,I)+ι)^α)/N\n",
" ifrac=1.0-exp(-λ*δt)\n",
" rfrac=1.0-exp(-γ*δt)\n",
" infection=rand(Binomial(S,ifrac))\n",
" recovery=rand(Binomial(I,rfrac))\n",
" du[1]=S-infection\n",
" du[2]=I+infection-recovery\n",
" du[3]=R+recovery\n",
" du[4]=Y+infection\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false,
"deletable": true,
"editable": true
},
"outputs": [],
"source": [
"γ=7.0/13.0\n",
"α=1.0\n",
"ι=0.5\n",
"N=403.0\n",
"M=10\n",
"R₀=5.0\n",
"δt=1.0/convert(Float64,M)\n",
"β=R₀*(1.0-exp(-γ*δt))/δt;"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false,
"deletable": true,
"editable": true
},
"outputs": [
{
"data": {
"text/plain": [
"(2.62110617498984, 0.5384615384615384, 1.0, 0.5, 403.0, 0.1)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u0=[60,1,342,0]\n",
"tspan=(0,540)\n",
"parms=(β,γ,α,ι,N,δt)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false,
"deletable": true,
"editable": true
},
"outputs": [
{
"data": {
"text/plain": [
"DiffEqBase.DiscreteProblem with uType Array{Int64,1} and tType Int64. In-place: true\n",
"timespan: (0, 540)\n",
"u0: [60, 1, 342, 0]"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = (t,u,du)->sir(t,u,du,parms)\n",
"prob=DiscreteProblem(f,u0,tspan)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[60, 1, 342, 0]\n",
"[60, 1, 342, 0]\n",
"[60, 1, 342, 0]\n",
"[60, 1, 342, 0]\n",
"[60, 1, 342, 0]\n",
"[60, 1, 342, 0]\n",
"[60, 1, 342, 0]\n",
"[60, 1, 342, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[60, 0, 343, 0]\n",
"[59, 1, 343, 1]\n",
"[59, 1, 343, 1]\n",
"[59, 1, 343, 1]\n",
"[59, 0, 344, 1]\n",
"[59, 0, 344, 1]\n",
"[59, 0, 344, 1]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 1, 344, 2]\n",
"[58, 0, 345, 2]\n",
"[58, 0, 345, 2]\n",
"[58, 0, 345, 2]\n",
"[58, 0, 345, 2]\n",
"[58, 0, 345, 2]\n",
"[58, 0, 345, 2]\n",
"[57, 1, 345, 3]\n",
"[57, 1, 345, 3]\n",
"[57, 1, 345, 3]\n",
"[57, 1, 345, 3]\n",
"[57, 1, 345, 3]\n",
"[57, 1, 345, 3]\n",
"[57, 1, 345, 3]\n",
"[57, 1, 345, 3]\n",
"[57, 1, 345, 3]\n",
"[57, 1, 345, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[57, 0, 346, 3]\n",
"[56, 1, 346, 4]\n",
"[56, 1, 346, 4]\n",
"[56, 1, 346, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[56, 0, 347, 4]\n",
"[55, 1, 347, 5]\n",
"[55, 1, 347, 5]\n",
"[55, 1, 347, 5]\n",
"[55, 1, 347, 5]\n",
"[54, 2, 347, 6]\n",
"[54, 2, 347, 6]\n",
"[54, 2, 347, 6]\n",
"[54, 2, 347, 6]\n",
"[54, 2, 347, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[54, 1, 348, 6]\n",
"[53, 1, 349, 7]\n",
"[53, 1, 349, 7]\n",
"[53, 1, 349, 7]\n",
"[53, 1, 349, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[53, 0, 350, 7]\n",
"[52, 1, 350, 8]\n",
"[52, 1, 350, 8]\n",
"[52, 1, 350, 8]\n",
"[52, 1, 350, 8]\n",
"[52, 1, 350, 8]\n",
"[52, 1, 350, 8]\n",
"[52, 1, 350, 8]\n",
"[52, 1, 350, 8]\n",
"[52, 1, 350, 8]\n",
"[52, 1, 350, 8]\n",
"[52, 1, 350, 8]\n",
"[50, 3, 350, 10]\n",
"[50, 2, 351, 10]\n",
"[50, 2, 351, 10]\n",
"[50, 2, 351, 10]\n",
"[50, 2, 351, 10]\n",
"[50, 2, 351, 10]\n",
"[50, 2, 351, 10]\n",
"[50, 2, 351, 10]\n",
"[50, 2, 351, 10]\n",
"[50, 2, 351, 10]\n",
"[50, 2, 351, 10]\n",
"[49, 3, 351, 11]\n",
"[49, 3, 351, 11]\n",
"[49, 3, 351, 11]\n",
"[49, 3, 351, 11]\n",
"[49, 3, 351, 11]\n",
"[49, 3, 351, 11]\n",
"[49, 3, 351, 11]\n",
"[49, 3, 351, 11]\n",
"[49, 3, 351, 11]\n",
"[49, 3, 351, 11]\n",
"[49, 3, 351, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 2, 352, 11]\n",
"[49, 1, 353, 11]\n",
"[49, 1, 353, 11]\n",
"[49, 1, 353, 11]\n",
"[49, 1, 353, 11]\n",
"[49, 1, 353, 11]\n",
"[49, 1, 353, 11]\n",
"[49, 1, 353, 11]\n",
"[49, 1, 353, 11]\n",
"[49, 1, 353, 11]\n",
"[49, 1, 353, 11]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[48, 2, 353, 12]\n",
"[47, 3, 353, 13]\n",
"[47, 3, 353, 13]\n",
"[47, 2, 354, 13]\n",
"[47, 2, 354, 13]\n",
"[47, 2, 354, 13]\n",
"[47, 2, 354, 13]\n",
"[47, 1, 355, 13]\n",
"[47, 1, 355, 13]\n",
"[47, 1, 355, 13]\n",
"[47, 1, 355, 13]\n",
"[47, 1, 355, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n",
"[47, 0, 356, 13]\n"
]
}
],
"source": [
"u=copy(u0)\n",
"du=u\n",
"for i in 1:tspan[2]\n",
" f(i,u,du)\n",
" u=du\n",
" println(u)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false,
"deletable": true,
"editable": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.205316 seconds (631.00 k allocations: 17.975 MiB, 3.42% gc time)\n"
]
}
],
"source": [
"@time begin\n",
"for j in 1:1000\n",
" u=copy(u0)\n",
" du=u\n",
" for i in 1:tspan[2]\n",
" f(i,u,du)\n",
" u=du\n",
" end\n",
"end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false,
"deletable": true,
"editable": true
},
"outputs": [
{
"ename": "LoadError",
"evalue": "\u001b[91mInexactError()\u001b[39m",
"output_type": "error",
"traceback": [
"\u001b[91mInexactError()\u001b[39m",
"",
"Stacktrace:",
" [1] \u001b[1mconvert\u001b[22m\u001b[22m at \u001b[1m./rational.jl:85\u001b[22m\u001b[22m [inlined]",
" [2] \u001b[1mInt64\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Rational{Int64}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./sysimg.jl:24\u001b[22m\u001b[22m",
" [3] \u001b[1m(::DiffEqBase.#kw##init)\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.DiscreteProblem{Array{Int64,1},Int64,true,##1#2,Void}, ::OrdinaryDiffEq.Discrete{true,false}, ::Array{Array{Int64,1},1}, ::Array{Int64,1}, ::Array{Any,1}, ::Type{Val{true}}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./<missing>:0\u001b[22m\u001b[22m (repeats 2 times)",
" [4] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:515\u001b[22m\u001b[22m"
]
}
],
"source": [
"integrator=init(prob,FunctionMap(scale_by_time=false);dt=δt)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.6.0",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
@sdwfrost
Copy link
Author

@ChrisRackauckas any idea what's going wrong here?

@ChrisRackauckas
Copy link 

using OrdinaryDiffEq
using Distributions

function sir(t,u,du,parms)
    (S,I,R,Y)=u
    (β,γ,α,ι,N,δt)=parms
    λ=β*((convert(Float64,I)+ι)^α)/N
    ifrac=1.0-exp(-λ*δt)
    rfrac=1.0-exp(-γ*δt)
    infection=rand(Binomial(S,ifrac))
    recovery=rand(Binomial(I,rfrac))
    du[1]=S-infection
    du[2]=I+infection-recovery
    du[3]=R+recovery
    du[4]=Y+infection
end

γ=7.0/13.0
α=1.0
ι=0.5
N=403.0
M=10
R₀=5.0
δt=1.0/convert(Float64,M)
β=R₀*(1.0-exp(-γ*δt))/δt;

u0=[60,1,342,0]
tspan=(0,540.0)
parms=(β,γ,α,ι,N,δt)

f = (t,u,du)->sir(t,u,du,parms)
prob=DiscreteProblem(f,u0,tspan)
@time solve(prob,FunctionMap(scale_by_time=false);dt=δt)

The issue was tspan=(0,540) meant time was integers when you should've been using Float64 (or Rationals to be exact). I think this is the kind of "quick" problem where you'll measure our keyword argument + startup time though,

0.001946 seconds (10.97 k allocations: 1.105 MiB)

is definitely too much and mostly due to kwarg + splat penalty for the big type. Hopefully that's fixed in 1.0 though!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment