notebooks/kr08/2_6/julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program 2.6: SEIR model in Julia using DifferentialEquations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Author: Lloyd Chapman @LloydChapman\n",
    "\n",
    "Date: 2018-10-01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DifferentialEquations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "seir_ode (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function seir_ode(dY,Y,p,t)\n",
    " dY[1] = p[4]-p[1]*Y[1]*Y[3]-p[4]*Y[1]\n",
    " dY[2] = p[1]*Y[1]*Y[3]-(p[2]+p[4])*Y[2]\n",
    " dY[3] = p[2]*Y[2] - (p[3]+p[4])*Y[3]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 365.0)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "par=[520/365,1/60,1/30,774835/(65640000*365)]\n",
    "init=[0.8,0.1,0.1]\n",
    "tspan=(0.0,365.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DiffEqBase.ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true\n",
       "timespan: (0.0, 365.0)\n",
       "u0: [0.8, 0.1, 0.1]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seir_prob = ODEProblem(seir_ode,init,tspan,par)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol=solve(seir_prob);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "R=ones(1,size(sol,2))-sum(sol,1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip2600\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip2600)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip2601\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip2600)\" points=\"\n",
       "62.7871,360.121 580.315,360.121 580.315,11.811 62.7871,11.811 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip2602\">\n",
       "    <rect x=\"62\" y=\"11\" width=\"519\" height=\"349\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  77.4341,360.121 77.4341,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  211.197,360.121 211.197,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  344.96,360.121 344.96,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  478.722,360.121 478.722,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  62.7871,350.263 580.315,350.263 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  62.7871,267.386 580.315,267.386 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  62.7871,184.51 580.315,184.51 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  62.7871,101.633 580.315,101.633 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  62.7871,18.7565 580.315,18.7565 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  62.7871,360.121 580.315,360.121 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  62.7871,360.121 62.7871,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  77.4341,360.121 77.4341,354.896 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  211.197,360.121 211.197,354.896 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  344.96,360.121 344.96,354.896 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  478.722,360.121 478.722,354.896 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  62.7871,350.263 70.55,350.263 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  62.7871,267.386 70.55,267.386 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  62.7871,184.51 70.55,184.51 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  62.7871,101.633 70.55,101.633 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  62.7871,18.7565 70.55,18.7565 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 77.4341, 373.921)\" x=\"77.4341\" y=\"373.921\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 211.197, 373.921)\" x=\"211.197\" y=\"373.921\">100</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 344.96, 373.921)\" x=\"344.96\" y=\"373.921\">200</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 478.722, 373.921)\" x=\"478.722\" y=\"373.921\">300</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 56.7871, 354.763)\" x=\"56.7871\" y=\"354.763\">0.00</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 56.7871, 271.886)\" x=\"56.7871\" y=\"271.886\">0.25</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 56.7871, 189.01)\" x=\"56.7871\" y=\"189.01\">0.50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 56.7871, 106.133)\" x=\"56.7871\" y=\"106.133\">0.75</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 56.7871, 23.2565)\" x=\"56.7871\" y=\"23.2565\">1.00</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 321.551, 397.6)\" x=\"321.551\" y=\"397.6\">Time</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(-90, 14.4, 185.966)\" x=\"14.4\" y=\"185.966\">Proportion</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  77.4341,85.0577 77.5801,89.1438 78.476,112.789 80.1123,150.581 82.2278,191.482 84.9658,234.328 88.2871,274.152 92.1507,306.607 96.3066,328.221 100.528,340.177 \n",
       "  104.698,345.944 108.781,348.486 112.777,349.542 116.721,349.966 120.685,350.133 124.785,350.197 129.172,350.221 134.027,350.229 139.573,350.231 146.128,350.232 \n",
       "  154.178,350.231 164.534,350.229 178.225,350.227 193.423,350.222 209.436,350.216 226.888,350.208 245.94,350.198 267.026,350.183 290.631,350.162 317.551,350.131 \n",
       "  348.985,350.084 386.791,350.008 432.222,349.888 480.97,349.718 532.385,349.496 565.668,349.329 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#e26f46; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  77.4341,317.112 77.5801,313.09 78.476,289.991 80.1123,253.807 82.2278,215.964 84.9658,178.37 88.2871,146.365 92.1507,124.31 96.3066,114.691 100.528,115.142 \n",
       "  104.698,121.435 108.781,130.304 112.777,139.951 116.721,149.613 120.685,159.107 124.785,168.551 129.172,178.18 134.027,188.253 139.573,199.045 146.128,210.876 \n",
       "  154.178,224.144 164.534,239.361 178.225,256.687 193.423,272.752 209.436,286.685 226.888,299.014 245.94,309.736 267.026,318.98 290.631,326.819 317.551,333.354 \n",
       "  348.985,338.673 386.791,342.853 432.222,345.876 480.97,347.712 532.385,348.783 565.668,349.205 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#3da44d; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  77.4341,317.112 77.5801,317.169 78.476,317.359 80.1123,317.096 82.2278,315.815 84.9658,313.002 88.2871,308.449 92.1507,302.363 96.3066,295.674 100.528,289.347 \n",
       "  104.698,283.904 108.781,279.474 112.777,275.997 116.721,273.342 120.685,271.376 124.785,269.997 129.172,269.165 134.027,268.906 139.573,269.322 146.128,270.601 \n",
       "  154.178,273.06 164.534,277.214 178.225,283.68 193.423,291.325 209.436,299.244 226.888,307.246 245.94,314.967 267.026,322.191 290.631,328.728 317.551,334.463 \n",
       "  348.985,339.321 386.791,343.251 432.222,346.148 480.97,347.918 532.385,348.944 565.668,349.342 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2602)\" style=\"stroke:#c271d2; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  77.4341,350.263 77.5801,350.142 78.476,349.406 80.1123,348.062 82.2278,346.284 84.9658,343.844 88.2871,340.579 92.1507,336.265 96.3066,330.958 100.528,324.878 \n",
       "  104.698,318.261 108.781,311.28 112.777,304.055 116.721,296.624 120.685,288.929 124.785,280.8 129.172,271.98 134.027,262.157 139.573,250.946 146.128,237.836 \n",
       "  154.178,222.11 164.534,202.74 178.225,178.952 193.423,155.245 209.436,133.4 226.888,113.076 245.94,94.6447 267.026,78.1909 290.631,63.8362 317.551,51.5971 \n",
       "  348.985,41.4671 386.791,33.4325 432.222,27.6333 480.97,24.1967 532.385,22.3229 565.668,21.6688 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip2600)\" points=\"\n",
       "489.608,108.291 562.315,108.291 562.315,32.691 489.608,32.691 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  489.608,108.291 562.315,108.291 562.315,32.691 489.608,32.691 489.608,108.291 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  495.608,47.811 531.608,47.811 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 537.608, 52.311)\" x=\"537.608\" y=\"52.311\">y1</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#e26f46; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  495.608,62.931 531.608,62.931 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 537.608, 67.431)\" x=\"537.608\" y=\"67.431\">y2</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#3da44d; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  495.608,78.051 531.608,78.051 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 537.608, 82.551)\" x=\"537.608\" y=\"82.551\">y3</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip2600)\" style=\"stroke:#c271d2; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  495.608,93.171 531.608,93.171 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip2600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:start;\" transform=\"rotate(0, 537.608, 97.671)\" x=\"537.608\" y=\"97.671\">y4</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(sol.t,[sol',R'],xlabel=\"Time\",ylabel=\"Proportion\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.3",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}