Skip to content

Instantly share code, notes, and snippets.

@sdwfrost

sdwfrost/ModelingWithJulia.ipynb

Last active August 17, 2017 11:48
Show Gist options
  • Save sdwfrost/8454788 to your computer and use it in GitHub Desktop.
Save sdwfrost/8454788 to your computer and use it in GitHub Desktop.
Download ZIP
Some examples of mathematical modeling and simple MCMC in Julia
Display the source blob
Display the rendered blob
Raw
ModelingWithJulia.ipynb
{
"metadata": {
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Ordinary differential equations"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"using ODE"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function SIR(t,x,p)\n",
" (S,I,R)=x\n",
" (beta,gamma)=p\n",
" dS=-beta*S*I\n",
" dI=beta*S*I-gamma*I\n",
" dR=gamma*I\n",
" return([dS,dI,dR])\n",
"end\n",
"\n",
"# Initialise model\n",
"t = 0.0:0.01:200\n",
"inits=[0.99,0.01,0.0];\n",
"p=[0.1,0.05];\n",
"\n",
"# Run model\n",
"result=ode4((t,x)->SIR(t,x,p),t,inits);"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Dataframes"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"using DataFrames;"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"df=DataFrame();\n",
"df[\"t\"]=result[1];\n",
"df[\"I\"]=result[2][:,2];"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"head(df)"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 2,
"text": [
"6x2 DataFrame:\n",
" t I\n",
"[1,] 0.0 0.01\n",
"[2,] 0.01 0.0100049\n",
"[3,] 0.02 0.0100098\n",
"[4,] 0.03 0.0100147\n",
"[5,] 0.04 0.0100196\n",
"[6,] 0.05 0.0100245\n"
]
}
],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Plotting"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"using Winston;\n",
"wp = plot(df[:,\"t\"],df[:,\"I\"])\n",
"setattr(wp,\"xlabel\",\"t\")\n",
"setattr(wp,\"ylabel\",\"I\")\n",
"wp"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEsCAIAAADfNCTgAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3de1yM6f8/8Cs11Zw6OC2jzWBtbBRhrVisj1Sy48yUfJNDiJRo1bRLIqEoVLaEZBm+265D6WRXtHa31War3SwqTTlE55oO0xya3x/z/fbzRUzmcM3c837+8XmYu6trXtcO78/c933d16UnlUoRAACA99UHdwAAANBuUEYBAEAhUEYBAEAhUEYBAEAhUEYBAEAhUEYBAEAhUEYBAEAhUEYBAEAhUEYBAEAhUEYBAEAhUEYBAEAhUEYBAEAhUEYBAEAhBrgDqERTU9Mvv/yCOwUAQNNZWVl9/PHHCnZCzG+jjo6Ojx496uzsFIlE8rSPj4+/evWq/P1XVFRcu3ZN/vYSieT48ePyt0cIxcTEqLR9UFAQj8eTv/3t27eLiorkb19VVZWamtqrSKoesre3d6/aX7hwoa6uTv72v/322927d+Vv//Tp08uXL/cqkqqH3Nv+v//++5qaGvnb//HHH3/++af87aurq3/88cdeRZJzyBKJpLOzs7Oz08nJqVf9vxExv41KJJK7d+9WVFQMHTp069at72yflZX10Ucfffnll3L2X1RUpK+vL397sVicnZ0tf3uEUGpqqkrbR0REzJgxw8bGRs72HR0dffv2nT17tpztS0pKxGKxRg35q6++6lX7v/76y8HBwcLCQs72QqGQRqM5OjrK2f7BgwdtbW0aNeTe9v/333//5z//YTKZcrbv6uoyMDBwcXGRs31ZWVlTU5MqhpyRkXHz5k2EkFgslr/znhCzjBoZGQUFBY0aNQp3EACAJnJ2dnZ2dkYIZWVlKaE7KRFxOJyamhr52587dy4nJ0dlcaRSqbSrq0uj2n/zzTfV1dW9+hVVU/WQ161b16v2GkjVQ+5t/xqot0NWyt8KPNdG+Xw+m82mUqkMBiM2NvaNbWJiYuzs7Egk0tq1a1/5UXJyspWVlZGR0ciRI3///XfV51UCPT09jWqvgXRwyL2l6iHr4H9SpcBzUu/v719bW8vj8UpLS52dna2trWfOnPlKGwaDERIS8sMPP7xyPC0tLTAw8OTJk5MnT378+LGpqamaQgMAwJvoSdW+pZ1IJDI3N8/IyPj8888RQhs3buzo6EhKSnpjYz8/v9bW1sTExO4jdnZ2vr6+Hh4eb3mL4OBgPz+/AQMGyBlJ9h9Bp/6vuKurq08fYs7T6AkMWRf0dsheXl4JCQkKvimG/8Q8Hq+trc3W1lb20tbWtqSkRM7f7ezsLCwsfPHihYWFxeDBgzdv3tzR0aF4JD09PZ2qoQghXfvXhWDIugHLkDGc1Le2tiKE6HS67KWpqSmfz5fzd58+fSqVSlNTUwsKCoRCIYvF2rt3b1hY2CvNqqqqpk6dqq+vT6fTx40bhxBiMpkcDkd5gwAAaJ8rV67IZnxXVlZWVVUhhMhksuLdYiijNBoNIcTn801MTBBCzc3N3SX1nWRj3r59+wcffIAQ8vPzi46Ofr2MWlpaHj58WP6TegCALpg/f/78+fNfPuLl5aV4txjKKJPJpFAoxcXF06ZNQwgVFRVZW1vL+buDBw/u37+/rp2A67iysrKCgoLffvtNIpE8ePCAyWSSSKQxY8ZMnTrVxsYG/jIA7DCUURKJ5OrqGhoaevHixbKyMi6X2/1IXGJiIoPBmDt3LkJILBaLxWKJRCKRSAQCgYGBgYGBAULI09MzMjLS3t5eJBIdPXp03rx56h8CUIOWlpbExMQ//vijo6Nj7ty5CxYsGDJkSP/+/RsaGqqrqx89erR///6KiooZM2Zs2LBh2LBhuPMCHab41NP30NLSsmzZMgqFMmjQoJiYmO7jjo6OO3bskP05ODj45Zzbtm2THRcIBOvWrTMxMRk4cKCPj09HR8fr/fd2+j3QKAKBIDw8nM1mR0dHNzY2vqVlR0fH+fPnnZycVq1aVVFRoa6AgDiUMv0enmICmiU/P3/x4sXh4eFv/D/InmRkZHz66af79+8XCoWqywaIR4ufYgLgdVKp9MCBA1u3bj1w4EBgYKCxsbH8v+vk5HT79m2RSDRjxoyysjLVhQTgdRim36vBpk2b2trajI2Nhw0btmPHDtxxwLsJBAIPD4+xY8cGBgbKLoK/n4KCgo0bN3711VdLlixRYjxAMKmpqbKZTzweLzMzU8HeiFlGg4ODt2zZMmDAAB2cV6+N+Hz+8uXLv/zyy40bNyreW2Njo7u7+6xZs/z9/eHTB28kOxlHCG3YsEErn2JSjz59+vTp0wf+FWm+5ubmuXPnenl5KaWGIoTMzc2vXr1aUVHh5eUlkUiU0icgGD09PVmJUEpvhC2jQCvw+XwWi+Xr67tgwQIldquvrx8TEzNkyJCAgAClrMsLwFtAGQXYSCQSd3f39evXq+g6ZkhIyODBg93c3KCSApWCMgqw8ff3nzhxopubm+reIiAgwNbWlsPhdHV1qe5dgI6DMgrwOHv2rEgk+vrrr1X9RsHBwXp6etu2bVP1GwGdRcy9mNrb2y9cuECn0wcMGCD//llAbQoKCmJjY2/cuKGee4D79+/38PCIi4vr7U6ZgKgKCwsLCwsRQvX19Yr3RswyamBgMGLECDMzM/nXjgJq09LSsm7dupSUFAqFop531NPTO3ny5NKlSz/88MNebTMJiKpfv36y7el79ZRHT4hZRg0NDSdNmgQL5Wkmb2/v0NDQ4cOHq/NNSSTS6dOnZ82aNWLEiE8++USdbw000Icffvjhhx8ihHrad6NX4NooUCsul0sikbCsy2Vubs7lcvft29fU1KT+dwcEBmUUqM/Tp08PHDhw5MgRXAFGjRq1ZMmSlStXwo17oERQRoH6+Pv7R0dHy3Y9wGXBggXW1tb79+/HmAEQDJRRoCbnzp0zMzN7fSdt9du7d++DBw9u376NOwggCGLeYhKLxWVlZQ0NDWQy2dLSEnccgBoaGiIiInJycnAHQQghAwODsLCwefPm3bhxo2/fvrjjAAzq6+vr6uoQQgKBQPHeiFlGBQJBdnY2jUaTPQuIOw5AkZGRISEh5ubmuIP8DwsLi127dm3evPn8+fO4swAMHjx48PvvvyOEmpubFe+NmGWURqN5e3vDhCcNcefOnYKCgn379uEO8n8sXLjw9u3bZ86c8fDwwJ0FqJu9vb29vT1C6MGDB4r3BtdGgWpJpVLZnSXcQd4gNDT02LFjlZWVuIMA7YanjPL5fDabTaVSGQxGbGzsG9vExMTY2dmRSKS1a9e+/tOqqioajTZ79mwVJwWK4nK5n3766ejRo3EHeQMqlXrkyJH169cTcvFyoDZ4Tur9/f1ra2t5PF5paamzs7O1tfXrN3AZDEZISMgPP/zwxh62bNliZ2en8qBAMW1tbQcPHtSQO0tvNHXq1ClTphw7dmzLli24swBtheHbqEgk4nK5ISEhAwYMsLe3d3Nze+PzWIsWLWKxWG+8KZGamioSiVgslsqzAsVERESsXbtWc+4svdGOHTvOnj1bXl6OOwjQVhjKKI/Ha2trs7W1lb20tbUtKSmR/9fb29sDAgKOHj2qmnRAaZ4/f56amrp+/XrcQd7B2Ng4JiYGTu3Be8NwUt/a2ooQ6l57ydTUlM/ny//roaGhy5cvHzFixFvaVFVVTZ06VV9fn06njxs3DiHEZDI5HI4CqUGvffvttxwOh0Qi4Q7ybpMnT54xY0ZCQoLmF32giCtXrsg2BK2srKyqqkIIkclkxbvFUEZpNBpCiM/nyx4KbG5uln85u3v37l2+fFm2UOBbWFpaHj58GCY8YVReXv7TTz/t2rULdxB5+fv7T5s27csvv2QwGLizAFWZP3/+/PnzXz7i5eWleLcYyiiTyaRQKMXFxdOmTUMIFRUVWVtby/m7ubm5jx8/ZjKZCKH29vbOzk4LC4snT56oLi14PyEhISEhIVq0MyuVSt2/f/+WLVtSUlJwZwFaBsO1URKJ5OrqGhoa2tjYmJ+fz+VyV61aJftRYmJienq67M9isVggEEgkEolEIhAIZLuSeXh4lJeXy1au3rp16+TJkwsKCtQ/BPB29+7dq66u1rrpaI6OjsOHD09NTcUdBGgZPPNGo6KizM3NLSwsWCxWWFhY92ynlJSU3Nxc2Z9DQkLIZHJMTExSUhKZTA4MDEQIkcnkQf+LTqcbGhp+8MEHWIYA3uLUqVO7d+/GneJ9+Pv7f/PNN+3t7biDAG2iR8i7k8HBwX5+fnBtFIvi4uKAgICsrCzcQd5TXFzc48ePw8PDcQcB6uDl5ZWQkKBgJ/AwKFCyPXv2aNGdpdd5eXn9+eefDx8+xB0EaA1iLk3S1NS0bds2Y2PjYcOGBQUF4Y6jQ/7555/m5mbZog9aysDAYN++fT4+Ptr7hRq809WrV9PS0hBCSllRgZhl1MzMTPaUFO4gOic2NlYNW8+r2qRJk5hM5n//938vW7YMdxagEiwWS/YYpFImPMFJPVCa0tLS+/fvT58+HXcQJdi3b9++ffva2tpwBwFaAMooUJoDBw7s2LEDdwrl6Nevn7+//6FDh3AHAVoAyihQjqdPn5aUlDg5OeEOojQrVqy4du1aRUUF7iBA00EZBcpx+PDhrVu34k6hTPr6+hEREdu2bcMdBGg6KKNACRobG8vLyxcvXow7iJJNnz6dTCbfunULdxCg0Yh5p57P5x86dIhKpQ4ZMmT16tW44xBffHy8k5OTvr4+7iDKt2/fvqVLl/72228GBsT8x6Kbbt26JXtgsrq6WvHeiPk3g0KhsNnsfv36GRkZ4c5CfJ2dnVwuNy8vD3cQlRg6dKijo2NCQoK3tzfuLEBpxo8fP3z4cISQUpbrJmYZ1dfXHzJkCMwbVY9z584tXrxYKes2aqbAwEDZNg1mZma4swDlMDExkS3UaWhoqHhvcG0UKEQqlV65cmXjxo24g6gQlUoNDAyEDRdAT6CMAoVkZWUNHjyY8F/8ly9fnpmZWVZWhjsI0ERQRoFCoqOjfX19cadQuT59+oSHh8tWawTgFVBGwfsrKSnR19fXzD3olW7GjBlisfj27du4gwCNQ8xbTJ2dnbm5uaampqamppMmTcIdh7COHDmiC19Fu4WHh/v4+Fy/fl2LNkcBb/To0aNHjx4hhHq1n2ZPiPlttKurq7GxsaGhoaWlBXcWwqqvr+/o6HBwcMAdRH1Gjx79ySefXLhwAXcQoKiOjo6GhoaGhgbZ7kQKIua3UTKZPH/+fMLf98DrxIkTU6dO1bXvZTt37nRwcFi0aBFMSdZq1tbWsp00f/rpJ8V7I+a3UaBqYrGYy+WuXLkSdxB169+///Lly+Pj43EHARoETxnl8/lsNptKpTIYjNjY2De2iYmJsbOzI5FIa9eu7T4oFAq9vb1HjBhBoVBsbGwuXbqkrsjg/7hy5cqcOXOoVCruIBj4+fmdOnWqoaEBdxCgKfCUUX9//9raWh6Pl5KSwuFwbt68+XobBoMREhLi5ub28kGhUGhgYHDx4sUnT574+vqy2ez79++rKTR4SVxcHLGn3L+FsbHx1q1b9+3bhzsI0BQYyqhIJOJyubJNPmTP2CUlJb3ebNGiRSwWy9zc/OWDNBrt6NGjEydO7Nu375o1a4YOHQr71Kvf33//PXbsWNkjybrJ3d39l19+4fF4uIMAjYChjPJ4vLa2NltbW9lLW1vbkpKS9+jnxYsXlZWVNjY2r/+oq6urrq6upqamvr6+ra2tra1NIBAoFBq8JC4ubu7cubhT4KSvrx8aGnrs2DHcQUDviEQiWUFoamqqqampqanR1jv1ra2tCCE6nS57aWpq+h5Tt4RC4YoVKzw9PceOHfv6T589e7Z8+XIDAwMTE5Px48cjhCwtLQm2qDAuTU1Nd+7ciYuLwx0EM0dHx8OHD//111+yv2BAK/z000/Z2dkIIR6PJ9sTVCnrH2IoozQaDSHE5/NlK6w0Nzd3l1Q5icViNpttbm7e0+0pCwuLyMhImPCkCmfPnl25cqWuzXN6o/Dw8KCgoMzMTNxBgLycnZ2dnZ1fPqKtO4MymUwKhVJcXCx7WVRUJJvAJSeJROLm5iYSic6fP0/IdYI1mVQqPXPmzKpVq3AH0Qh2dnb9+/e/fv067iAAMwxllEQiubq6hoaGNjY25ufnc7nc7n+WiYmJ6enpsj+LxWKBQCCRSCQSiUAgkF3CkEgkK1eurKmp+e6772THJRKJ+oegs3JychwdHWHZzW579+7lcrldXV24gwCc8Ex4ioqKMjc3t7CwYLFYYWFhM2fOlB1PSUmRreyPEAoJCSGTyTExMUlJSWQyWba4zuPHj7lc7q1bt8zMzMhkMplMjoqKwjIE3RQXF7dkyRLcKTQIk8k0NTWFx0N1nJ5UKsWdQfmCg4P9/Pzg2qhyyW7c/fLLL7iDaJa6uro5c+bk5eUpZR11oGZeXl4JCQkKdgIPgwJ5nThx4uUnyoBM//79lyxZ8u233+IOArAh5rfRTZs28fl8Y2PjYcOGBQUF4Y5DBBKJxM7OLi8vj8B7Lr23trY2d3f3M2fOyCafAM139erVtLQ0hFBlZWVWVpaCvRFzhSczMzPZU1K4gxDHtWvXli1bBjX0jahUqrOzc2RkZGhoKO4sQC4sFovFYiHtnfAEtNHx48fh5tJbrF69Oj09/fnz57iDAAygjIJ3q6ioEIvFVlZWuINoLgMDAw6Hs2fPHtxBAAZQRsG7JSQkrFu3DncKTbdw4cKamhrYPVQHQRkF7yAUCtPT0xcuXIg7iKbT09Pz8fEJDg7GHQSoG5RR8A5paWkrVqwgkUi4g2iB6dOnt7e35+fn4w4C1AomPIF3mDNnzrfffqvLq4v2SnFx8TfffHPlyhXcQcDbwISnd4MJT8pSWlqqp6cHNVR+NjY2/fv3z8zMdHJywp0F9AgmPAH1SUxMVMrfM52yc+fO3bt3w3olugPKKOiRUCjMzMyU/Z82kN/QoUOnTJkC65XoDiijoEfp6emurq5wc+k9cDicpKSkzs5O3EGAOkAZBT2KjY1dtmwZ7hRaqX///jNnzoT1SnQElFHwZmVlZfr6+nBz6b35+vqePn26ubkZdxCgcsS8U9/e3s7lcul0+oABA+bNm4c7jlaCZfEURKVSN2/eHBERsXfvXtxZwKv++uuvwsJChFB9fb3ivRHz26iBgcHIkSNHjRrFZDJxZ9FKcHNJKTw9Pe/fv19dXY07CHjVgAEDRo0aNWrUKGNjY8V7I2YZNTQ0nDhx4pQpU8aMGYM7i1bKzMxcsGABLOeuIH19/RUrVoSEhOAOAl5lYWExZcqUKVOmUKlUxXsjZhkFCoqLi1u5ciXuFESwcOHCBw8e/Pvvv7iDABXSrDLK5/PZbDaVSmUwGD3tQR8TE2NnZ0cikeDKnYo8evRILBZ/9NFHuIMQRFhYGIfDwZ0CqJBmlVF/f//a2loej5eSksLhcG7evPl6GwaDERIS4ubmpvZ0uuLkyZNr1qzBnYI4pk6dOnjw4F9//RV3EKAqGlRGRSIRl8uVPQtvb2/v5uaWlJT0erNFixaxWCxzc3O1B9QJYrE4NTV10aJFuIMQiq+vb1BQECGXAQJIoyY88Xi8trY2W1tb2UtbW9uTJ0++X1cCgeD69eumpqZGRkZDhgxBCJHJZLhrL49r1645ODgYGRnhDkIoVlZWY8aM+fHHHxcvXow7i06rra2tq6tDCDU2NjY2NiKElDKxV4PKaGtrK0KITqfLXpqamvL5/PfrqqWl5ezZs0ZGRnQ63c7ODiH0wQcfQBmVR0JCwoEDB3CnIKCdO3e6uLiwWCx4uBaj0tLSP/74AyHE4/EqKysRQg0NDYp3q0FllEajIYT4fL5sl9rm5ubuktpbAwcO3LdvHyyU11tVVVUdHR0wS0wVBg0atGrVqpMnT27YsAF3Ft1lb29vb2//8hGiLZTHZDIpFEpxcbHsZVFRkbW1Nd5IuubkyZOrV6/GnYKwVq9effz48ZaWFtxBgJJpUBklkUiurq6hoaGNjY35+flcLnfVqlWyHyUmJqanp8v+LBaLBQKBRCKRSCQCgUAsFmNLTCwSieTSpUuwi7LqyB4P3b9/P+4gQMk0qIwihKKioszNzS0sLFgsVlhY2MyZM2XHU1JScnNzZX8OCQkhk8kxMTFJSUlkMjkwMBBbXGJJT0+fNWuWUp6NAz1ZvXr19evXHz9+jDsIUCZi7sUUHBzs5+cH10Z7hc1m79y585NPPsEdhOCysrLOnz9/5swZ3EEAQgh5eXklJCQo2IlmfRsFuDx58uTZs2dQQ9XA0dGxurr67t27uIMApdGgO/VK1NTUtG3bNtgZVH5wc0mdDh48GBAQcP36ddxBdNfLO4Mq3hsxyyjsDNorEonkxx9//P3333EH0RXjxo0bPXr0tWvXXFxccGfRUbAzKFCy7Ozs6dOnUygU3EF0yI4dO3bt2iUSiXAHAUoAZRSgq1evwi7KajZkyJB58+bBZk3EAGVU1z19+rSkpGTs2LG4g+icgICAxMTEpqYm3EGAoqCM6rpTp07BzSUsqFQqh8OJjo7GHQQoCsqoTpNIJCkpKbCLMi5Lly7NysoqLS3FHQQohJh36mHCk5wyMzNnzJgBN5dw6dOnz4EDB7766qtLly7hzqJbYMLTu8GEJzllZWWtX78edwqdNn369JiYmJycnC+++AJ3Fh0CE56Acjx58gSW0dIEBw8e5HA4EokEdxDwnqCM6q4TJ07AnkuagMlkzpo168SJE7iDgPcEZVRHicXiy5cvw80lDREUFHT8+HHZthZA60AZ1VFpaWmzZ8+GZfE0BI1G27FjR1xcHO4g4H1AGdVRmZmZcHNJo7i6umZmZv7777+4g4BeI+adej6fHxERQaVShwwZsnbtWtxxNE55efmjR48+/vhj3EHA/6enp3f48GE/P7+srCzcWYjv5s2bt27dQghVV1cr3psSyugbd5OX6d4FRM2oVKqrq2u/fv0MDQ2xBNBwCQkJ8BC9Bpo0aZKFhcWlS5cWLlyIOwvBTZgwYeTIkQghHo+neG9KKKNveZoNVxnt06fP4MGDYd7oGwkEgszMzL179+IOAt4gPDx82bJlTk5OZDIZdxYio9Ppso2HlbLftRLKaGFhoeKdALVJSUmZP38+7JaumQYOHLho0aIDBw6EhITgzgLkhecWE5/PZ7PZVCqVwWDExsb2qk1paamDg4OJicngwYP9/f1h0nJv/fzzz3BGr8m8vb3T09MrKipwBwHywlNG/f39a2treTxeSkoKh8O5efOm/G08PDyGDh1aXV19+/btK1euKL4dlU4pKipqbm62sLDAHQT0yMDAIDIyMioqCncQIC8MZVQkEnG5XNkz7/b29m5ubq/fpHpLm5KSkv/6r/+iUqkjRoyYM2dOSUmJmvNrtdjY2E2bNuFOAd5h+vTpjY2NqampuIMAuWCY8MTj8dra2mxtbWUvbW1tT548KX8bFxeX5OTkCRMmPH/+/Pr160eOHHn9LVpaWo4cOUKlUqlU6qhRoxBCpqamkydPVtWQtERTU1NBQUF8fDzuIODdIiIinJycZs+eDfealKisrOzRo0cIoerq6mfPniGEZP+rIAzfRltbWxFCsttkCCFTU1M+ny9/m8jIyF9//ZVGo3300UcODg5v3BRMIpE0NjbW19c3NDS0tLS0tLS0tbWpaDha5MyZMx4eHnp6eriDgHcbNGjQ6tWrDx48iDsIoXR2dsoKgqw+1NfXK2U7LAzfRmk0GkKIz+ebmJgghJqbm7vL5TvbiESi//znP66urgEBAfX19a6uriEhIa/f0zQ3N/fz84MJTy+TSqV5eXmw+Y8W2bRp09y5cx88eGBlZYU7C0FYW1u/sqSZti6Ux2QyKRRKcXGx7OUb12rrqQ2Px7t//76vry+ZTLawsFixYkVmZqY6w2uvrKwsU1NTU1NT3EGAvPT19ffs2ePj44M7CHgHDGWURCK5urqGhoY2Njbm5+dzudzuWfqJiYnp6elvaWNpaWlmZhYbGysUCl+8eMHlcm1sbNQ/BG0UGxsL/yC1zqeffvrRRx+dO3cOdxDwNngmPEVFRZmbm1tYWLBYrLCwsJkzZ8qOp6Sk5ObmvqWNkZHR5cuXr1692q9fvzFjxlhYWOzfvx/LELRLeXl5R0cHrNCsjfbt25eWltbQ0IA7COgRnqVJ6HT6xYsXXz/+8hl6T21mzJiRl5enwnBEFBcXt3nzZtwpwPswMzObP39+YGAgTJHWWMRcKE8qlba2tvL5/Pb2dtxZ8GttbX327NmXX36JOwh4T2w2+8mTJ90nakBxQqGQz+fz+fyuri7FeyPyQnlkMtnS0tLX1xd3HMzOnj1rZ2enr6+POwh4f7GxsQEBAZMnTzYyMsKdhQhycnKys7ORkuaN6kmlUsV70TTBwcEw4UlGKpVOnDjx+vXrffv2xZ0FKOTQoUMtLS27d+/GHYRQvLy8FL9aQsyTetAtOzt70qRJUEMJwNfX9/r16/D0swaCMkpwR44cgcsaxGBgYBAbG7t9+3alXM4DSgRllMju378/fPjw0aNH4w4ClGP8+PETJ048duwY7iDg/4AySmTR0dHz5s3DnQIoE4fDOXPmDKxGqlGgjBJWfX19fn6+o6Mj7iBAmchk8tGjR9evX0/Im8NaipgTnpqamrZt22ZsbDxs2LCgoCDccfCIj49fv349rOdEPNOmTZswYcKJEydgF4P3dvXq1bS0NIRQZWWl4r3BhCdiEgqFU6ZMuX37NqxWSUhtbW2ff/755cuXLS0tcWfRbjDhCfTo/PnzCxYsgBpKVFQqNSoqCk7tNQSUUQKSSqWxsbHr16/HHQSo0IwZM+zt7eFBe00AZZSAMjMzx48fP3DgQNxBgGpt3749Pj6+vLwcdxBdB2WUgKKiovz9/XGnACpHJpPj4uK8vLxgm3G8oIwSzd27dz/++GPZRn6A8D777LO5c+ceOnQIdxCdBhOeiObAgQN+fn64UwD12bJlyzSDnJoAABh+SURBVPTp0x0dHbt30gXvpNwJT8Qso2ZmZrI97nEHUbfy8vKampopU6bgDgLUh0QinTp1ysPDIzc319jYGHcc7cBisVgsFtLeLe2A6kRGRm7fvh13CqBuo0eP9vHx2bFjB+4gOgpPGeXz+Ww2m0qlMhiM2NjY3rZJTk62srIyMjIaOXLk77//rpbIWqC6urqwsHDu3Lm4gwAMVq5c+ejRo4yMDNxBdBGek3p/f//a2loej1daWurs7Gxtbd29q90726SlpQUGBp48eXLy5MmPHz+GHYO7HT582MfHB57+1FknT56cM2fOhAkTYK6buknVTigUUqnU3Nxc2csNGzZ4eHjI32b8+PFJSUlvfwsOh1NTU6PEzJqvrq5u3LhxYrEYdxCA088//zxv3ryuri7cQbTGunXrFO8Ew0k9j8dra2vrvqtoa2v7+oLePbXp7OwsLCx88eKFhYXF4MGDN2/e3NHRoc7wGis6OtrHxwc2XNJxs2bNsrKyioqKwh1Et2A4qW9tbUUI0el02UtTU1M+ny9nm6dPn0ql0tTU1IKCAqFQyGKx9u7dGxYW9sqv19fX+/j4GBsb0+n0cePGIYQGDBgguzFHSM3NzVeuXCkoKMAdBOAXHh7+xRdfTJ8+feLEibizaJyCgoLCwkKEUGVlZVVVFULo8ePHineLoYzSaDSEEJ/PNzExQQg1Nzd3l8t3tpGttbF9+/YPPvgAIeTn5xcdHf16GTU2Np4yZQqdTieTycOHD+/ukKhOnDixfft2EomEOwjAj0QinTt3buHChTdu3DAzM8MdR7N88MEHY8aMkf2ByWQihGSzRxWEoYwymUwKhVJcXDxt2jSEUFFRkbW1tZxtBg8e3L9//3feRaFSqW5ubjoyb7S1tTU5OTk/Px93EKAphg4dGhwc7OnpeenSJdxZNIuFhYWFhcXLR/Ly8hTvFsO1URKJ5OrqGhoa2tjYmJ+fz+VyV61aJftRYmJienr629t4enpGRkbW1dVVV1cfPXoUNsmIjY1du3YtbF8OXrZ48eJhw4YlJibiDqIbFL9L9R5aWlqWLVtGoVAGDRoUExPTfdzR0XHHjh1vbyMQCNatW2diYjJw4EAfH5+Ojo7X+9edO/V8Pn/MmDHt7e24gwCN09nZOW3atF9//RV3EI2mlDv1sPq9douOjqbRaGvXrsUdBGiiqqoqFouVnZ0NM0l7Aqvf67qWlpbTp0+7u7vjDgI0lKWlZWRkpL+/v1gsxp2FyIhZRsVi8YMHD+7du0fsfWiPHDmyfv16WI0CvMXs2bOtrKyCg4NxB9EstbW19+7du3fvnkAgULw3Yq7w1NnZefPmTRqNNmjQoGHDhuGOoxKNjY0//PDDnTt3cAcBmu7rr79euHDh999/v3TpUtxZNEV5ebnsHn1LS4vivRGzjFKp1PXr1xP72mhcXNyOHTsMDQ1xBwGaTk9P7+zZs3PmzLGysrKxscEdRyN89tlnn332GULo3r17ivdGzJN6wquurr506RJ8uQByotPpSUlJq1evbmhowJ2FgKCMaqU9e/YEBQUZGBDzZAKogpWV1Z49ewIDA0UiEe4sRANlVPuUlpb+/fffixYtwh0EaBlnZ+eRI0du2bIFdxCigTKqfQ4fPnzw4EFYVxS8h4CAAIFAALvbKxeUUS2Tl5f39OlT2G0JvLdvv/02LS0N1slXImJeXCPwzqCBgYE9bbsCgDyMjIxOnz7t4OBgYWExduxY3HHwgJ1B342oO4P++OOPo0aNen1BLAB6pV+/fhcvXvzmm28OHTo0ZMgQ3HEwgJ1BdZRQKIyIiNi9ezfuIIAIRo4c6ePjs2jRoubmZtxZtB6UUa1x9OhRFoslW68aAMVNnTr1q6++2rRpk1AoxJ1Fu0EZ1Q41NTXnzp3bunUr7iCAUBYvXjxlyhQPD4+uri7cWbQYlFHtcPTo0Z07d8IqJEDpNm3aNHz48G3btuEOosWgjGqBP//8Mz8/f+HChbiDAGLau3cviUSCy+7vjZjLNm/atInP5xNjwlNXV9eMGTPi4+M/+eQT3FkAYUkkEjab/fnnn+vIM04vT3jKyspSsDdillEirX6fmJj48OHDgwcP4g4CCK6zs3Pr1q2ffvpp975nukApq98Tc94oYdTV1Z0/f/7q1au4gwDiMzIyioyMdHFxoVAoy5Ytwx1Hm+C5Nsrn89lsNpVKZTAYPT2T8/Y2VVVVNBpt9uzZqg+L0/bt2zdt2kSj0XAHATqBQqHIznavXLmCO4s2wfNt1N/fv7a2lsfjlZaWOjs7W1tbz5w5s1dttmzZYmdnp87M6peTk9PQ0LB48WLcQYAOodPpsn3L9fT0ZM/5gHfC8G1UJBJxuVzZw5r29vZubm5JSUm9apOamioSiYj9GQsEghMnTsTExOAOAnSOmZlZamrqmTNnrl27hjuLdsBQRnk8Xltbm62treylra1tSUmJ/G3a29sDAgKOHj2qtsBYhIaGTpo0ydLSEncQoIvMzc0TExP37t0LZ/fywHBS39raihCi0+myl6ampnw+X/42oaGhy5cvHzFixFveorq6esGCBSQSycTERHbuz2AwlLIGgXoUFhb+8ssvt27dwh0E6C5zc/OMjIw1a9YIhULCbFdz48aN3NxchBCPx+PxeAghpTwIi6GMym6Y8Pl8ExMThFBzc3N3uXxnm3v37l2+fLmwsPDtb9G/f393d3dzc3MSidS3b1+EkBZt/SYSiaKjo+Pj4/v0gYcjAE5mZmanT59esGBBe3u7h4cH7jhKMGnSpFGjRiGEWltb29raEEKHDh1SvFsMZZTJZFIolOLi4mnTpiGEioqKXl/5rac2ubm5jx8/ZjKZCKH29vbOzk4LC4snT5688uskEmns2LFaOm80PDz8448/hsn2QBOYmJikpqauWbOmra3N29sbdxxF0en0V760USgUxbvF8H2HRCK5urqGhoY2Njbm5+dzudzu6b6JiYnp6elvaePh4VFeXl5YWFhYWLh169bJkycXFBSofwiqU1RUdP369a+++gp3EAD+B5VKTU5OzsnJ2bVrF+4sGgrPaWNUVJS5ubmFhQWLxQoLC+ueyZSSkiK7ctFTGzKZPOh/0el0Q0NDIi0c19nZGRYWlpiYCFt+Ao1iaGh48eLFpqamgIAAiUSCO47GgYdBNYifn99HH320efNm3EEAeAOpVBoSEvL333+fP3+eMIuNKeVhULiJoSmys7NLS0s3bdqEOwgAb6anp7d79+65c+cuWbKkvr4edxwNQsyTx87Ozhs3bpiampqZmX322We447zb8+fP4+PjT506BdsmAw23du1aBoPh4OBw4cKFjz/+GHec91RWVlZWVoYQamlpUbw3Yn4blUqlHR0d7e3tAoEAd5Z36+rqWr169fr164l0nRcQ2Ny5c0+dOrVy5UrtndosEona29vb29uVsuw/Mb+NGhsbu7i4aMu10f3799va2s6ZMwd3EADkNW7cuEuXLi1ZssTb29vd3R13nF4bPXr06NGjEUKZmZmK90bMb6NaJDc3t7S0dM+ePbiDANA7DAbjp59+Sk1N9fX1FYvFuOPgBGUUp2fPnvn4+ISFhcEMJ6CNKBTKhQsXBgwYsGHDhrq6OtxxsIEyio1QKHR3dz9y5AiDwcCdBYD3pKen9/XXXy9cuHD27NkEexZGflBGsdm6deuyZcteX2gVAK3j4uLy/fff7969W/E5mNqImGVUKpW2tLQ0NzfLVh/QQPHx8Xw+f8OGDbiDAKAcI0eOvHDhwq+//hoYGChboU2TdXZ2Njc3Nzc3K+WhLGJekuPz+YcPH6ZQKJaWlj4+PrjjvConJyc7O/u7777DHQQAZaJQKGfOnDl9+vTMmTOTkpLGjBmDO1GPbt26df36dYTQ8+fPFe8NHgZVt4cPHy5fvjwjI2PQoEG4swCgEv/888++ffs+//zzDRs2aPgTJfAwqPZ58eKFq6vr2bNnoYYCAhszZsypU6dKSkq8vLxqa2txx1E5KKPq09raunDhwqioKE0+2QFAKYyNjWNiYpYsWTJr1izC70QCZVRNhELh4sWLfXx8pk+fjjsLAGri6Oh448aNmzdvuru7NzY24o6jKlBG1UEikXh6ei5btszV1RV3FgDUasCAAVFRUfPmzfP09Lx06RLuOCoBZVTlurq61q1bZ2VltWbNGtxZAMCDzWYnJCRcuHBh5cqV1dXVuOMoGTHv1G/atInP5xsbGw8bNiwoKAhjEqlUumPHDjKZvHv3bowxANAQaWlpHA7H29t73bp1+vr6uGJcvXo1LS0NIVRZWZmVlaVgb8Qso5oz4cnb29vQ0DA6Ohp3EAA0RWtra2Rk5E8//XTs2LHx48fjDQMTnjRaV1dXYGAglUqFGgrAy2g0WkhISFxcXEhIyObNmwlw6wlPGeXz+Ww2m0qlMhiM2NhY+dsIhUJvb+8RI0ZQKBQbGxuNvWItFotXrlzZp0+fiIgI3FkA0EQ2NjaXL1/+7LPPHBwcjh8/rtU75eEpo/7+/rW1tTweLyUlhcPh3Lx5U842QqHQwMDg4sWLT5488fX1ZbPZ9+/fV3P4d+ro6Ni4ceP48eP37duHOwsAmktPT8/d3T0nJ+fp06eTJ0/OyMjAneh9SdVOKBRSqdTc3FzZyw0bNnh4eLxHG6lUOnLkyO++++714xwOp6amRomZ5VdXVzd16tRTp05heXcAtBSPx3N3d3dwcLh7964633fdunWKd4Lh2yiPx2tra7O1tZW9tLW1LSkpeY82L168qKystLGxUXVg+ZWWlq5Zs4bD4Xh6euLOAoA2GTp06NmzZ8PDw7/55pulS5c+fPgQd6JewLDCk2wRLTqdLntpamrK5/N720YoFK5YscLT03Ps2LGvv0VVVdXUqVP19fXpdPq4ceMQQkwmk8PhKHso/8fNmze3bNmSnJwse0cAQG9NmDAhLS0tJydn8+bNQ4YM2blz57Bhw5TY/5UrV65du4YQqqysrKqqQgiRyWTFu8VQRmk0GkKIz+ebmJgghJqbm7vLpZxtxGIxm802Nzfv6faUpaXl4cOH1TnhKS4uLiMjIzMzE5ayB0BBX3zxxRdffJGRkeHr6ztgwAAOhzNixAil9Dx//vz58+e/fMTLy0vxbjGc1DOZTAqFUlxcLHtZVFRkbW0tfxuJROLm5iYSic6fP49x+m43gUAQFBSUl5f3/fffQw0FQFmcnZ2vXr26dOnSr7/+esWKFd3VQANhKKMkEsnV1TU0NLSxsTE/P5/L5a5atUr2o8TExPT09Le0kUgkK1eurKmp+e677yQSiUAgwDtPory8/Isvvhg6dGhycrKxsTHGJAAQkpOTE5fL3bBhw9GjR11cXG7duoU70ZsofpfqPbS0tCxbtoxCoQwaNCgmJqb7uKOj444dO97SpqKi4pX8ERERr/evnjv1Fy5cYLPZf/zxh6rfCAAglUoLCwu3bt366aefnj9/XiAQKKVPpdypx7OJCJ1Ov3jx4uvHMzMz396GyWRKNeDp1aamptDQ0Orq6vj4eNnVWwCAqtna2h4+fPjZs2fJycl2dnYrVqzw8PAYMmQI7lzwMGjvZWVlTZ8+fdKkSVwuF2ooAGrGYDACAwMLCgqGDx++dOnStWvXZmdnd3V1YYxEzKVJ/Pz8LC0t6XT6wIEDX7kxp4iGhoZt27ZRKJTg4GC4mwSAJigoKEhISHj8+LG9vb2Hh8eHH34o52/dvXsXIZSRkfHjjz8qmIGYO4MaGhqOGTPG3NxcNnFKKc6ePZuSkuLq6spms5XVJwBAQRMmTIiPj29tbb148aK/v/+qVatcXFze+VuDBw+Wze9Wyj0rYpZREok0fvx4Zc0bLSoq2rlzJ4PBSEpKMjc3V0qfAAAlotFoa9askX9ldAaDITuhpFAoir87Mcuosjx//nz37t337t2Lioqys7PDHQcAoIngFtObtbS0hISEODk5zZs37+bNm1BDAQA9gTL6qvb29gMHDkybNm3o0KH5+fkuLi56enq4QwEANBeU0f+vvr5+7969U6ZMMTExycvL8/T0JJFIuEMBADQdMa+NisXi+/fvv3jxgkKhDB8+/J3tW1paOBxOXl7exo0b8/PzDQ0N1RASAIBLTU1NTU0NQqijo0Px3ohZRjs7O3/55RcajTZo0CB5yiiFQnFycjp69GifPvD1HADiq6io+OOPPxBCr6/S+R6IWUapVOq6devkn/BkYGAwb948lUYCAGiOyZMnT548GSH0zz//KN4bfPkCAACFQBkFAACFQBlFCKHm5ub29nbcKdSqtrZWLBbjTqFW1dXVuCOoGwxZPaCMIoTQtWvX7ty5gzuFWh07dqyurg53CrXatWsX7gjqBkNWDyijAACgEGKW0erqan9/fy8vr/DwcHnad3R0dHZ2yt9/Z2dnY2NjryI9f/5co9q3t7f3av8VPp/fq+seQqFQ0/4T9fa6TV1dXa+ue/D5/La2NvnbC4XChoaGXkVS9ZB72399fb1IJJK/fWtrq2zTXzmJRKLe/ieSc8hXr1718vLy8vJSyk7OxCyjDx48CA4OTkhICAoKkqf9X3/99e+//8rf//3798+fPy9/e7FYHBoaKn97hNDOnTtV2v7OnTv19fXyt8/IyPjtt9/kb19aWpqcnNyrSKoeckFBQa/ax8XF9aqsZGdn3759W/72FRUVp0+f7lUkVQ+5t/3Hx8c/ffpU/vY///xzrxamq6ysTExM7FUkOYfMYrESEhISEhJevHjRq/7fiJhlFAAA1IaY0+9NTEySk5P79u0rZ/uHDx/26s51bW3tkydP5H+MrKurq7y8PDIyUs72CKGKigqVtm9qajp37lx2drac7R88eGBsbFxYWChn+7q6usrKyl6d8al6yJ2dnb1qX1hYeOLECTqdLmf7hw8fkkikv//+W872DQ0Nqv5b0dsh97b/u3fvdnR0mJqaytm+rKysT58+8p/5NTU1PXjwQKVDNjIykr9xT4i5iYhAIKiqqsKdAgCg6fr169evXz8FOyFmGQUAALWBa6MAAKAQKKMAAKAQKKMAAKAQKKMAAKAQKKMAAKAQKKMAAKAQKKMAAKAQXS+jfD6fzWZTqVQGgxEbG4s7jqrMmzdP73/RaLTu4wQbfkxMjJ2dHYlEWrt27cvHexqmtg+/p/ES9eMWCoXe3t4jRoygUCg2NjaXLl3q/hHej5iYD4PKz9/fv7a2lsfjlZaWOjs7W1tbz5w5E3colTh+/PiqVasQQnp6et0HCTZ8BoMREhLyww8/vHK8p2Fq+/B7Gi8i6MctFAoNDAwuXrw4fPjwS5cusdnsoqKiUaNGIewfsVSHCYVCKpWam5sre7lhwwYPDw+siVTFxcXlxIkTrxwk6vB9fX3XrFnT/bKnYRJm+K+MV6ozH/fIkSO/++47qQZ8xDp9Us/j8dra2mxtbWUvbW1tS0pK8EZSnd27d3/44YezZs36+eefZUd0ZPg9DZPYwyf8x/3ixYvKykobGxukAR+xTp/Uy1aQ7V7Cx9TUVCmbVmugjRs3Dho0iEqlpqamuri45OXljRs3TkeG39MwCTx8wn/cQqFwxYoVnp6eY8eORRrwEet0GZVdfefz+SYmJgih5uZm+VdF0y4uLi6yP4waNer27ds//PDDuHHjdGT4PQ2TwMMn9sctFovZbLa5uXn3LSPsH7FOn9QzmUwKhVJcXCx7WVRUZG1tjTeSGhgaGsq2D9GR4fc0TB0ZPsE+bolE4ubmJhKJzp8/r6+vLzuI/yNWxQVXLbJmzRoHB4eGhoY7d+6Ymprm5OTgTqR8ra2t586de/r0aV1d3alTpwwNDX/99VfZjwg2fJFI1NHRsXnz5lWrVnV0dIhEItnxnoap7cN/43gJ/HGLxWJXV9cZM2Y0NTV1dHR0dHSIxWLZj/B+xLpeRltaWpYtW0ahUAYNGhQTE4M7jkrw+fzp06ebmppSKJTx48enpKR0/4hgww8ODn75K8K2bdtkx3saprYP/43jJfDHXVFR8cq3wIiICNmP8H7EsGwzAAAoRKevjQIAgOKgjAIAgEKgjAIAgEKgjAIAgEKgjAIAgEKgjAIAgEKgjAKdYGFh8eeff+JOAYgJyigAACgEyiggPm9v7+fPn7NYLCaTmZycjDsOIBp4ignoBAsLi8uXL0+cOBF3EEBA8G0UAAAUAmUUAAAUAmUU6IQ+feCvOlAV+LsFdMLAgQPLyspwpwDEBGUU6AQOhxMQEGBubn7ixAncWQDRwJ16AABQCHwbBQAAhUAZBQAAhUAZBQAAhUAZBQAAhUAZBQAAhUAZBQAAhUAZBQAAhUAZBQAAhUAZBQAAhUAZBQAAhUAZBQAAhUAZBQAAhUAZBQAAhUAZBQAAhfw/TIxW5YI5TBAAAAAASUVORK5CYII=",
"prompt_number": 3,
"text": [
"FramedPlot(...)"
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Stochastic simulations"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"using DataFrames\n",
"using DataArrays\n",
"using Distributions;"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"prompt_number": 3
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function sir(beta,gamma,N,S0,I0,R0,tf)\n",
" (t,S,I,R) = (0,S0,I0,R0)\n",
" ta=DataArray(Float64,0)\n",
" Sa=DataArray(Float64,0)\n",
" Ia=DataArray(Float64,0)\n",
" Ra=DataArray(Float64,0)\n",
" while t < tf\n",
" push!(ta,t)\n",
" push!(Sa,S)\n",
" push!(Ia,I)\n",
" push!(Ra,R)\n",
" pf1 = beta*S*I\n",
" pf2 = gamma*I\n",
" pf = pf1+pf2\n",
" dt = rand(Distributions.Exponential(1/pf))\n",
" t = t+dt\n",
" if t>tf\n",
" break\n",
" end\n",
" ru = rand()\n",
" if ru<(pf1/pf)\n",
" S=S-1\n",
" I=I+1\n",
" else\n",
" I=I-1\n",
" R=R+1\n",
" end\n",
" end\n",
" results = DataFrame(t=ta,S=Sa,I=Ia,R=Ra)\n",
" return(results)\n",
"end"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 4,
"text": [
"sir (generic function with 1 method)"
]
}
],
"prompt_number": 4
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Simulations"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"N = 100\n",
"@elapsed [sir(0.1/10000,0.05,10000,9999,1,0,1000) for i in 1:N]"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 5,
"text": [
"0.37622139"
]
}
],
"prompt_number": 5
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"srand(1)\n",
"s = sir(0.1/10000,0.05,10000,9999,1,0,1000)\n",
"wp = plot(s[:,\"t\"],s[:,\"I\"])\n",
"setattr(wp,\"xlabel\",\"t\")\n",
"setattr(wp,\"ylabel\",\"I\")\n",
"wp"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEsCAIAAADfNCTgAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dd0BT9/4//ncgkMHeYWkYigiCFMSNCBQHiqU4UcSCg/Zav7S91raiV6la9bbWq22tFq2tWwStOGvcuIpVhoAUCEOGyA6B5GT+/jj3w/UnQokZJzl5Pf6Sd/I+eR6EF2e8z/tNkcvlCAAAwJsyIDoAAADoNiijAACgFCijAACgFCijAACgFCijAACgFCijAACgFCijAACgFCijAACgFCijAACgFCijAACgFCijAACgFCijAACgFCijAACgFCijAACgFCijAACgFALK6I4dO/z9/alUampq6isv7d2718PDg06n+/j4lJWV4Y0dHR2xsbFMJpPNZh8/frznzX21AwCAJlE1/5Gurq6bN28+cODAK+0ZGRlbtmz5+eefAwICysvLLSws8PaUlBQMwxoaGh48eBATExMQEODl5dVPOwAAaBKFqNnvlyxZ4uLismnTpp4WX1/f9evXz5079+W3YRhmZWV169atoKAghNCsWbNGjhy5cePGvto1vBcAAKAt10a7urqKioqePHnCYrFcXV3Xrl0rk8kQQlwuVyAQ+Pn54W/z8/MrKirqpx0AADSMgJP616qrq0MI3bt3r6ioqLW1dcqUKYMGDVqxYgWfz6fRaMbGxvjbzM3N+Xw+Qqiv9t5Wr179559/IoQcHBzMzMwQQl5eXlZWVn0l6ejo6LmeMEBCoZBOp5Oji2Z2H8MwGo2m1i5yuVwkEinapbOz09zcXN2fIhaLe3501dRFJpNJJBKFukgkEgzDTExM1PopMplMKpUaGRmptYtUKpXL5VTqq/WtsbGxsrISIdTe3t7c3IwQCg8PX7t27cC3/FraUkYZDAZCaPXq1TY2NjY2NkuXLr106dKKFStMTU0xDBOJRPh/FY/HMzU1RQj11d7b48ePAwMDEUIsFgsvo25ubv1Uivnz5yt6w2rr1q2fffYZObpoYPelUumOHTtWr1498C4ikei77777+OOPB96lu7s7PT191apVA+9SX1+/c+fO/fv3D7xLW1tbRkbG8uXLB96lqanp7NmzSUlJA+/y/Pnzy5cvJyQkDLxLbW3tzZs3Fy5cOPAud+7cOXv27LZt2wbepaKi4vHjx7Nnzx54l9LS0pKSknfeeWfgXYqKiqqrq2fOnDnwLvn5+Q0NDdOmTXulncVi2djYIITa2tpaWloQQjk5OQPfbF+0pYy6uLhYWlr2bnd3d6fT6QUFBfg10MLCQn9//37ae/P09Pz3v/898CQMBiM4OFih8Pb29qTpooHdl0gkDg4OCnURCoUsFkuhLp2dnefPn1eoS2VlpampqUJdXrx4cefOHYW61NXVPXr0SKEu1dXVRUVFCnWxtrbmcrkKdWlqarKyslKoC4PBaG1tVagLlUoVCASK/ozJ5XKFuojFYjqdPpAuycnJCiV5LQKujUokEqFQKJVKpVIp/g+EEIVCWbx48TfffNPe3l5ZWXngwIGoqCiEEI1GmzdvXlpaGo/H43A4HA4nLi6un3blKXp+ihCys7MjTRcN7D6FQrG1tVWoi4GBwRt0wY87Bs7Q0BA/K1Koi7W1tUJdqFRqP9eUXsvIyOgNurz2uKT/LoruvrGxsaKXgGg0mkKXTRBCdDr9Dbrgp54aIte4NWvWvBxg9+7deHt3d3dCQoKZmZmTk9O6devwqxtyubytrS0mJobBYAwaNOjYsWM92+mr/RUrVqxQKJ6i7ycZfd791tbWzz77jOgUhCkuLt65cyfRKTRNJT/wBJzUb926devWrb3bGQzGwYMHDx48+Eq7paVlVlZW7/f31Q4AAJqkLQOetIezszPREYikz7tvaGjIYrGITkEYOp2u6JUTgCNs+L3GREVFjR07FiH0zjvv+Pr6Eh0HAEC8+/fvczgchFBubu5vv/2m5Na05U69+jg4OCxZsgQhpOitAAAAWfn6+rq4uCCEuFyu8lsjfxk1NjbGv1+ABORyOYVCQQgJBAJFbysD0MPU1BQfaa7QswN9gWujQGeIxeKCggIej4cQ+uOPP4iOA8B/QRkFOqOgoODYsWNXrlxpaWnBH3IjOhEACOnDST0gjYqKivPnz9+7dy8vL+/w4cMdHR3btm1TdJA5ACoHR6NAB+DjSa5du1ZaWtrc3Hzy5MnGxkYmk1lTUyMSifA5JgAgCvmPRsvLy/HHZhMTExV9khdoiadPn3Z1dWVnZ9PpdD6fX1tba2Rk5O/v39LSUlZWlpubiw/GAGCAOBzOqVOnEELV1dXKb438ZdTT0/PHH38kOgVQyuXLl/H5ZDdv3vzNN98ghBYvXmxnZ5eamrpy5crTp0/HxsZq9BlqoOMiIiIiIiKQiqYmIX8ZBbpOKBR++umntra2R48eZbFYdXV1N27cWL16NZfL3bt375kzZ/Ly8l68eAFlFBAFyijQditXrhSLxRiGhYaGyuXy8PBwmUw2ZMgQNzc3Gxubmzdv8vn8CxcufPjhh0QnBXoKbjEBrcbn80tKSoYOHYrPQEyhUJydnT/55BOEEJVKTUpKcnNzMzMzy8jIaG1tlUgkbW1tREcGegeORoH2am1t5XK5tbW1X3/9dUBAAN44dOjQnsUhxo0bN3r06I6Ojvj4+PXr148dO9bS0hKfqRYAjYEyCrTXuXPnUlNTLS0tZ8yY0fPo58sL7Pj5+dFoNKlUamlpmZ+ff/v27ZiYGCijQMPIX0ZhwJPuWrt2bVtbW2xsbF+Pz/dMiv7+++/v2rVLLBbn5eUhhMrKyoYMGaK5oEDXwIAnxcCAJ93V3Nw8Y8aMgdw7mj59+pYtW4RC4d27d69evVpeXi4UCn19ffF5TAB4hWoHPBFwi2nHjh3+/v5UKjU1NbX3q8XFxcbGxi+Ppu7o6IiNjWUymWw2++VFK/tqB+QgkUgcHR1nz57t7u7+t2/29fVdvnz53Llz7ezs3n333f/85z9fffXVrVu3NJATAAKORl1dXTdv3nzgwIHXvvqPf/zjlVPvlJQUDMMaGhoePHgQExMTEBDg5eXVTzsgh/Lyck9Pz3Hjxg3kzQYGBqmpqX/++aednV1NTU1WVtbz588pFAqbzR48eLC6owI9R8DR6Jw5c2bMmPHapf6OHDni7Ow8YcKEnhYMw06cOLFhwwYLC4vIyMiIiIijR4/20w5Io6qq6q233nJ1dR3g+y0tLUNDQwMDA4OCgoyNjTs6Oh4/fpyVlVVVVaXOmABo07hRHo+3cePG7du3v9zI5XIFAoGfnx/+pZ+fH/5QYF/tgDS+/PJLRZduNjQ0dHR0DAwMdHd3l8lk9fX1t2/fhjIK1E2Lyuj69euXLVvm5OT0ciOfz6fRaD0zVJubm/P5/H7aAQnIZLKysrKmpqZZs2Yp2tfOzm7YsGH//ve/nZycOjs7c3Nz6+rq1BESgB7acqe+oKCAw+E8fvz4lXZTU1MMw0QiEV4xeTwePvV/X+295eTk+Pj4IIS8vLzs7e0RQsnJySNHjlTr7gBlNDU1HTp0iM1me3p6Ktp30KBBFhYWLS0tSUlJBQUF169fv379+rvvvgsrjgCE0K1bt/Crf3V1dfgqTCr5wdCWMnrr1q3Kykr8Qhifz5dKpaWlpffu3XN3d6fT6QUFBUFBQQihwsJCf39/hFBf7b1NmDABBjzpkOLi4i+//PLhw4fTp09/g+74LM6DBw+eNWvW0KFDLSwsbt26dfPmzalTp6o6KdA9ISEhISEhL7fo6oAniUQiFAqlUqlUKsX/gRBKTEysqKjIy8vLy8t77733oqOjs7OzEUI0Gm3evHlpaWk8Ho/D4XA4nLi4uH7aga7bsWPHmTNnuFzuggUL3ngj1tbW/v7+s2bNCgwM5PF4J0+eVGFCAF5BQBlNTU1lMBiHDx/eunUrg8HYs2cPQojJZLL+j4mJCYPBsLW1xd+/c+dOKpXKYrGSkpLS09N7RjX11Q50WktLy7p16yIjIwd+j/61qFSqmZlZVFSUjY1NTk6OquIB0BsBJ/Vbt27dunVr/294+UtLS8usrKzeb+urHei0Fy9ejB492s3NzdHRUfmteXh4UCiU2tpamUxmYKBFN1QBmcAPFtAulZWVFhYWISEhqqp6H374oVgsvn79ukqengagNyijQIvk5OS0tLSYm5s7Ozurapv29vb29vbbt2+HBzSAmmjLnXr1gRmedMjDhw+XLFnCYrFUuM3hw4fPmjXr0KFDKizNQNfBDE+KgRmedMiVK1dWrVr12geF39jQoUNdXFz4fH5OTg6Xyx3IRCeA9HR+hicAXksikfz5558ODg6q3SyFQnnvvfcQQlVVVb/88otqNw4AgjIKtMe1a9f4fL63t7fKt+zg4DBx4kQKhZKZmdnc3Kzy7QM9B2UUaItDhw6Fh4fTaDSVb9nAwGDz5s2urq5FRUVr165tbW1V+UcAfQZlFGgFDMOuX79uaGiopu2PHz9+4sSJVCp13759mZmZavoUoJ/IX0YxDKuurq6uru7q6iI6C+hTbW2ti4vL22+/rabtGxgYhIaG4kM17t27p6ZPAbqCx+PhZUEkEim/NfLfqX/x4gU+YDA6Ohqf6glooYsXL5qZmb08Y7fKzZkzRy6X379/H5+uAeizp0+fXr16FSGkkis85C+jrq6un3/+OdEpwN9oaWmZOnWqWv/OMZnM0NDQlStXfvfddxUVFR4eHur7LKDlgoOD8VMTGPAEyKOysnLkyJHqfuzdxsYmNTWVQqG8ssgCAMqAMgqIJ5fLHz16pNpR969lZmaGT49/8uRJmUym7o8DegLKKCBeQUFBdXW1ap8B7cfGjRs7Ojrq6+s183GA9KCMAuJdv36dQqEoOcHowEVFRYWFheXm5mrm4wDpkf8WU0dHx40bNxBC3t7eKn/QEKhEdXV1WFiYxj6OTqc7OTkdOnQoJiZGYx8KtMqzZ88qKioQQp2dncpvjfxHoyKR6Pnz58+fPxcKhURnAa9RVVV1584dDVc0c3Nzlfz+AB3V3d2NlwWxWKz81sh/NGpnZzd//nyiU4A+paWllZWVKbPy0hvw9/c/cuSIRCKhUsn/KwB68/Lywpcdwk9VlUTA0eiOHTv8/f2pVGpqampPo1AoXLFihZubm4mJyahRo17et46OjtjYWCaTyWazjx8//rftQIdkZGT89ttvdDpdw+Vs/Pjx7e3tLS0tmvxQQFYElFFXV9fNmzdHR0e/3IhhmImJSVZWVk1NzZw5c6Kjo3ueLkhJScEwrKGhYd++fUlJSaWlpf23Ax1y5MgRqVSq+VVdnZycGAzGixcvNPy5gJQIKKNz5syZMWPGK4MELSwsduzYERAQYGNj8+mnnyKESkpKEEIYhp04cWLDhg0WFhaRkZERERH4k519tQPd8vTp0+jo6EWLFmn4c83NzSdOnFheXq7hzwWkpI23mEpKSjAMGzZsGEKIy+UKBAI/Pz/8JT8/v6Kion7agQ7Jzs5ubGz88MMPNf9cpoGBgbGx8TfffKPhzwWkpHXX17u6uuLj49etW2djY4MQ4vP5NBrN2NgYf9Xc3JzP5/fT3ltJScnSpUsRQl5eXnZ2dgihsLCwQYMGaWBfQD9aWlrWr1+fkpISGBhIyNLHY8eO7X+hb0A+paWl+PxedXV1lZWViJRrMWEY9s477wQGBvbcfTI1NcUwTCQS4RWTx+OZmpr2096bpaXlmDFjEEIuLi6WlpYIITMzM83sDuhHXl5eQ0PD+PHjiVo+furUqWvXrsUwTB0TRQPtZGVlhZ/m2tra4qPI29ralN+sFpVRsVg8e/Zse3v7PXv29DS6u7vT6fSCgoKgoCCEUGFhob+/fz/tvTk6OuJHo0Cr1NXVjRo1avjw4UQFGDp0KEKIx+Ph5yhAH+Crbb/c8vDhQ+U3S8CBgEQiEQqFUqlUKpXi/0AISaXSBQsWSKXSvXv3ikQioVCIzxxBo9HmzZuXlpbG4/E4HA6Hw8Hv6vbVDnRFfX39kCFDnJyciApgYmJiaGhYVlZGVABAGgSU0dTUVAaDcfjw4a1btzIYDPzYs7KyMjMzE5+7l8FgMBiMrKws/P07d+6kUqksFispKSk9PR0fNNtPO9AJubm5kZGRBAagUChubm4cDofADIAcKHK5nOgM6pWcnAzr1GshPz+/GzduWFtbE5hh0aJFjx494nA4BB4UA2KppD5o44AnQHpyudzAwIDYGooQWrlyZXNz8969e0l/MAHUSotuMalJdXX1v/71L4RQbGxszzhTQKy//vpLA5M0/63Ro0e3t7dfu3YtLi7O3d3dyMiI6ERAQ+7evXv58mWEkEqmnSV/GXVycvrggw8QQhYWFkRnAf/17NkzbZi0kEKhuLu7NzQ0JCQk7N+/H1Y81B8jR47EH/pYu3at8lsj/0m9kZGRg4ODg4MDnU4nOgtACCG5XP7gwQNtKKMIobFjx7a2tv7555/wFJxeYTKZeFlQyZw45C+jQNuIRKKSkhJ82CbhPDw8urq67O3ti4uLic4CdBWUUaBpQqGwoKCgr8clNGz8+PE+Pj4wTQlQBpRRoGkYhj179iwkJIToIAghNGbMmFWrVq1Zs+bZs2fPnj0jOg7QSeS/xQS0DZ/Pt7W1pVAoRAdBCCEGg7Fw4UIjI6Pa2toHDx5obFk9QCbkL6Pl5eXJyckIocTExODgYKLjACQWi8eOHUt0iv/Bxzk1Nzenp6fPnj2b6DhAEzgczqlTpxApZ3hSB09PT3iKSavU1tba2toSneJVTCazvLycz+f3NVUYIJOIiIiIiAiEEH6MpSS4Ngo07a+//tLChy+DgoLq6uquXLlCdBCge6CMAk0rKCiYOHEi0Sle9fnnnwuFQnztcgAUAmUUaFR9fX1mZqYWrj7g5ORkZmaWmZlJdBCge6CMAo0qLi42NDRksVhEB3kVm82ePn16bm5ud3c30VmAjoEyCjTqwYMHgYGBWjLa6RWff/65VCrNz88nOgjQMeS/Uw8DnrRKY2Pj3LlziU7xej4+PiwWKzs7W6vGYwF1UO2AJ5i2GWhOW1ubh4dHbW0tk8kkOsvrbdu27ciRI48fPzY0NCQ6C9AEXZ22eceOHf7+/lQqtWf5T1xHR0dsbCyTyWSz2cePH3/jdqC1nj59GhUVpbU1FCE0c+bM0tLSb7/9luggQJcQUEZdXV03b94cHR39SntKSgqGYQ0NDfv27UtKSiotLX2zdqC1Pv74Y2dnZ6JT9Gfw4MEWFhZnz57FMIzoLEBnEFBG58yZM2PGjFcmP8cw7MSJExs2bLCwsIiMjIyIiDh69OgbtAOt1dDQcP/+fW9vb6KD9IfJZEZHR+fl5W3cuJHoLEBnaMudei6XKxAIehb58PPzw6fRVbQdaK3jx49TKJRZs2YRHaQ/FAolMjLSwMBg27ZtlZWVRMcBukFbyiifz6fRaMbGxviX5ubmfD7/DdqBdhKJRJmZmSEhIZaWlkRn+RseHh5WVlYymQz+MIMB0pYBT6amphiGiUQivDLyeDx8hghF23vLycnB19jx8vKyt7dHCCUnJ48cOVJTewYQQqiqqio/P7+kpIToIH9v5MiRISEhVVVVRUVFM2bMIDoOUKVbt27hV//q6uq4XC5CiMFgKL9ZbSmj7u7udDq9oKAgKCgIIVRYWIjPjq5oe28TJkyAAU+Ea2pq6u7u1sKHl3ozNDSMj48/dOhQQUEB0VmAioWEhLwyX7iuzvAkkUiEQqFUKpVKpfg/EEI0Gm3evHlpaWk8Ho/D4XA4nLi4uDdoB9rp3r17VlZWKlk+TAP8/f0nTJjw8OFDmPAJDIhc49asWfNygN27d+PtbW1tMTExDAZj0KBBx44d63m/ou2vWLFihfr2BQyETCYzMzOLiIggOogC8DOY9957j+ggQL1UUh/gKSagXnK5vLGx0d/f/9dff50yZQrRcQaqpqZm+PDh1tbWhw8f1pJlo4A66OpTTEB/SCSSH3/88eeffx49erSWLAU6QCwWa968eQ0NDbt374Y5n0D/yF9G5XK5RCKRSCSkP+7WQg0NDUeOHPnqq68WLlyoE/eXehgbG48bN47FYhUXF8O9JvJRbVnQjUv+yuByuSkpKQihhISEUaNGER1Hv9TW1lZVVfn4+EyePJnoLAobO3asl5fXtWvXLl26FBwcbGBA/mMO/XHt2rXTp08jhFSyqjb5y6iHh8d3331HdAo9lZubW19fHxYWho/Y1S3Dhg0bNmxYXl7eo0eP8vPzAwICiE4EVCY8PDw8PBzp7oAnoD8ePnw4duzYadOmER3kTRgYGEyePHnSpEl3794tLCwkOg7QXlBGgRoVFBR88sknvWfz0hUTJ04MDQ21tLQ8ffq0Ss7+AClBGQXq0tTU1NjYOHXqVBMTE6KzvCF7e/t58+YFBgYWFRX98ssvRMcBWgrKKFCX0tJSLy8vbZ6keSDs7e3HjRtXVlZWVlZGdBagpch/i0koFJaXlyOEWCxWX9OXAHU4e/asDo2378e4ceMsLCzy8/O5XK67uzvRcYAKtLe3Nzc3I4RUMj83+ctoS0vLmTNnEEJRUVFaPmcwyZw6derSpUtEp1ABf39/NptdWVkJZZQ0ysvLb9y4gRBqb29XfmvkL6POzs7//Oc/iU6hd+RyeWtrq6enJ9FBVMDY2NjBwcHc3Pynn36aOHEijUYjOhFQVlBQED45HH6qqiS4NgrUoq2tzd7enjRD1keOHBkeHn769OlHjx4RnQVoHZL8lANtU1tbS6ZnxhYsWGBnZycWix8+fEh0FqB1yH9SDzRPJpNt2rSJzWYTHURlRo4cyWAwqFRqXl4e0VmA1oGjUaB6NTU1Fy9eHDp0KNFBVGno0KG+vr5//PHHX3/9RXQWoF3IX0bb29uvXLly5cqV58+fE51FXzQ0NIjFYg8PD6KDqBKFQlm4cOGTJ0/OnTtHdBagrOrqarws8Hg85bdG/jIqkUh4PB6PxxOLxURn0RcnT550cXEhWRlFCM2ePRshdOrUKXzlG6C7RCIRXhZU8l9J/mujtra2sbGxRKfQL9nZ2Y6OjoMGDSI6iIrRaDRPT8+ysrJVq1Z9//33RMcBb27IkCFDhgxBCKlkuS0tOhotKSmZPHmymZmZi4vLpk2beto7OjpiY2OZTCabzT5+/PjftgNiSaXSqqoq/GeUZJhM5ieffNLa2pqenk50FqBFtKiMLl68eOTIkS0tLVeuXNm1a9f58+fx9pSUFAzDGhoa9u3bl5SUVFpa2n87IFZNTQ2TyVy4cCHRQVTPwsJiwoQJJiYmIpGouLiY6DhAW2hRGS0uLo6Pjzc2Nvb29h4/fjz+Y4ph2IkTJzZs2GBhYREZGRkREXH06NF+2gHhGhsb3d3dQ0NDiQ6iFg4ODoGBgQihvXv3Ep0FaAstKqNRUVEHDx4UCoWFhYV//PHH22+/jRDicrkCgcDPzw9/j5+fX1FRUT/tgHB//PHH+PHjDQ0NiQ6iFnZ2dt9//72Dg8OJEydgdS+A06JbTNu3b4+IiNi9ezdCaMOGDSNHjkQI8fl8Go1mbGyMv8fc3JzP5/fT3ltRUVF8fDxCyMvLC1/KIjIykkwjw7XN6dOnIyIiiE6hRsOHD//4448/++yzwsLCnj/kQCeUlJTcvn0bIVRbW1tZWYkQqq+vV36z2nI0imFYeHj4P/7xD3xeu1OnTuFX8U1NTTEME4lE+Nt4PB4+2V1f7b3Z2dlNmzZt2rRpEyZMwOcjsLKy0sg+6anGxsawsDCiU6hXYmKiXC5PT09//Pgx0VmAAuzs7PAiEBISgpcFa2tr5TergqPRU6dO9fUSPs5uILhcbmVl5cqVK42MjDw8PGJjY69cubJ06VJ3d3c6nV5QUIBPx1JYWIgvd95Xe2/29vZxcXEK7xV4IzKZrL6+fvDgwUQHUS8Gg+Hi4nLs2DF7e3tY6k6H2Nra2travtxy69YtFWxXrjSvvg18I11dXebm5rt27RKLxVVVVSNGjNi0aRP+UkJCwsyZMzs6Oq5cuWJiYvL06dP+21+xYsUKJXcQDNzhw4ddXV3b2tqIDqJ2+HhDT09PooMApaikPqjgaPTp06fKb4TJZGZlZa1Zs+bzzz83MTGJiYnpmSR0586diYmJLBbLzs4uPT3dy8ur/3ZAoLKysuHDh+vDjJwjRoygUChcLrekpASmA9dzWnSLKTw8/LWzkFlaWmZlZQ28HRCopKRk0KBB+lBGbW1tly9fvnfv3suXL0MZ1XPacosJkMPVq1cdHR1JM1tzPwwNDXft2kWj0eC5D6BFR6NqUlVV9cUXXyCE5s6diw+iAuqDYRhZB973ZmxsvHLlSnKsN6VvcnJyLly4gBCqq6tTfmvkL6MuLi74ZVZYFlTdpFKpiYnJ2LFjiQ6iOREREf/5z39aW1tVMm4GaExQUNDw4cMRQp9++qnyWyP/yReVSrW2tra2tu4Zqw/U5MKFC1ZWVnQ6neggmhMQEODt7X337l2igwDF0Ol0vCxQqSo4lCR/GQUaU1FRgf+F1x8ODg4eHh6nTp1SyXLnQEdBGQWqIRaLd+3aNXr0aKKDaFpYWNipU6f++OMPooMAwkAZBaqRkZEhlUotLS2JDqJpixYt6u7uvnTpUmdnJ9FZADGgjAIVKCkpuXbtmqenJz4vl16xsrIaPHjwgQMHYAl7vUX+MlpeXp6cnJycnAynXeqTlZWFz2/g6upKdBYCxMbGPn/+PDc3l+ggYKA4HA5eFqqrq5XfGvkHPHl6ev74449EpyC5lpYWLpe7aNEildz31Dne3t6mpqY3b97seYIZaLmIiAh8Osfk5GTlt0b+o1GgARiG8fl80k/s1JdJkyZNmjQpPz9fJpMRnQUQAMooUIHu7m4nJyfSTzPaF09PzwkTJjQ2Nl6+fJnoLIAAUEaBsng8XlFR0bx588zMzIjOQphJkyY5Ozt//fXXRAcBBIAyCpRVWVlZVVW1aNEiooMQKTg4eNmyZfX19TDsSQ9BGQVKEYvF12d5W6sAACAASURBVK9fNzU11c979D0MDQ29vLxqamr27dtHdBagaeS/r4oPeEIIJSYmBgcHEx2HbK5du7Z//35fX19zc3OisxDMz8/PyMjowIEDH3/8MYVCIToO6A+Hw8FXP4IBTwMCA57Uqqur68mTJ5s3byY6CPE8PT3d3Nzy8vJevHhhbW1tZGREdCLQJzIPeNq7d6+HhwedTvfx8SkrK8MbOzo6YmNjmUwmm80+fvx4z5v7ageahC/OOnXqVKKDaIW4uDhzc/PTp0/DUHy9okVHoxkZGVu2bPn5558DAgLKy8stLCzw9pSUFAzDGhoaHjx4EBMTExAQgC+71Fc70KTs7GwqlQqTEOJcXV3feuuttWvXTps2bdy4cUTHARqiRUejGzdu/Pe//x0WFmZlZTVq1Ch7e3uEEIZhJ06c2LBhg4WFRWRkZERExNGjR/tpBxqWl5enDysvDVBoaOiwYcNaW1tfvHhBdBagOdpSRru6uoqKip48ecJisVxdXdeuXYs/EMLlcgUCgZ+fH/42Pz+/oqKiftqBJonF4pKSkvfff5/oINqCxWI5OTkNHjxYJpPdvHmT6DhAQ7SljOIroty7d6+oqOjatWvHjh376aefEEJ8Pp9Go/WcM5qbm/P5/H7agcZ0dHSUlpbK5fK1a9cSnUWLvPfee3v27GlqasrMzCQ6C9AQbbk2ymAwEEKrV6+2sbGxsbFZunTppUuXVqxYYWpqimGYSCTCKyaPx8OXVOqrvbecnBwfHx+EkJeXF36hIDk5Gda2U95ff/117tw5R0dHPZxjtB8uLi7W1taNjY3l5eVEZwGvunXrFn71r66ujsvlov+rPErSljLq4uLy2t9Gd3d3Op2OT8KGECosLPT39++nvbcJEybAgCd1EAqF2dnZc+bMITqI1mEymba2tmVlZQKBQCW/pUBVQkJCQkJCXm4h1YAnCoWyePHib775pr29vbKy8sCBA1FRUQghGo02b968tLQ0Ho/H4XA4HE5cXFw/7UBjMAwrKyvTn+WUFTJo0KAXL15kZ2cLBAKiswC105YyihDaunWro6PjoEGDJkyYEBcXl5iYiLfv3LmTSqWyWKykpKT09PSeUU19tQPNaG5uFolEcEb/WikpKaGhoZcvX7516xbMnkd62nJSjxBiMBgHDx48ePDgK+2WlpZZWVm9399XO9AAsVj8008/hYWFwVOPrzVx4sSMjIzMzMy6ujovLy82m010IqBGWnQ0CnTIzZs3//rrrw8//NDb25voLNqIwWAMGzasra3t+vXre/bsIToOUC/yl1GZTCYQCAQCgVQqJToLeRQWFtbV1QUEBDg4OBCdRUuNGTNm0KBBIpHowIED8LOnbSQSCV4WVHLJRYtO6tWkqqrqiy++QAgtWrQoMDCQ6Dgk8eLFC2dnZ6ih/XB2dp49e/ahQ4eamppOnjy5YMECohOB/8nJyfntt9/Q/41YVxL5y6i7u/u3335LdAqy4fP53t7eBgbkP5t5Y4MHDw4PD58xY0ZYWNiRI0egjGqV0NBQfJCJSgY8kb+MApUTCARXrlyZPn060UG0GoVCGTt2LJ1Od3Z2hgmfyA2OJoDCmpqa+Hw+LCb8t6ysrBgMRkpKSnd3d0NDA9FxgLpAGQUKa2lpcXNzc3JyIjqIbpgyZUpwcHBGRgaGYURnAWoBZRQoLD8/38bGhugUOsPX19ff3//s2bM5OTlEZwFqQf4y2t3dXVxcXFxc3NHRQXQWksjPz4d79ANHoVCioqJu3Lhx69YtorOA/2ppacHLglAoVH5r5L/F1N7ezuFwEEJTpkzpmVEfKOPu3bvLli0jOoUuGTt2rFwuh+frtUdNTc3t27cRQjweT/mtkb+MOjk5rVq1iugU5CGTyRobG1+ZJgf0D18xrKKigugg4L8CAgICAgIQQsXFxcpvjfwn9UC1GhoaXF1dhw4dSnQQHbNw4cLHjx8fPnyY6CBA9aCMAsWUlJTY2toSnUL3rFu3TiAQ7N27l+ggQPWgjALFXL16ta+FBkA/jIyM3NzcCgsLa2tric4CVAzKKFBAVVXVoUOHnJ2diQ6ik8aMGSMSieCJJvIh/y2mtra2CxcuIIRGjhwJI8aVdOPGDR6PN3fuXKKD6KRly5bV19dzOJzo6GhDQ0Oi4+i1ysrKkpIShJBKxkGSv4zKZDKJRIL/g+gsOu/rr7+2sLCABQHfjJeX1xdffDFhwgQMw3bt2sVkMolOpL96yoJcLld+a+QvozY2NtHR0USnIIOurq7i4uKJEyfCxE5vxsDAwM/Pj06n79+/f9SoUfHx8VBJieLh4eHh4YEQwk9VlaR1vw/FxcXGxsZLlizpaeno6IiNjcVH3h0/fvxv24Ga7N69Wy6XT5gwgeggum3ZsmXGxsa7du26fv060VmAamhdGf3HP/4RHBz8cktKSgqGYQ0NDfv27UtKSiotLe2/HajJ6dOnGQzG6NGjiQ6i2+Lj4z/55JPa2tpz584RnQWohnaV0SNHjjg7O798vINh2IkTJzZs2GBhYREZGRkREXH06NF+2oGaFBQUtLe3z549e/z48URn0W34FVIKhfL777/X19cTHQeogBaVUR6Pt3Hjxu3bt7/cyOVyBQKBn58f/qWfn19RUVE/7UBN8vPz29vbfX19YW4nJVEoFFNT07Fjx3K53GXLllVXVxOdCChLi24xrV+/ftmyZa+MSeLz+TQazdjYGP/S3Nycz+f3097bkydP5s2bhxAaNmyYvb09Qmj69Olubm7q2xFSunnzZmtra1xcHNFBSGLVqlWXLl26cOHCb7/9BnM+aExRUdGNGzcQQrW1tVwuFyHU1NSk/Ga1pYwWFBRwOJzHjx+/0m5qaophmEgkwismj8fDH6Hpq703FouFl1EWi2VmZoYQguMpRRUXFx86dMjDw8PFxYXoLCTh7e3t4+NTVFSUlpY2fPjwiIgIohPpBUdHR3xWnfb2dryAnjhxQvnNaksZvXXrVmVlpaurK0KIz+dLpdLS0tJ79+65u7vT6fSCgoKgoCCEUGFhob+/P0Kor/bebG1t3333XQ3uCgnl5OSIRKLZs2cTHYQ8LCwsMjIyfH1929rakpOT7969i58qAbWytra2trZ+ueX3339XfrPacm00MTGxoqIiLy8vLy/vvffei46Ozs7ORgjRaLR58+alpaXxeDwOh8PhcPDzyr7agcpJJJKamhpHR8fExESis5CHlZWVt7d3YGAglUqtrq4+c+ZMd3c30aHAG9KWMspkMln/x8TEhMFg9EwjtHPnTiqVymKxkpKS0tPTvby8+m8HqvX9999fv379X//6l7u7O9FZyObnn3+ePHmyiYnJn3/++d133xEdB7whikqehdJmycnJP/74I9EpdFVVVdWkSZNaW1vr6urMzc2JjkNC58+f/+233w4fPmxqanr37l2BQDBixIjGxkZYpkUzVFIftOVoFGgVmUwmEAjq6+s/++yzmpoaBoMBNVRNgoODw8PDTU1Nm5qa5s6d++2334rF4nPnznV1dREdDQyUttxiUh8ul4uvqB4XF/fWW28RHUc3yOXyu3fv5uTknDhxwsDAAJYMUR87O7tJkyb5+vreuXPn8ePHeXl5kZGRv/zyC5vNDg8PJzodad28eRO/+6Ka6V/lZLd06dLOzs7Ozk6xWEx0Fp3R3d2dlJRkZGSEEAoJCcnMzCQ6EcllZWUlJCTg9wN8fHwQQitXriQ6FJmJRCK8LCxbtkz5rZH/aNTQ0BBma1fU4cOHDx06JJfLEUKrVq2CUY3qNmvWrNDQUJlMlp+fX1hYSKPRfv/9d5lMBpNpqYmRkRF+lKCS7zD5yyhQSEVFRWtr665du0QiUXBwcFVVlaurK1wYVTcDAwMrK6uFCxcGBQUdP37cz8/v6NGjL168oFKpNjY2FAqF6ICgP1BGwf8IBIILFy6kpaXx+XwKhbJw4cLhw4fDkiEaExoaymazAwICTExMurq69u7d29LS8sEHHwwbNozoaKA/UEbB/3A4nF9//bW5udnU1HT58uUJCQlmZmZwXqkxNBrNy8sLHwFtb2+flpZmamrq5+c3ePBgAwMDGo1GdEDwelBGwX81Nzd/9tln+LStPj4+aWlpFhYWRIfSX4sWLbp37979+/e3bNnCZDINDAzmz59PdCjweuQvo+Xl5cnJyQihxMTEVyaEBi/7/fffKysrmUxmV1fXyJEjoYYSKyAg4KuvvlqyZEllZeXq1avt7OxGjRqFr3sBlMfhcE6dOoUQUslEheQvo56envAU09+Sy+VZWVk0Gi0gIMDa2nrNmjVEJwIoJCRk7ty5O3furK+vF4lE9+7dgzKqKhEREfj4E/wYS0lw2QsghNAvv/xy7ty5YcOGrVu37u2334b5WLUBhUJZvnz57t273d3dMQz7+uuv8/LyiA4FXoP8R6Pgb5WUlPzzn//EMCw+Pn7y5MmwfrL28PDwMDQ0dHV13bdv35kzZ9asWfP99997enoSnQv8/0AZ1XdCofDdd99taWmZPHnyiBEjEEJWVlZEhwL/4+TkNHjwYBcXl6ampvv37//444/btm3j8/lw8Vp7wEm9vtu+ffvTp08RQkuWLMEfQwRaxdjYmEKhjBgxIiMjIykpKTMzs6CgoKSkhOhc4H+gjOq1nJyc3NxcCoUSGhoaGhr6ysTgQKs4OjpaWFg0NTXNnj377NmzRMcB/0P++UYjIiLwa0kw4OkVra2tgwYNQgiFhISkpqaOHTsWHjrUcmfOnMnPzz969Cifz7969So83fTGXh7wdPHiRWU3p/zsJlpuxYoVREfQUqtXr0YIubu7X7x4USaTER0H/L3W1laZTFZZWYkQ2rJlC9FxyEAl9UFbTuqFQuGKFSvc3NxMTExGjRqFL4KK6+joiI2NZTKZbDb7+PHjf9sO+icWi58/f97a2nr+/HmEUFRUVGRkJByH6gQrKysKhcJmsydPnszhcOrr64lOBBDSnmujGIaZmJhkZWXV1NTMmTMnOjq6tbUVfyklJQXDsIaGhn379iUlJeFPK/bTDvpXXl6+fv36rVu3VlRU0Gi08ePHw1PzOmfx4sWlpaWRkZEdHR0ymYzoOHpP+QNadTAzM8vJyZHL5UKhkMFg5Obm4u3R0dHr16/vp703OKl/xcOHD6lUKkKIzWavXbu2oqKC6ERAYc3Nzd9//z1CKD09PT8/n+g4OoxUJ/UvKykpwTAMv3zO5XIFAoGfnx/+kp+fX1FRUT/t4G89fPhQIpG4u7t7e3vPnz8f1vvURTY2Nu+//76fn9+ZM2du374NB6TE0rrh911dXfHx8evWrbOxsUEI8fl8Go1mbGyMv2pubs7n8/tpB/2Ty+W//PLL8uXL2Wy2l5cXDBTVXRQK5Ycffli4cKGlpaWzs3NUVBQ+nTvQPO0qoxiGvfPOO4GBgampqXiLqakphmEikQivmDweD18RpK/23nJycvBi4eXlZW9vjxBKTk7Wz+cdi4uLT5069ejRo4sXLz59+nT06NFEJwJKGT9+/IIFC/bu3Xv16tVTp06NGzeO6ETa7tatW0ePHkUI1dXVcblchBCDwVB+s1pURsVi8ezZs+3t7ffs2dPT6O7uTqfTCwoKgoKCEEKFhYX+/v79tPc2YcIEmOEJt3379gcPHiQkJFhYWAQGBhIdB6jAl19+mZeXd/v27cOHD48ZM4ZCoVAoFLlcDkMvXiskJOSVZW5JNcOTVCpdsGCBVCrdu3evSCQSCoX45R4ajTZv3ry0tDQej8fhcDgcTlxcXD/toC8ymSw3N7euru67775DCOF3mYCuo1Kpx44dCw4OPnDgwNSpU+fMmbNlyxb8OAtojLb8LlVWVmZmZiKEzMzM8JaMjIzZs2cjhHbu3JmYmMhisezs7NLT0/ElFvppB691/vz58vJyKysruIJGMpaWlikpKW5ubjk5OTKZ7Pbt23Q6PSEhwcjICNYi1AxtKaOenp7yPh5LtbS0zMrKGng76E0qlW7btk0kEuF/mQDJREdHu7u7+/r6JiQkhIeH79ix4/r161KpNDMzs6ury9bWluiAJKctZRSo1d27d+/fv+/g4PDpp58SnQWoha+v7+DBg+l0+o8//nj+/Pkvv/zSwMAgJiZm4sSJH3/8sUpupIC+kL+MSqVSHo+HEGIwGPp5Pvvs2TMulztmzJjVq1fjc5EAUsIviI0ePXr06NG+vr6XLl06cuTIjRs37Ozsli9fjmEYrC3aA78BgxCSSqXKb438ZbSmpmbTpk0IoQULFgQEBBAdR3Pw0WDd3d0HDx7EBxhOnz6d6FBAQ+bOnevg4GBmZnb58uU9e/YkJSXt37//gw8+IDqXtrh///65c+cQQg0NDSrYnPIPQmk5vX0Y9MqVK+3t7QcPHjQzM6NSqXv27CE6EdAomUzW3d3d1dW1aNEiNzc3T0/PwsJCokNpHdI+DAqUx+PxNm7cuHv37rNnz3Z2dgYFBcHclPqGQqEwGAwmk/nrr7/OnTuXwWAEBQV9++23ROciIfKf1OunioqKnJycu3fvMplMV1fXDRs2TJo0iehQgBgUCmXixIlhYWEXL1789NNPHR0dZ86cyWAw6urqXF1diU5HBlBGyenYsWPTpk1rb2+PiooaMWKEt7c3PNaiz4YOHcpkMr29vcvLy5cvX/7TTz+NGTOmtrYWyqhKQBklodbW1mPHjv3www/nzp1bunSpmZkZk8kkOhQgEpvNlsvlVCr14MGDubm577//fmpqKiy9pSrkL6NdXV35+fkIocGDB1taWhIdR+1EItH+/fuXL18+ZcoUd3d3BwcHohMB4vUM9bOxsXn77be/+uqrlJSUiIgIBoMxZcoUqVTa3t6uV6P0m5qa8LUDBAKB8lsjfxnl8Xh3795FCDGZTH0oo48fP960adOdO3eMjY1hHjzQm6Gh4fz58+fMmfPzzz9/9NFHUVFRIpFIJBKlpqa6uLgQnU5D6uvr8bKgkgk2yV9GHR0d33//faJTqF1bW5tcLjc2Nr527ZqJiQkUUNA/Q0PDpUuXTpw48cKFC9XV1c3NzXFxcWw2e//+/frwlIq/vz8+Jxx+qqok8pdRPXH58uU7d+5QqdSampqTJ0/CDSUwEF5eXj1z+uTl5X3yyScBAQHz589PTU3tmcwX/C0YN0oGV65cuXTp0tGjR0+dOjVkyJAJEyYQnQjoHjc3t6ysrA8++GDPnj2jR4/29PTE144FfwuORskgMzPz5s2beXl5IpHIzc2N6DhAJ1lYWCCEPvjgg8TERAqFkpGRcfDgwS+//DI+Pn7u3Ll2dnZEB9ReUEZ1DF4rg4ODpVJpR0dHU1NTdnb2+fPnORwOjAEEKkGn0xFCixYtWrRoUVFR0fjx40+fPp2cnPzuu+/CWtyvRf4y2tLS8ttvvyGEAgMDdehGZElJibe3d8+X7e3txsbGRUVFH3/8sUgkWrJkSU5OzoULF/z8/Nra2vLz862srAhMC8jKx8fn9u3bP/zww5EjRzZu3Dh8+PCAgIA5c+ZYWVnx+XzdnTCsoqLiyZMnCKH29nblt0b+MmpgYIBfKdetP6T37t1jsVh4cezo6Pjoo4/w//j4+HgbG5uNGzfGxsY+fPjQw8Ojs7OzZ8kAAFRuxIgR+PJoXC63paXlzp07EydOtLa2lkgk8+fPZ7PZw4YNe+utt3TrflRPWVDJzViKvI8550kjOTlZt5a06+7urqqq+uc//2lra2ttbS2TyXJycpqamtasWTNz5szBgwcjhLhcrpubG9yOB4R4+vTpsGHDDh48KJfLm5qanjx5kpubm56ePn78eKKjKUwl9UGHj0Y7OjoSExMvXrxob2+/devW+fPnE51IKXl5eWw2+4cffvj222+DgoKCg4OLi4vNzMzEYvHixYuDg4NfXj7X3d2dwKhAz+GzhS1ZsqSn5cKFC2vWrPHy8vL29g4PDw8ICCgtLa2vr3d0dCwoKGCz2cHBwYTFVT8dLqMpKSkYhjU0NDx48CAmJiYgIEBXVrW7cOGCo6NjaWlpfn7+s2fPnj17hu8Ig8FwdXWtrq5++RF4mUxWVlamK7sG9NP06dOnTZv2+PHjmzdvbtiwobi42NHRkUajyWSyoKCgc+fOnTt3LjAwcOHChR4eHvb29iT7edbVk3oMw6ysrG7duoWvUz9r1qyRI0du3Lix9zsJOamXSCRNTU0VFRVyufzhw4fHjx/HMMzc3Bw/Jb9z546FhYWTk1NCQoKPj8+gQYPu37/v7OwsEok8PT1NTU01nBYAdWtubm5vb8/MzGxoaGhsbCwpKTEyMjI0NHzrrbfCwsI8PDwsLCwIOcfS65N6LpcrEAj8/PzwL/38/IqKitTxQR0dHTKZzMTEpPcVdPxJZB6P19jYWFNTU1tb29XVVV9fn5+fX1lZyWQyx4wZY21t7eLicvnyZQsLC5lM9uLFCwzDqFSqi4tLR0cHPlIPIfT222+rIzwAWsLW1tbW1nbNmjU9LRKJRCQS3b9/v6io6PLly3l5eY6OjjY2NsbGxlQqFb9rik9OZmJi4uLigk+yY2BggP/WGBgYODg4aMlSfbpaRvl8Po1G6ylt5ubmfU0xUFBQEBMTgxDy9va2t7dHCM2cOdPDw6OvLf/nP//Jzs5ub2+Xy+VSqdTa2trQ0LCrq0skEslkspdv98tkMkNDQysrKzc3NyqVOmTIEEdHx+DgYHxm3N73fwwNDR0dHXu+7KmhAOghKpVKpVLDwsLCwsLwls7OzsrKSrFYLJFIOjs7EUI8Hg9fB+XOnTvNzc3opRUqZTLZ8+fPOzs7qVQqhUKh0Wj4pbD33nsvLi6urw8tLCy8evUqQqi2traiogIh1NbWpoJ9UX4ThDA1NcUwrOexXx6P19e5sIuLC34t3N7eHn9P/3PHLVu2bPHixTAMEwANMzMz6zm/VJRQKMSnvOv/KqWrq+vUqVMRQh0dHS0tLQihQ4cOvdknvkxXy6i7uzudTi8oKMCvjRYWFuLztfRmbW2t0IqYTCYTJjkGQLfQ6XT84av+WVpavjJb5tmzZ5X/dF0akf4yGo02b968tLQ0Ho/H4XA4HE4/R/IAAKA+ulpGEUI7d+6kUqksFispKSk9PZ1kQygAALpCh8uopaVlVlZWd3d3dXW1CsfenzhxQlWb0kX6vPsCgUAlp3g6qrGx8caNG0Sn0Ek6XEbV5Pr160RHIJI+775QKLx37x7RKQjT2tqqkqng9ZCu3mIauIqKiv/3//4fQig+Ph6/HwUA0HPXr18/c+YMQujZs2fKb438R6Nubm7bt2/fvn37W2+9NZD3//XXX4p+xK+//kqaLhrYfZlMduTIEYW6SCSSY8eOKdRFKBSePHlSoS4CgaC0tFShLp2dnfhv48C1t7dnZ2cr1KWlpUXRiehfvHhx6dIlhbo0NTUp+r9fV1d37do1hbpUV1ffvHlToS6VlZU5OTkKdSkvL+//xCIkJAQvCyqZPJP8ZdTAwIBGo9FotAFOlNfY2KjoR+BLDJKjiwZ2XyaTKXruLJFIHjx4oFAXsVicm5urUBd8ZgOFuggEgkePHinUpaurKy8vT6EufD6/oKBAoS48Hg+fT3PgOjs7Fd391tbWkpIShbq0tLQo+rfqDep7Y2NjeXl5P28wNDTEy4JKpkkjfxkFAAC10tWpSQZuy5YtNTU1A39/eXm5p6enQh/R3Nxsa2tLji7auftyuby1tdXGxmbgXWQyWVtbm0JdJBJJbW0tm80eeBf82USFnnl7gy4SiYTP578ybvxvu3R1dSn0wHF3d3dbW5uzs/PAu4hEIqFQaG5urlAXDMMUmmgcf15RoS5CoVAikQxklh9vb2/83okyyF9GAQBAreCkHgAAlAJlFAAAlAJlFAAAlAJlFAAAlAJlFAAAlAJlFAAAlAJl9L86OjpiY2OZTCabzT5+/DjRcdRlx44d/v7+VCo1NTX15fa+dp9M3xahULhixQo3NzcTE5NRo0a9PJuRPuw+QmjFihVOTk40Go3NZm/atKmnXU92X43kQC6Xy+VLliyJiopqb2+/fPkyk8l8+vQp0YnU4uTJk9nZ2TExMWvXrn25va/dJ9O3pb29/aOPPnr06FFzc/O2bdvMzMxaWlrwl/Rh9+Vy+b1792pra3k83u3bt83NzS9fvoy368nuqw+UUblcLhcKhQwGIzc3F/8yOjp6/fr1xEZSq4SEhJfLaF+7T+5vi5mZWU5Ojlz/dl8ikeTl5dnb2+fl5cn1b/fVAU7qEdLgcs3aqa/dJ/G3paSkBMOwYcOGIT3b/aVLl1Kp1KCgoM2bN+PLl+nV7qsJlFGEFFmumZT62n2yflu6urri4+PXrVuHP3GvV7ufnp6OTyH40Ucf4bNG6dXuqwmUUYReWq4Z/7Kf5ZpJqa/dJ+W3BcOwd955JzAwsOcmm17tPkKIRqPFxMRMmTIlMzMT6d/uqwOUUYReWq4Z/7KwsNDHx4fYSJrU1+6T79siFotnz55tb2+/Z8+enkb92f2XiUQifCFx/dx9FSP64qy2SEhImDlzZkdHx5UrV0xMTMh6U1IsFgsEgkWLFn322WcCgUAikeDtfe0+mb4tEokkNjZ22rRpnZ2dAoFAIBBIpVL8JX3Y/ba2tl27dlVWVjY1NR04cIDBYBQWFuIv6cPuqxWU0f9qa2uLiYlhMBiDBg06duwY0XHUZc2aNS//Ed29ezfe3tfuk+nbUlZW9soxREZGBv6SPuw+j8eLjIy0srLCh81euHCh5yV92H21gvlGAQBAKXBtFAAAlAJlFAAAlAJlFAAAlAJlFAAAlAJlFAAAlAJlFAAAlAJlFOgRFov19OlTolMAsoEyCgAASoEyCvTF8uXLm5qawsPD2Wx2dnY20XEAecBTTECPsFisGzdu4NOMAqAqcDQKAABKgTIKAABKgTIK9IiBAfzAA9WDnyqgR+zt7cvLy4lOAcgGyijQI198gXhKYgAAAHZJREFU8UVycrKVlVVWVhbRWQB5wJ16AABQChyNAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUqCMAgCAUv4/9kjdoQAC+eAAAAAASUVORK5CYII=",
"prompt_number": 6,
"text": [
"FramedPlot(...)"
]
}
],
"prompt_number": 6
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 'Individual based' simulations"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"using SimJulia\n",
"using Distributions\n",
"using DataFrames;"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"prompt_number": 6
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"type SIRPerson\n",
" state::Char\n",
"end\n",
"\n",
"type SIRModel\n",
" sim::Simulation\n",
" parray::Array{Process}\n",
" beta::Float64\n",
" c::Float64\n",
" gamma::Float64\n",
" ta::Array{Float64}\n",
" Sa::Array{Int64}\n",
" Ia::Array{Int64}\n",
" Ra::Array{Int64}\n",
" allIndividuals::Array{SIRPerson}\n",
"end"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"prompt_number": 7
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function increment(a::Array{Int64})\n",
" push!(a,a[length(a)]+1)\n",
"end\n",
"\n",
"function decrement(a::Array{Int64})\n",
" push!(a,a[length(a)]-1)\n",
"end\n",
"\n",
"function carryover(a::Array{Int64})\n",
" push!(a,a[length(a)])\n",
"end"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 8,
"text": [
"carryover (generic function with 1 method)"
]
}
],
"prompt_number": 8
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function SIRModel(heapsize::Int64,beta::Float64,c::Float64,gamma::Float64,S::Int64,I::Int64,R::Int64)\n",
" sim = Simulation(uint(heapsize))\n",
" N=S+I+R\n",
" parray = [Process(sim,string(i)) for i in 1:N]\n",
" states = [fill!(Array(Char,S),'S'),\n",
" fill!(Array(Char,I),'I'),\n",
" fill!(Array(Char,R),'R')]\n",
" allIndividuals=[SIRPerson(state) for state in states]\n",
" ta=Array(Float64,0)\n",
" push!(ta,0.0)\n",
" Sa=Array(Int64,0)\n",
" push!(Sa,S)\n",
" Ia=Array(Int64,0)\n",
" push!(Ia,I)\n",
" Ra=Array(Int64,0)\n",
" push!(Ra,R)\n",
" SIRModel(sim,parray,beta,c,gamma,ta,Sa,Ia,Ra,allIndividuals)\n",
"end"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 9,
"text": [
"SIRModel (constructor with 2 methods)"
]
}
],
"prompt_number": 9
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Main loop"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function live(p::Process,individual::SIRPerson,s::SIRModel)\n",
" while individual.state=='S'\n",
" # Wait until next contact\n",
" hold(p,rand(Exponential(1/s.c)))\n",
" # Choose random alter\n",
" alter=individual\n",
" while alter==individual\n",
" N=length(s.allIndividuals)\n",
" index=rand(DiscreteUniform(1,N))\n",
" alter=s.allIndividuals[index]\n",
" end\n",
" # If alter is infected\n",
" if alter.state=='I'\n",
" infect = rand(Uniform(0,1))\n",
" if infect < s.beta\n",
" individual.state='I'\n",
" push!(s.ta,now(p))\n",
" decrement(s.Sa)\n",
" increment(s.Ia)\n",
" carryover(s.Ra)\n",
" end\n",
" end\n",
" end\n",
" if individual.state=='I'\n",
" # Wait until recovery\n",
" hold(p,rand(Exponential(1/s.gamma)))\n",
" individual.state='R'\n",
" push!(s.ta,now(p))\n",
" carryover(s.Sa)\n",
" decrement(s.Ia)\n",
" increment(s.Ra)\n",
" end\n",
"end"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 10,
"text": [
"live (generic function with 1 method)"
]
}
],
"prompt_number": 10
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function activate(s::SIRModel)\n",
" [SimJulia.activate(s.parray[i],0.0,live,s.allIndividuals[i],s) for i in 1:length(s.parray)]\n",
"end\n",
"\n",
"function run(s::SIRModel,tf::Float64)\n",
" SimJulia.run(s.sim,tf)\n",
"end\n",
"\n",
"function out(s::SIRModel)\n",
" result = DataFrame(t=s.ta,S=s.Sa,I=s.Ia,R=s.Ra)\n",
" return(result)\n",
"end"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 11,
"text": [
"out (generic function with 1 method)"
]
}
],
"prompt_number": 11
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Running the model"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Set parameters\n",
"# All Float64 so use decimal points\n",
"beta = 0.1\n",
"c = 1.0\n",
"gamma = 0.05\n",
"# Set initial conditions\n",
"# All Int64\n",
"S0 = 99\n",
"I0 = 1\n",
"R0 = 0\n",
"N0 = S0+I0+R0;"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"prompt_number": 11
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Initialise and run\n",
"srand(1)\n",
"sirmodel = SIRModel(N0,beta,c,gamma,S0,I0,R0);\n",
"activate(sirmodel);\n",
"@elapsed run(sirmodel,1000.0)"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 12,
"text": [
"0.352742373"
]
}
],
"prompt_number": 12
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Results"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"r=out(sirmodel) # my function\n",
"wp = plot(r[:,\"t\"],r[:,\"I\"])\n",
"setattr(wp,\"xlabel\",\"t\")\n",
"setattr(wp,\"ylabel\",\"I\")\n",
"wp"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEsCAIAAADfNCTgAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVxU1f8/8DPADAww7CCC7AgKIiqKGCJpmWkkYqh8XOkjWVqppX0q9ZNmZmVpapbWA0wNQcNdVCgX3EBQVFZ3QUYGkGWYYYYBZvv9cb7Nbz4sw8DcmXsvvJ9/+BjuHC7vGYcX5y7nHIZSqUQAAAB6y4jsAgAAgN4gRgEAQCcQowAAoBOIUQAA0AnEKAAA6ARiFAAAdAIxCgAAOoEYBQAAnUCMAgCATiBGAQBAJxCjAACgE4hRAADQCcQoAADoxITsAnrpwoULYrGY7CoAAPTm7+/v7++v405o2RstKir68ccfr1+/fv369ebm5l27dnX7LZ9++imXy9Xc5sKFC/fu3dPcprS09MKFC5rbNDY2Jicnd1uSNmUvXbqUkP0cPHiQz+drbnPp0qWSkhLNbe7fv//3339rbsPj8f7zn/90W5I2ZRPVJjU1tb6+XnObPXv2nDp1SnObR48eZWZmam4jEon27dvXbUlEvTRtPiGHDx+ura3V3ObKlSsFBQWa2zx+/DgjI0Nzm+bm5r1793ZbElEfbG3apKWl1dTUtNv47NkzHCDLly/vdg/dU9JQYWHh8uXLMzIyMjIyRCLR4sWLu/2WcePGlZSUaG6zd+/ea9euaW5z9erV33//XXOb6urqNWvWdFuSNmX7+fkRsp9169bxeDzNbfbv33/58mXNbbKzsxMTEzW3efDgQVhYWLclaVM2UW3Wr1/P5XI1t1m6dOm2bds0t8nNzf311181t6mvr//kk0+6LYmol6bNJ+TLL7989uyZ5jYpKSnnz5/X3ObmzZu7d+/W3KaxsfHjjz/utiSiPtjatNm0adPTp0/bbXzw4AEOkFmzZnW7h27R9aDe29t7ypQpZFcBAKAlPz8/Pz8/hFBaWhoBu9M9iQ2vsLBw+/btqi9bW1u7/ZZPP/20pqZGcxuZTCaXyzW3kcvlMpms2x/X1tbWbRttyn777bcJ2Y829RD18mtra7XpjmlTNlFt2traFAqF5jYHDhy4cOGC5jYG/t8n6hMilUq7ffla/u9LpdJuf5w2L5+oD7Y2bTS/fG36s92ia29UHYvF6raNkVH3Z4GNjY0J2Q9CiMlkdtuGqLK12Y829RD48okqm6g22rx8IyMjBoPRbZtu96PljyPqpWlTkolJ97/mWv7va/PjtHy3u21D1FukzcvXES0vMQEAAHUwlDRc0q6oqOjixYsrVqzQ/luam5vZbHa33Q2qEYlElpaWZFfRM0qlUiKRmJubk11Iz7S2thoZGWnTk6IUOn5CEJXKTkhISExM1HEnfeGgXhu0+63GKPJR6xEGg0HHd9vU1JTsEnqDjp8QRNuyuwIH9f2aQCAguwQAaA9itF87ceKESCQiuwoA6I2uB/Xp6elFRUUIoY0bN7q4uJBdDl01NjaePXt29uzZZBcCgEEdOXIEj8iqqKjQfW90jdGoqKgeXWICnZJKpX/88QfEKOhvYmNjY2NjEUIJCQm67w0O6vs1qVT6999/NzQ0kF0IADQGMdqvyWSykSNHHj58mOxCAKAxiNF+TSqVLliwIDU1lexCAKAxiNF+TSqVent7Gxsbl5WVkV0LAHQFMdqvSaVSExOTf/3rX4cOHSK7FgDoCmK0X5PJZEwmc9asWX/++SfZtQBAVxCj/ZpUKmUymba2tu7u7t1Ofg4A6BTEaL+GYxQhNHfu3JSUFLLLAYCW6Hr7PYxiIgQ+N4oQioqK+uGHH8guBwADgVFMCMEoJoK0trbimY0sLCykUinZ5QBgIDCKCRBGFaMAgF6DGO3X2traIEYB0BHEaL+m3hul3dIAAFAExGi/1traqloUzNjYWC6Xk1sPAHRETozKZLL333/f19fX3Nw8MDBQfUy3WCyeP38+h8NxdnbesmULKeX1H+q9UXNz8+bmZnLrAYCOyLlSL5VKjYyMUlNTfXx8MjIy4uPjBw8ePHr0aITQqlWrKioqysvLuVzu5MmT/fz8ZsyYQUqR/UFbW5uqN8pms5ubmzkcDrklAUA75PRG2Wz2Tz/9NGbMGDs7u7lz5w4fPjwvLw8hJJPJkpOT169fb29vP2LEiISEhH379pFSYT+hVCpVK4ZDbxSA3iH/vlE+n3///v3g4GCEUFlZmVgsHjFiBH4qODj4yJEjnX5XaWkpfsrT0xP3pwYPHsxmsw1VdR8EMQr6PB6PV1dXhxCqq6vDs5Xz+Xzdd0tyjMpksvnz50dHR4eHhyOE8PJq1tbW+Flra+umpqZOv7GwsBC/HWPGjMHL+To7O0OMam/q1KkHDx5Uv6ZkZmbW0tJCYkkA6NujR4/w3BFPnz7lcrkIodraWt13S2aMyuXyBQsWKBSKpKQkvAWvXi0QCOzt7fGDrk7VxcXFwSgmXVRXVz969MjV1VW1BWIU9HmRkZGRkZHqW+g9ikkuly9atKiuru748eOqi8VeXl4WFhaqqYYKCgoCAwPJqrBvE4vFd+7cGTp0qGqLmZlZa2sriSUBQFPkxKhCoYiPjy8vL8fTXLa0tMhkMoSQiYnJvHnzNm7cWF9fX1hYmJSUFB8fT0qFfV7HGDU1NYUYBaAXyIlRHo+XnJx8/fp1Ozs7NpvNZrM3bNiAn9q6daurq6uHh8fkyZNXrVoFdzvpSXNz8+3btyFGAdAdOedGBw0apFQqO33K0tLy4MGDBq6nH2pubi4pKRkyZIhqC8QoAL0Dg0H7I6lU2tbWxuFwHBwcVBvhEhMAvQMx2h+JxWKEkPoRPYLeKAC9BTHaH+Hb7DvGKPRGAegF8kcx9U52djbuOi1evBjfZAq0JxaLLS0t8SQGKmZmZoWFhUKh0MrKiqzCADCMy5cv5+bmIoRqamp03xtdY9Tf3x9fxIepNHqhqqoqPj5+8eLF6htNTU137NgxderUKVOmkFUYAIYRGBg4cOBAhBBe0k1HdI1Re3t7Pz8/squgq4qKCjc3t3Yb8SAImHIU9AcODg74+iohqz/AudH+iMvluru7t9sIMQpA70CM9kcQowAQCGK0P4KDegAIBDHaH/F4PBcXl3YbIUYB6B2I0f5ILpcbGxu32wgxCkDvQIz2O13dGQoxCkDvQIz2OxUVFR2vLyGIUQB6i673jaanp+P7Zjdu3NjxNB/QgMfjqU96r2JtbT158mSIUdAfHDlyJCMjAyFUUVGh+97oGqNRUVGwiEjvCIVC1WpX6iwtLWNjY2F2EtAfxMbGxsbGIrovIgLIIhKJLCwsOn0KJnkCoBcgRvsdPC9Jp09BjALQCxCj/Q70RgEgFsRovyMWiyFGASAQxGi/Awf1ABALYrTfgYN6AIgFMdovNDY2Pn/+HD/WcFDPYrHa2toMWBcAfQFdY1Qul7e1tbW1tXW1UDNQ19DQcOLECfxYJBJ1dVAPi4OCfoLYAKHr7fd//fXX48ePEUL//e9/8WIAQAOBQJCenv7BBx8oFIpr166ZmZl12szCwgIvGgpA33b69Om//voLIaQ6StMFXWN06tSpMIpJe3w+/+LFi3w+39zcvL6+vuP0TtiAAQMIWeELAIqbMWMGXswNRjEBbQkEAiaT+ddff8lkMoQQg8HotJmlpaVIJDJsaQDQHsRov8Dn82fNmnX8+HGpVIq6jlHNTwEAOgUx2i/w+fxJkyaVlJTgU59GRl3+v5uZmUkkEgOWBgDtQYz2CwKBwMbGJjIy8vz580hjlxNOjwLQUxCj/UJjY6ONjU10dPSxY8eQxt4oxCgAPQUx2i/gGJ0wYUJubi7S2Bt1cnJ68eKFAUsDgPYgRvsFHKOmpqaRkZEIDuoBIBRd7xt9+PDhyZMnEUKTJk3icDhkl0N1OEYRQtHR0X/++afmg/oHDx4YsDQASFBaWvro0SOEkEAg0H1vNO6NGhkZGRkZwQ062hCJRPiPzbRp0xYvXmxubt5VS1dX18rKSgOWBgAJGAwGDhBC9kbX3qifn9+bb75JdhW0oVQq8d8bGxub1157TUNLDw+PZ8+eGaouAMgxdOjQoUOHIoTwQa2OaNwbBVpSKBTqf3WnTp2qobGzs3N1dbX+iwKg74AY7fv4fL6tra3qS82nko2MjGCNZQB6BGK076urq3N0dNS+va2tbWNjo/7qAaCPIS1Gd+3aFRISwmKx4uLi1LfPmDGDoQZPpQF0UVtb26MYdXd3r6io0F89APQxpMWoi4vLF198ER8f3/Gp7du3N/3DxISuF8Goo66uzsHBQfv2bm5uXC5X9SXM+QSAZqTF6MyZM6Ojo+3s7Do+ZWpqavkPwxfW9/SiN6q6WK9UKq9evaqfugDoI6h4bnTz5s1ubm6RkZFnz54lu5a+oLa2tke9UXd3d1VvtLy8nJDpwQHowyh3yLxkyRInJycOh5OZmRkTE3Px4sXw8PCOzXbu3Pnbb78hhEaPHo1vJv/ss888PDwMXS4d1NXVjRs3Tvv2bm5uFRUVf/3112uvvXb79m1YVgT0GYcPH87KykIIlZWV4b6Cqamp7rulXIxOmzYNP/D398/NzT18+HCnMbp06dJly5YhhJhMJt4CZ1G70tMr9d7e3k+ePFmxYkVRUdGtW7fg1AroM2bNmjVz5kyEkEKhwDf2ffDBB7rvltLRw2KxurpSz2QyNYxoBOp6eomJzWY3NzcjhIRC4a1bt8aMGaO30gAwqI4DQAkZD0rauVGZTNbS0iKXyxUKRUtLC45LiUSSnJxcWVnZ0NCQkpKSmpqKl50Cuqirq7O3t+/Rt1hZWQkEAqFQmJ+fjyMVANAV0mJ006ZNbDb7hx9+SEtLY7PZ7733Ht6emJg4bNiwQYMGbdmyJSkpSfMAcKANmUymOvWhJU9Pz+rq6sLCQnNzczg3CoBmpB3Ub9iwYcOGDe02stlsfAIYkMvHx0epVF64cOHll1+GGAVAMyre8AQI1NLSYmZm1tPv8vT0RAhdvnx5/PjxcFAPgGYQo31cfX19j64vYd7e3gihoqKiCRMmUKE3qlAoyC4BgC5BjPZxPR3ChOEYdXNz8/LyIiVGb968qf7l9evXDV8DAFqCGO3j6uvre3qZHiHk4uLCYrFGjRplZmbW0tKij8I0O3HihPqX6enphq8BAC3RNUbT09MTEhISEhJ4PB7ZtVBaL+52QggZGxsvWbIkJCSEwWAolUp9FKZZZmamSCRSXYQ8ffo0zIIKCHTkyBEcIIRMZkbp2+81iIqKWrFiBdlV0EBjY6O1tXUvvjEmJoasWQqbmpru3r378OHDAwcO4CQtLy+/ffs2DAQARImNjY2NjUUIJSQk6L43uvZGgZaEQqGVlVUvvpHD4YwcOZLwerSRn58vl8sLCgpwR1gqlUokksuXL5NSDADdghjt43odo0OGDOnFtSlCZGdnDx8+XBWjAoHAx8fn0qVLpBQDQLcgRvu4Xh/Ua16ySa+uX78+e/bs27dv4xgVCoXBwcGVlZVSqZSskgDQAGK0j+t1b5QsSqWyoqIiLCysqKhI1Ru1srIaNWrU7du3ya4OgE5AjPZxuseoUqmUSCRE1dOtBw8e+Pn5ubi4NDY2qmLU2tp64sSJcFwPqAlitI8TCAQ2Nja67MHExKSkpISoerp17dq18PBwFxcXhJB6jL766qsXL140WBkAaA9itI/TvTfKZDKLioqIqqdbOTk5L730krW1tYWFhXqMDhw4sLa2trW11WCVAKAlut43Wl9f/+DBA4SQl5cXi8UiuxzqkslkxsbGuuyBxWKVlpYSVU+37ty5M2LECITQwIED8TUlfG4UITRmzJj8/PyXXnrJYMWAvqq2trahoQEhRMgfZrr2Rh8+fHjq1KlTp041NTWRXQulMRgMHffAZDKLi4sJKaZbZWVlDg4OeEqqrKwse3v71tbW5uZmvJBJYGAg/tsJgI5KS0txgDQ2Nuq+N7r2RseNGwejmAzDxMTEYAf1SUlJb7/9Nn7s6uqKr863trbidceMjIxgqidAiMjIyMjISIQQIX+Y6dobBQbDZDIrKyvr6+v1/YOkUunJkyffeust1Zbw8PDr16+rxygpA/wB0AxiFHQDL7lqgIv1Z8+efeWVV9TPdIeFhd24caO1tRVvJGueFAA0gxjt43TPHSaT6eXlVVhYSEg9Gvz6669Lly5V3+Lv7//gwYPW1lZ8thRiFFATxCjohomJydixY/XdGy0vLxeLxf7+/uobGQyGl5fX/fv3oTcKqAxitI/T/Uq9iYnJmDFj9B2jBw4cWLRoUcftYWFhOTk5+NwoxCigJohR6tJ99Q6lUklIjLq6uvL5fB3301F1dTV+IJfL09LS5syZ07HNuHHjxGIxxCigMohR6jp27JiOqSGTyfAFIl0wmUxLS0t3d3cul6vjrto5fvw4fnDu3LkJEyZYWFh0bBMWFsZisVQxCjc8AQqia4zm5eXt2LFjx44deChCn1RWVtZuZbeeIiRGnZ2d7ezsAgICiD2ub25uVo2R379/v+p20XbYbPawYcPwuVG44QkQ5fr16zhAamtrdd8bXW+/9/T0nDRpEkKo0y5M39DW1paWlhYaGtrrPVRUVLi5uelYxpdffokQunHjRk1NjY67UpeZmamaOOrx48chISFdtXzppZdwb5TFYhEy5gQAHx8fPMI4Ly9P973RtTfq5OQUFBQUFBSEf8H6JKlUeujQIV36XwUFBcOHDyekGA6HQ+y422PHjqkeaz6HGx4ejv+XLS0tRSIRgTWAfsvZ2RkHCJvN1n1vdI3R/qCtrU2hUNy6davXeyguLg4KCiKkGGIjTCaTZWVladkYnx5FCJmbmzc3NxNVAwBEgRilLqlUGhcXl5aW1us9FBUVUTNGr169quEovh1PT0/cG7WwsND97gUACAcxSl1tbW0xMTFnz57t9XF9eXm5p6cnIcUQG6PHjh2LiYnBj6VSqZbXweCgHlATxCh1SaVSS0vLkSNH9u4suFAo5HA4ut83ihEYYUqlMisrC8+vgxASiURaLp8HvVFATRCj1NXW1sZisWbNmtW743oCj+gRoZeY8vPzg4ODmUwm/lIkEuHpRLsFMQqoCWKUuqRSKZPJnDJlSmZmZi+O64mNUQJ7oydPnoyOjlZ9CTEK6A5ilLpwb9TU1HTEiBG9uA+fsjGakZExdepU1ZdNTU1wUA9ojZjb7xMTE7t6KiEhgZAf0U56ejqekn3jxo14Fck+IykpafHixeif3ihCKDo6+syZMz29D//evXtDhw4lqipLS0tCIqyqqsra2hrvDS+1JBaLtRxDYWJiIpfLda8BgCNHjmRkZCCEKioqdN8bMTG6a9eurp7SU4xGRUX11UVEPv/8cxyjqqGcw4YNO3LkSE/309zcTOAQL2NjY0IirKamZtCgQQihAQMG8Hg8hJBUKtV+UUKirpiBfi42NjY2NhYRFFDExOjdu3cJ2Q8QCARCoRA/VigURkZGCCF3d/de/M0kZEC9OkIiTH3BZw8Pj2fPnsnlcu3XLsXLMeG3BQCKgI8jtVRVVak6fap86d3onR7FkzZMTU11X41WKBRaW1vjx6GhoXl5eT2Ke1tbW31M2QeALiBGqYXH43WMUYSQqalpS0sLeXUhhJC9vb3uC9vhu1nx49DQ0Js3b8pkMu3j3tHRkZApeQAgEMQotfB4PKVSiW9vUo9RNze358+fk1oasre3r6ur03En6gf1oaGht27dksvl2vdGIUYBBZEWo7t27QoJCWGxWHFxcerbxWLx/PnzORyOs7Pzli1byCqPLFVVVQgh3CFVPwno7u7+7NkzMitDyMHBgZDeqCpG7ezs+Hx+W1ub9r1RBwcH3aMcAGKRNt+oi4vLF198cebMGdUVFWzVqlUVFRXl5eVcLnfy5Ml+fn4zZswgq0jDwzGKTxeq90Z7d5WJWHZ2drrHaFNTkypGEUJDhgwpKSmxtbXV8tshRgEFkdYbnTlzZnR0tJ2dnfpGmUyWnJy8fv16e3v7ESNGJCQk7Nu3j6QCyVFVVWVnZ4d7o+ox6uHhQXqMtuuNtrW19WInAoFAPUZDQ0NzcnLgoB7QGrXOjZaVlYnF4hEjRuAvg4OD9b0gJdXweLxBgwZ1PKjH9waRWlr7S0xnz57txela9YN69M9Vph5dYoLeKKAaai0igocbqm6Isba27mo6jO3bt//yyy8IoZCQEHNzc4TQunXriJoUjkRCodDe3l4ulwsEgidPnqhu1Rw4cCA+3ieRo6Pj1atXVV/W1dV9//33O3bsaNfsl19+CQkJGTt2bKc7qaysHDhwoOrLwMBAkUikfW/Uzc2trKysh4UD8H9SUlLwCmBlZWW4E4DTQ0fUilE8RYVAILC3t8cPuhptvXLlyr46igmPF2psbFQffMlms0m/4cnV1VW9+9nU1JSYmLhmzZoBAwaoN8vKynJ1de10D21tbfX19ert8f+4mZmZljV4e3s/ffoU7sAHvTN37ty5c+eqbyFkFBO1PoteXl4WFhYFBQX4y4KCgsDAQHJLMjwTExOZTNbY2Kg+aoiQW991NGjQoMrKStWXTU1N06ZN+/7779s1q6iokMlkne7hxo0bYWFh6ltMTExMTU17tKBWcHCw6hMCABWQFqMymaylpUUulysUipaWFvyLZ2JiMm/evI0bN9bX1xcWFiYlJcXHx5NVoeG1traampri3qhAIKBajJqZman3iJuamhISEi5evNjumg+Xy8VzjnR08eLFiRMntttobm6ufW8UIRQZGan9Ok4AGABpMbpp0yY2m/3DDz+kpaWx2ez33nsPb9+6daurq6uHh8fkyZNXrVrVr+52EggE1tbWnfZG8Vhy7XdF+EhQDNeGHwuFQltb2w8++GDnzp2qBm1tbdXV1V31RrOysl5++eV2G83NzXvUG42MjLx8+XLP6gZAn0g7N7phw4YNGzZ03G5paXnw4EGDl0MJOEblcjnujepy+g/PVUpgbRi+0uXm5ob+mSd0wYIFY8aMWbVqlY2NDUKIx+MpFIpOY7S5ubmpqcnR0bHd9p72Rn19fcvLy+H0KKAO+CBSCI7RTg/qe0o1VymxXF1dVadH8az1TCbz3//+92+//YY3crlcJpOpHqOqeftzc3PHjRvXcZ9sNrtHvVGEUGBgYHFxcW9eAAB6ADFKpnb37jQ2NtrY2HR6UI/hCTq1oafeqHqMqmatf+edd1JSUiQSCUKIy+W6ubmpzo2KRCLV7a6XLl3qeGIU9fygHsHpUUAxEKNkunbtmvoMeHiEj6o32vGgFc/XLRaLu7qHVHWPlJ5idNCgQap7nlQxymaz586dm5SUhBDicrne3t6q3ujt27cFAgF+rL4aqLqeHtQjhF5++WWIUUAddI3RJ0+eZGRkZGRk0Hrh8vr6+lu3bqm+7Pag/ty5cwghoVC4e/fujnurqKhQjZ3V00G9+j1P6lexli5dmpSUJJVKcYw2NTXhDmleXh5OdrFYLJFIHBwcOu6zF71RPz+/R48e9eiaGwDqHj58iAOk3ZwevUPXGG1raxOJRCKRiNa/S0KhMCcnR/XlnTt3hgwZ4ujo+OLFC6FQmJ+fr97YyMgoJyenvr5eLBb/+uuvHe9/+uOPPy5duoQf83g89cFCRHF1deVyuQihtrY29ZjmcDgjR468e/duXV3dwIED//vf/+JVT/Ly8vDfuUuXLk2YMKHTfZ4+fdrPz6+nlTg5OfXh+ZtPnTr18ccfk11FX9ba2kpggNA1RocOHYpXU1EfoE07QqHwxo0bqi+vXLkSGRk5dOjQe/fuicXidsPULCwsampqrl+/LhaL6+vrU1NT2+3tyJEjpaWl+JJOUVHRsGHDCC9YNdHUnTt3Ro4cqf5UUFBQUVFRS0sLh8NRKBTV1dVILUZPnTo1ffp0AivhcDhdDRTuA+RyOeljf/u2oKAgHCD4DhMd0TVG+wb13mhJSYmvr6+pqakqRtstSMfhcORy+dWrV8Vi8fTp0/fs2aP+bE5OTmBg4LBhw0pLSxFCRUVFw4cPJ7xg1TLL2dnZ7S67Dx8+vKioSCKR4BOm1dXVNTU1z549E4vFCoUiOzs7PDycwEr6fIzW1NSQXQXQFsQomfA1JXy9Pj09PSoqCiGEY7Tjup4cDicgICA7O1ssFvv6+g4ePPjKlSuqZ/fv379o0SLVrenELlKvjslkSqXSjjGKe6MSiQQfH9TU1Ny8edPe3l4kEuXl5YWEhBC7vp4q0PskuVyOu/OAFiBGySQUCl977TXcIT179uy0adMQQjY2No2Nje1OPiKEOByOi4sL7qdYWFgsX75cNbtSS0tLdnb2K6+8oroTCN87pY+a8T1PDx48aHdC08nJqba2Fh/UI4Sqq6vz8vIiIyPFYjHhR/QIeqOASiBGySQUCqdMmXLjxo36+nqZTKaa+sje3r7j9RMOh8PhcMLCwi5cuGBhYTFmzJgXL16Ul5cjhE6dOvXmm28aGRkFBgaWlpZyuVy8Frw+uLu7X7lyxcPDo+ONBAMGDOByuebm5ra2tnV1dXl5eRMnThSJRJmZmVOmTCG2jD4fo3w+n/RZFICWIEbJJBaLIyMj8/Lyzp07N3XqVNX2oUOHqtYHVeFwOJaWlhERERkZGfh4/4MPPvj5558RQvv371+4cCFCiMFg+Pv7//nnn/o4MYq5ubkdOnSo0/FIw4cPF4lETCYTT5laWVk5ZMiQgoICZ2dnPCEegfr8Qb1SqXzx4gXZhQCtQIySzMrKSiaTpaWl4ROj2NChQ9udGEX/xOj48ePxQT1CaObMmRkZGY8fP25qavL398fNIiMjd+/eracTowghd3f38+fPdxqjQUFBZmZmTCbTzs6OwWD4+PhYWFicO3cuOjqa8DJ60RulUeziP6JwXE8XEKPkGz169N27d4ODg1VbNMTogAED/Pz8cOeOyWT+61//euutt+bNm6dqFhkZ+eTJE73GqFKpDA0N7fhUUFAQm81mMpm2trbOzs6hoaGWlpZtbW3qfyGIoiFGu1okKj8/H9/0SgVCoVDDn7LMZfkAAB/7SURBVAG5XG5nZwcxShd0jdHMzMyVK1euXLmSdhc0c3JympqacnJyVHMUhYWFTZs2Tf1U49ChQzuubeDr64vDMSIiQhWyS5Ys4XK5c+bMUTULCgpydnZWdU4J5+7uPnz48I4pjxAKCAiwtLS0sbEZP368s7PzmDFj8DlcFxcXwsvo6qBeqVQuXLjw/v37x48fb/dUVVXV+fPnCa+kd06cOJGZmdnVs3w+39XVlXafbRo5ffo0DpBerCfWEbUWEdHexIkTly5dighaSsWQoqKisrKy1q9fv3fvXpwv48aNazd93MCBAzuuw/HKK6/gB+ox6uDgkJqaqn5R3sjI6J133tHHSFBVbV2NRzIzMwsKCvL29l63bt3atWtHjx4tlUr1cUSPELKysup0GN+mTZtOnjzp6el54MCB0NBQ9beRx+PdvHnz7bff1kc9PZWcnOzo6BgbG9vxqYcPHyYnJ69evRpWndKfKVOm4Ilyli9frvve6NobZbFYlpaWlpaWtJt0UiwWP3z48MWLF0+fPvX29kYI+fn5qSJSZfTo0V3tQT1GEUIdL4K/++67xNXbnrGxcae//JjqYD8iIsLW1tbS0pLwW50wa2vrjjF64sSJa9euLVu2rLy8fMqUKZ999pn6s9XV1VlZWaqJ+0hUW1srEAiKi4s7Tswqk8neeeedXbt2DRs2jDqnIPoeYgOEZhnUByiVyrt377548eLJkyc4RhkMRscpjjTEqJeXl4+Pj4Yf0dWKckR56aWXunpKFaP4T725ubk+xqQihKysrFRzR2HFxcUbN25MTU319fV99OhRaGioUCg8c+aMqgGeZ6CwsFAf9fTI4cOH4+LiIiIi1AdQYF988cWrr74aGRlJhSW1gZYgRg1NoVDgKTyePHmiIQ1HjRqlYSfkziSgYT5pVdk9nbSpp9od1NfV1S1cuDAlJcXOzs7T0/Px48dMJvOnn35au3YtnggVIcTj8RYsWECF06NpaWlxcXHTp08/deqU+vYLFy7cvHlz7dq1CCEHB4f6+nqSCgQ9AzFqaEZGRnfv3pVKpfn5+bg32ilra2tDVkWUTi896YN6jEql0ri4uK+//nrIkCEIIU9PT6FQyGQy3d3d582bt3nzZtyMz+fPnj37woULhqmwK+Xl5cbGxgMHDpw4caL6rKl1dXUff/zx77//rn6YSYVTEKBbEKOG5unpiSexz8vL0xCjQDP1xfVWrlw5bdo01fgFT09PBoOBL7KtWLHi7Nmzjx49wk+5urrW1NS0W7hUfblTAzh69OisWbMQQkwmc/DgwSUlJXj7smXLPv/8c/XhZ46Oju1WXQXUBDFqaObm5l5eXgghNpvNZrPJLof2du7cqVAo1GfnZLPZzs7OOEZZLNbu3buXLFmimnhw3Lhx2dnZ6nu4c+eOIQvGR/T4cUxMzLFjxxBCu3fvtrOzU23H4PQoXUCMGpRSqWQwGPi2UBymQBdXrlw5fvz49u3b22338PBQTSgVGhrq7++/ZcsWPI/1K6+80u64Xn0BAn0rLS11dHS0tbXFX77xxhsZGRmFhYW///57x1fh5uaGZ3cFFAcxalCqGB00aJDmq+2gW6amph999NGhQ4c6Xs7y9PRUv3P266+//u2335ydnVFny9wbsjd65MiR2bNnq760trY2NTVdtGjRnj17Ot6t4ebmBvc80QJdb79PT08vKipCCG3cuFEfg2T0RCwWm5mZxcXFsVgsfSw5168MHDhww4YNqmmx1CUkJKh39u3t7b/55hu8ipSdnR2e9cPJyQk/e+vWrXXr1m3atMkANZ85c6bdrQJvvvmmQqHo9MYMFxeX+/fvG6CqfujIkSN4gUhi+vtKGiosLNy+fTvZVfTGpUuX3n//faVSef369eTkZLLLoTcul6t9Y4VCUVRUhB/v2bPnu+++w49lMhmeS6WyspL4Ev9XbW1tZGRku418Pl+hUHTa/tKlSytXrtR3Vf3c4sWLdd8JHNQbVHZ2Nr533cnJCQ7qddSjOVUZDIZqIMDcuXMPHTqkVCoRQi9evJBKpXw+X3XFXH+ysrJefvnldhttbGy6ug/XxMSk3U0FgJogRg3qxo0beIo5iFES4dmv//77b4QQnplCqVQWFxfr++eeP3/+1Vdf1b69+k1dgMogRg2qrKwMn7OzsrJqNx0JMKTFixcnJSUhhHg8Hr5ChZcC1Ku8vLwxY8Zo357JZEKM0gLEqOFoHv0JDCkkJITL5b548aKysnL48OFWVlb37t3T60+srKx0dHTs0RhZ6I3SBV2v1NNRTk5OWFgY2VWA/7N48eLExESRSBQaGsrj8UQikVKp1DBdgI7Onz/fcR4vzSBG6QJ6o4aTl5c3duxYsqsA/ycuLu7PP//kcrkhISEMBkPfN2leunRp0qRJPfoWuMREFxCjvaeaOkhFoVA8f/68q0nLb9261aNTY0CvLCwsxo4de+bMmVGjRjEYjICAAL0e19+5c2fkyJE9+hbojdIFXWO0sbGxrKysrKyMrD/Xd+7cuXTpkvoWuVy+Y8cOmUy2b9++ju2bm5vlcjnhC2QCXSxZsqS5uTkgIMDExCQgIEB/9zw9efLE09PT2Ni4R98Fl5j0h8/n4wDpauWuHqFrjBYXF6ekpKSkpLSbu9dgTp48efPmTfUtmZmZx48fVyqVf/zxR8f2jx490t/6SKB3QkJCoqKimExmenp6REQEniWEcEqlcvXq1QsWLOjpN3I4nMbGRn2UBO7cuYMDpKGhQfe90fUS0/jx41esWEFiAadPn243DHHXrl1lZWUymay0tDQ3N7fdadDKykp9T0oPemH9+vUIoYCAAITQ4MGDz5w588YbbxD7I77++ms3NzcNK690xdbWlpBfctDRpEmT8KnqhIQE3fdG194ouR49ejRw4MCqqirlP7PqPnz4UCqVhoeH379/38HBoeNxPY/HgxilIPWVqL/88sv169crFAoC93/06NHr16//+OOPvft2CwsLDeswA4qAGO2NtLS02NhYPz+/hw8f4i07duxYuXLlmDFjbt68OWHChBs3bjQ3N6t/y/PnzyFGKc7d3X3s2LFHjx4laoelpaUbN248ePBgT8+Kqvj6+j59+pSoeoCeQIz2xunTp6dPnz527Njc3FyEkEAguHHjxrRp00aPHn3z5k0jI6Po6Oh2v41wUE8L69at++677+Ryue67EggE8+bN279/v52dXa934uvr++TJE92LAXpFuRidMWMGQw0Fr1Q+f/7c3Nzczs4O9z0RQgcOHJg7dy6DwRgxYsSdO3cYDMb8+fOTk5PVv4vH49FoQr9+a+DAgREREYcOHdJxP0qlctGiRatWrRoxYoQu+/H19X38+LGOxQB9o1yMIoS2b9/e9A/VHObUceTIkbfeegshNHr06Fu3bikUit9//33JkiUIIWtra1tbWwaD4evrK5PJ1PsR1dXVePZ1QHGfffbZ999/r+Pf7++++87V1XX+/Pk6FgO9UVqgYoyamppa/oPsWjpx8uTJmTNnIoTYbLZSqTx69Gh4eDiHw8HPhoSE4Afx8fHqdz7J5fJenyADhjRgwIBJkyalpqb2eg+XLl06ffr0tm3bdC/Gx8cHx2hrayufz9d9h0AfqBijmzdvdnNzi4yMPHv2LNm1tFdbW9vW1oaXo0AIjRgxYvXq1UuXLlU1wENiEEIxMTEnTpzAl/JbW1thrnsaWbt27ZYtW3p3Y/aTJ0+WL19+7NixHs1C0hVHR8e6ujqEkEAgUF+NGVAK5Q6ZlyxZ4uTkxOFwMjMzY2JiLl68GB4e3rFZWlrajRs3EEKjR4/GKz7GxsYaYOq5U6dOTZ8+XfXl2LFjnzx5gu86xEJCQvAJU0tLy5EjR167di0iIqKqqgpOjNKIvb39G2+8kZKSEh8f36NvlEgkCxcu/Omnnzpd2qR3WCxWa2urQCC4cuVKTEwMUbvtn65du4YXHyorK8NTKFRVVem+W8rF6LRp0/ADf3//3Nzcw4cPdxqjQ4YMwROJ+/j44I6ehYWFAcorLy8PDQ1VfTlu3Lh2vzAhISGqg/eYmJiMjIyIiAgul9ujqdoB6WJiYvbu3dujGJVKpXFxcQkJCR2nuNdFQEBAcXGxQqFoN/gY9IKHhwde1dzT07O+vh4hRMj9bZSLUXUsFqurM/1BQUG6n7/vBalUqr7kpJ+fn5+fn3oDS0tL1aDPsLCwHTt2ILhMT0Ourq54CTwtyWSyOXPmvP7662+//TaxlYSGhubl5Q0ePLiwsLC2thZm+9aFm5ubm5ub+hZCTpVQ69yoRCJJTk6urKxsaGhISUlJTU2dMWMG2UX9D6lUqn6W08jIyMio/Xuo6q46OTk1NDS0tbVxudx2/3mA4gYMGFBTU6N9+5UrVwYEBKifJSfKqFGj8vPzBQKBmZnZ1atXCd8/0B21YhQhlJiYOGzYsEGDBm3ZsiUpKem1114ju6L/0a432in1xXJHjBhRWFgII0Fph8lkaj952IoVK9hstp6WaB45cmRBQUFjY+OUKVMuXryojx8BdEStg3o2m03xy5HaxKj6YVdERMTVq1efP38O50Zph8lktrW1dXuLxebNm5ubm7dv366nMkxNTY2MjCorK6dOnfrrr7/q6acAXVCuN0pNqtGB2sSourCwsJycHBgJSkcDBw7sagZule3bt9+9e3fPnj36W30EITRy5MisrCxnZ2dLS8senWoAhgExqpXa2lr88dWme6Ju6NCh9+7dE4lEhrmRABAIT+KlocHRo0dPnjy5f/9+fQ+sCAkJycnJsba2njBhApwepSCIUa1IJBK8AK9UKu3R+FQGg+Hh4aGaTw/QiOaL9enp6bt27Tpz5gy+gUavQkNDW1pabGxsJk2adOHCBX3/ONBTEKNakUgkhYWFqMOVem1ERETAET0dubi48Hi8Tp+6fPny119/ffz4cTz0Q98CAwPZbLa1tfW4cePwqBNAKdS6xKS99PR0PBph48aNBrglUyKR4PXOenpQjxAaN27co0eP9FMX0CN3d/fDhw933C4SiZYtW3bu3DkbGxvDVGJiYhIcHGxtbW1mZmZjY1NYWDh8+HDD/Oi+6siRIxkZGQihiooK3fdG195oVFRUYmJiYmKiYW5rb25uxqldV1fn4ODQo+8dO3asr6+vfuoCejR58uT6+voTJ060215eXj5q1Ch3d3dDFjN27Fhra2uEUFJS0oIFCwj55e/PYmNjcYAQ8v9I1xg1MIlEUlxcrFQqJRJJT8+FmZqavvnmm3oqDOgPg8FITExct25du8yqqKgw/GCKiRMn4oEe3t7eP//8c0xMDKx2Rx0Qo1qRSCRCofD58+e9+/bAwEBi6wGG4eDgsHXr1nfeeUd9gSZS7gKeMGGC6vH48eM//PDDRYsWETJLP9AdxKhWJBIJk8nMz883zCUFQB1TpkwZM2aM+giliooKDw8PA5dha2ur/mV8fHxoaOi7775r4DJApyBGtdLc3BwUFHT+/HkCJ0ADdLFhw4bz58+rJliiyAwJa9euZTKZW7ZsIbsQADGqHYlEMmbMmIsXL0KM9kMmJiYHDx5cuXIlnlqNlHOjndq5c+eFCxdSUlLILqS/gxjVikQiGTly5IMHD2Casv7Jzc3tiy++wCtu8fn8dofYZGEymWlpaTt27MjOzia7ln4NYlQrzc3Njo6OXl5e0Bvtt9566y17e/tffvmF7EL+h5WV1bFjx5YtWwb3JpMIYlQr+D6nwMBAiNH+bNu2bXv27HFyciK7kP/h6uq6e/fuuXPnCgQCsmvRSktLi0gkUt+C1/OgL7rGaHp6ekJCQkJCQlfD9YhVXFzs5uYWEBBAtV8hYEiWlpb79u3z8fEhu5D2xo0bl5CQsHXrVrIL6Vx5efnp06c3b94cFxc3evRoX1/fU6dOqZ5taWmZOHGigc9LHDlyBAcIIQMZGHScNaOoqOjixYsrVqwwzI+7cOFCYmJiampqSkpKSEiIao0Q0D/1YiSbATQ0NEyfPv3atWtkF4KampqKioqKiooKCgqKi4tFIpGHh0dQUNDw4cOHDx/u4+Oza9cua2tr1TpX3377bVVVVU5OzrVr1wy/gG5CQkJiYqKOO6HrmHqDUSqV69evxyvODxs2DHqjgIIZihCys7NTKpUNDQ12dnaG/LkSiSQ/P7+0tLSkpCQ/P5/P53t7e4eEhAQGBr7//vt+fn4d5+fF82Hjxzwe78SJE9evX//xxx83b968YcMGQxZPFIjRbhw7diwoKMjLywsh5O/vD8vNA8qaNGnSpUuX3nrrLf39CKVSWVpaqgrN8vJye3t7f3//gICAWbNmbdiwQZt7GFgsVktLC368bt26devWGRsbf/TRRxEREbNnz1ZfrpwuIEY1aW1t/eqrr/7++2/8pampKbn1AKDB5MmTU1JSiI1RkUh09+5d3N+8d++eTCYbNGhQQEBASEjIwoUL3d3dezT9LsZisYRCIULo1q1bPB4vKioKIWRsbLx79+5333338uXLHZeJpDiIUU327Nkzc+ZMuFcU0MJLL720cuVKXfYgl8vv37+vOkh/+PCho6NjSEhISEjIkiVLBg8ebGVlpXudLBYLLxf46aefbtu2TbU9ODh47Nixv/3223vvvaf7TzEkiNEuCQSCxMREmCUX0IWJiYmrq+vTp0+9vb21/Jaampq8vDwcmnfv3jU2Ng4KCgoMDAwPD1++fLmeZqHE50aPHz/u5eUVHBys/tTGjRvDw8PffPNNes10DjHape+++27ZsmWwhhKgkcmTJ58/fx6PtuqopaWlpKQEn9YsLS2tq6sbNGhQYGBgQEBAVFTU4MGDDXPqn8ViNTU1ffnll5mZme2eMjc3/+abb1auXJmWlmaASogCMdq5ioqKzMzM3NxcsgsBoAcmTZq0adMmVYzy+fzCwsKioiL8b3Nzs4+PT1BQ0IQJEz788ENvb29SzkKyWKwDBw4sWbKk08Esr7/+emJiYkZGxuuvv2742nqHrjFaXl6OZ9wZO3Ys4ZPXKZXKL7744pNPPunF6XMASBQYGFhaWvr5558XFhZyuVxbW1t8w+a///3vYcOGWVpakl0gQgixWCwjI6P//Oc/XTXYuXPntGnTJkyYoL95KZ88eYJvvG83nqp36BoTTU1NePwSPldNoOLi4vfffz8wMHDOnDnE7hkAfWMwGF999ZWxsfF7771n+ElRtcRisb766isNme7i4jJr1qy9e/d+8MEHeqpBKBQSGCAwiun/EwqFn3zySUlJyc8//9zuzDcAgCi1tbV2dnbGxsYa2tTV1b3yyit37tzR92kHQkYx0ez+LP05cODA+PHjw8PDr169ChkKgP44OjpqzlCEkIODw4QJEzquJ0hNEKPo/v37r776an5+flZW1sKFCxkMBtkVAQDQxx9/vHPnTrKr0Apdz40SorW19Ztvvjlz5sy2bdsiIiLILgcA8P95eXnZ2trm5eWFhoaSXUs3+m9v9K+//ho3bpytrW1OTg5kKAAU9NFHH+3YsYPsKrrXH3ujNTU1H374oZmZ2enTp+k1WAKAfmXChAlr1qyhyBqCGvSv3qhcLt+xY8frr7++cOHCAwcOQIYCQHHLli2j2sItHfWjGL1582ZERASfz8/OzsaTygAAKG727Nnnzp0j5CZ5/aHrQf358+crKysRQqtWrdJmfaT8/PzVq1f/8ssvI0aM0H91AABimJiYzJkzZ/PmzV9//TWBd9GcPXs2KysLIYRjREd0vf3+7Nmz77zzDkLI2tq623vQEEJKpVKpVNJuHkMAQGtr65o1a0pKSvbu3UvUpFMSiUQikSCEVq1a9fvvv+u4N7r2Rs3MzHq0WAKDwYAbQgGgI1NT061bt964cSMqKuqjjz5asGCB7vtks9lsNhshpE0nrFvQOwMA0EBYWFhWVta1a9dmz57N5/PJLud/QIwCAOjBysrq119/nTVr1qRJky5fvkx2Of8fXQ/qAQD906xZs8aOHRsfHx8UFPT9999TYZVJyvVGxWLx/PnzORyOs7Pzli1byC4HAEA57u7uFy5c8Pb2njhx4oMHD8guh3q90VWrVlVUVJSXl3O53MmTJ/v5+c2YMYPsogAA1MJgMFasWBEeHh4XFxcfH798+XISryFTqzcqk8mSk5PXr19vb28/YsSIhISEffv2kV0UAICiRo8enZ2d/fTp02nTplVXV5NVBrVitKysTCwWq+6QDw4OLikpIbckAACVsdnsHTt2rFy58vXXX09PTyelBmod1OMhX9bW1vhLa2vrpqamTltu374dj7QNCQnBC7asW7fO09PTQIUCAKhkypQpAQEB33777bRp0zSMsklJSbl48SJCqKys7Pnz5wghQpZ7olaM4uVZBAKBvb09fsDhcDptuXLlSsIXEQEA0Jebm9vPP/+suc3cuXPnzp2rviUhIUH3H02tg3ovLy8LC4uCggL8ZUFBQWBgILklAQCAZtSKURMTk3nz5m3cuLG+vr6wsDApKSk+Pp7sogAAQBNqxShCaOvWra6urh4eHpMnT161ahXc7QQAoDjKxailpeXBgwdFIlFNTc2nn35K1G6//fbb+vp6ovZmMKtXrya7hB5raGjYvHkz2VX02LFjx7Kzs8muosfo+AlBtC27K5SLUT0RCoVyuZzsKnqssbGR7BJ6TKFQCIVCsqvoMYlE0traSnYVPUbHTwiibdld6S8xCgAAetIXYnTNmjXdtsnNze32D+DRo0dv376tuU1+fv6xY8c0t6mvr9+2bVu3JWlT9q1btwjZz48//lhbW6u5zfHjx2/evKm5zZ07d9LS0jS3aWxszM3N7bYkbcomqs3OnTu7HeJy+/bt0tJSzW0KCwsPHTqkuY1AIPjuu++6LYmol6bNJ2TXrl1VVVWa26Snp+fk5GhuU1xcnJqaqrmNSCT65ptvui2JqA+2Nm1++eUXfIuo/lDrvlHtiUQiHo+HEHJycnrx4kW37SUSiUwm09xGKBTiCbE176fbI1aZTKbNeVgtyyZkPw0NDd2+/Kampm5/XEtLS7cvX6FQtLS0dFuSNmUT1Uably+RSNra2jS30ebly+Vyov73ifqENDQ0SKVSzW2amprwNMYatLS0CAQCzW3kcnldXV23JRH1wdamDZ/P7/jyhUIhHuzT7TujDbrG6NWrV2tqahBCixYtampqKisr09y+tbX1+fPnmj8odXV1HA5H866qqqpqa2s1t6mrq2tsbOy2JG3KlkqlhOynsbGxoqJCc7rV1taamppq3hWPx6urq9Pc5vnz562trYSUTVSbxsZGLper+RdGKBQ2NDTo/vIbGxsFAoHBXpo2nxA+n8/lcjVfG6itrVUoFLq//KamJm1ePlEfbG3a4JffbuPVq1dxj7jbfro2aLkWU2tr644dO1S/FUKh0MrKSvO31NTUODg4aF4wQCKRmJiYMJlMDW2kUqlcLjczM9PQRqFQNDc34xFZGmhTNo/H63bxGW32IxKJzM3NNS9F1dLSYmxs3O3Ll8lkmv8a4f5It+sMalM2UW3EYjGbzdb88gUCgYmJiYWFhYY22rx8Av/3ifqEiMViMzMzzR/+lpYWIyMjzXN3avPylUqlWCzu9uUT9cHW8n9fw8sfOnTozJkzNe+hW7SMUQAAoI6+cIkJAABIBDEKAAA6gRgFAACdQIwCAIBOIEYBAEAnEKMAAKATiFEAANBJ349RGi18P2PGDIYa1fhFqr2EXbt2hYSEsFisuLg49e1d1UmR+rsqm8pvu0wme//99319fc3NzQMDA9WHtFP23dZQM5Xfal3QdTCo9ui18P327dsXL16MH5uY/N//DtVegouLyxdffHHmzJl2A8y7qpMi9XdVNqLw2y6VSo2MjFJTU318fDIyMuLj4wcPHjx69GgN5ZFetoaaEYXfap0o+zSpVGphYXH+/Hn85WeffRYdHU1uSRpER0fv3r273UbKvoRPP/10zpw5qi+7qpNq9bcrW0mrt3306NE///yzkj7vtlKtZiWt3uoe6eMH9bRb+H7z5s1ubm6RkZFnz57FW+jyErqqkxb10+Jt5/P59+/fDw4ORvR5t9VrxmjxVvdUHz+o137heypYsmSJk5MTh8PJzMyMiYm5ePFieHg4XV5CV3VSv35avO0ymWz+/PnR0dHh4eGIJu92u5oRTd7qXujjMar9wvdUMG3aNPzA398/Nzf38OHD4eHhdHkJXdVJ/fqp/7bL5fIFCxYoFIqkpCS8hfrvdseaER3e6t7p4wf19F34nsVi4euYdHkJXdVJl/oxCr7tcrl80aJFdXV1x48fNzU1xRsp/m53WnM7FHyre4/sk7N6t2TJksjIyLq6uoKCAkdHx+PHj5NdUeeam5v/+OOP58+f19fXHzx40NTUNDMzEz9FtZcglUolEsnq1atnzZolkUikUqnmOilSf6dlU/xtl8vl8+fPDw8Pb2hokEgktHi3u6qZ4m+1Lvp+jDY1Nc2dO9fCwsLJyenbb78lu5wuNTc3R0ZG2tjYsNns4ODg5ORk1VNUewnr169X/0u8ePFivL2rOilSf6dlU/xt7zht+9q1azWXR3rZXdVM8bdaFzBtMwAA6KSPnxsFAAB9gxgFAACdQIwCAIBOIEYBAEAnEKMAAKATiFEAANAJxCjomzw9PW/cuEF2FaBfgBgFAACdQIyCPmjZsmXPnz+PiYnx9PTcvXs32eWAPg5GMYG+ydPT89ChQ2FhYWQXAvo+6I0CAIBOIEYBAEAnEKOgbzIygs82MBD4qIG+acCAAY8fPya7CtAvQIyCvmnNmjVr1661sbHZvn072bWAPg6u1AMAgE6gNwoAADqBGAUAAJ1AjAIAgE4gRgEAQCcQowAAoBOIUQAA0AnEKAAA6ARiFAAAdAIxCgAAOoEYBQAAnUCMAgCATiBGAQBAJxCjAACgE4hRAADQCcQoAADo5P8Btw8XyFQwDfcAAAAASUVORK5CYII=",
"prompt_number": 13,
"text": [
"FramedPlot(...)"
]
}
],
"prompt_number": 13
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Simple MCMC"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function norm(n::Int64, alpha::Float64)\n",
" vec=Array(Float64, n)\n",
" x=0\n",
" vec[1]=x\n",
" for i in 2:n\n",
" innov = rand(Uniform(-alpha,alpha))\n",
" can = x + innov\n",
" aprob = min(1, pdf(Normal(),can)/pdf(Normal(),x))\n",
" u = rand()\n",
" if (u < aprob)\n",
" x=can\n",
" end\n",
" vec[i]= x\n",
" end\n",
" vec\n",
"end"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 14,
"text": [
"norm (generic function with 1 method)"
]
}
],
"prompt_number": 14
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"@elapsed norm(1000000,1.0)"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 15,
"text": [
"0.371801304"
]
}
],
"prompt_number": 15
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Parallelism"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Spawn"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"np=100\n",
"rrefs={}\n",
"for i in 1:np\n",
" push!(rrefs, @spawn sir(0.1/1000,0.0,1000,999,1,0,100) )\n",
"end"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"prompt_number": 16
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Fetch"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"results={}\n",
"while length(rrefs)>0\n",
" push!(results,fetch(pop!(rrefs)))\n",
"end"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"output_type": "stream",
"stream": "stderr",
"text": [
"Warning: ignoring conflicting import of Stats.minmax into DataFrames\n",
"Warning: using FS.rename in module DataFrames conflicts with an existing identifier.\n"
]
}
],
"prompt_number": 17
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Reduce"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x=[r[\"I\"][endof(r[\"I\"])] for r in results]\n",
"x=convert(Array{Float64,1},x);"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 17
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function auc(h::(Range{Float64},Array{Int64,1}))\n",
" freq=convert(Array{Float64,1},h[2])\n",
" e=convert(Array{Float64,1},h[1])\n",
" numbins=length(e)\n",
" deltax=e[2:numbins]-e[1:(numbins-1)]\n",
" sum(deltax.*freq)\n",
"end;"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"prompt_number": 17
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Compare with analytical results"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"xhist=hist(x,100)\n",
"nbins=length(xhist[1])\n",
"xauc=auc(xhist)\n",
"xdf=DataFrame(xmin=xhist[1][1:(nbins-1)],xmax=xhist[1][2:nbins],count=xhist[2])\n",
"xdf[\"dens\"]=xdf[\"count\"]/xauc;"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"prompt_number": 17
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function diagFun(b,N)\n",
" [-(b*(N-n)*n) for n in 0:N]\n",
"end\n",
"\n",
"function lDiagFun(b,N)\n",
" [b*(N-n+1)*(n-1) for n in 1:N]\n",
"end\n",
"\n",
"function Qmatrix(b,N)\n",
" d = diagFun(b,N)\n",
" l = lDiagFun(b,N)\n",
" Q = zeros(N+1,N+1)\n",
" Q[diagind(Q)] = d\n",
" Q[diagind(Q,-1)] = l\n",
" Q\n",
"end\n",
"\n",
"Qmat = Qmatrix(0.1/1000,1000)"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 18,
"text": [
"1001x1001 Array{Float64,2}:\n",
" -0.0 0.0 0.0 0.0 0.0 \u2026 0.0 0.0 0.0 0.0\n",
" 0.0 -0.0999 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0999 -0.1996 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.1996 -0.2991 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.2991 -0.3984 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.3984 \u2026 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 \u2026 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" \u22ee \u22f1 \u22ee \n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 \u2026 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 \u2026 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 -0.2991 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.2991 -0.1996 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.1996 -0.0999 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 \u2026 0.0 0.0 0.0999 -0.0"
]
}
],
"prompt_number": 18
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"pinit = zeros(1001)\n",
"pinit[2] = 1\n",
"pinit;"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"prompt_number": 18
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function siexpm(Q,p,tm)\n",
" expm(Q*tm)*p\n",
"end\n",
"@elapsed aresult = siexpm(Qmat,pinit,100)"
],
"language": "python",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 19,
"text": [
"3.385622524"
]
}
],
"prompt_number": 19
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"p = FramedPlot(title=\"\", xlabel=\"x\", ylabel=\"Density\")\n",
"xhist2=(xhist[1],xhist[2]./xauc)\n",
"add(p, Histogram(xhist2...))\n",
"add(p, Curve(0:1:1000,aresult,color=\"orange\"))\n",
"p"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEsCAIAAADfNCTgAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3de1xUdf4/8PfcZxgYZrh6AQEBQRRWLbTULhamVlSaptUv+yVba7mybma2a9/a2NWf+y3Tyoe6RfVNrLTSvt521SwzLxmQIYaiKBdT7ve5MtffHweJkMsZzgxnZng9/5o5vOdz3ofyxTlnzvkcgcPhIAAA6C8h3w0AAHg3xCgAACeIUQAAThCjAACcIEYBADhBjAIAcIIYBQDgBDEKAMAJYhQAgBPEKAAAJ4hRAABOEKMAAJwgRgEAOEGMAgBwghgFAOAEMQoAwAliFACAE8QoAAAniFEAAE4QowAAnCBGAQA4QYwCAHCCGAUA4AQxCgDACWIUAIATxCgAACeIUQAAThCjAACcIEYBADhBjAIAcIIYBQDgBDEKAMCJmO8G3Ovll1++fPlyaGgoy/qCgoKwsLBhw4axrK+srFQqlYGBgSzrq6urpVJpUFAQy/ra2lqhUBgSEsKy/tKlSzabLSEhgWV9Y2Oj2WweMmQIy/qWlha9Xs/+96PValtaWiIiIljW19XVVVRU3HzzzSzrjUZjTU1NdHQ0y/q2trarV6/GxsayrLdYLOXl5fHx8SzrzWbz999/f8cdd7CsdzgcxcXFo0ePZllPROfOnUtKSnJf/eHDh9PS0jynH3f3X1xcfOjQIfb13XP4tPnz5+/fv599/VNPPbV9+3b29du2bTt16hT7+s8///zo0aPs6/fs2XPo0CH29WvXrn355ZfZ1x89evTzzz9nX3/q1Klt27axry8oKMjOzmZfv3///vnz57OvLy0tffPNN9nXV1ZWrlmzhn19Q0PDq6++yr6+pqbmlltuYV9vtVqXLVvGvt7hcCxdutSt9b/73e/cOr6n9e9sfbdwUA8AwInob3/7G989uFFpaemECRPCw8NZ1ldUVIwdO5b9QWhkZOSwYcMkEgnL+mHDhg0fPlwqlbKsHzJkSGRkpEwmY1nf2NgYFBQ0duxYlvXBwcHR0dFyuZxlfVBQ0MiRIxUKBcv6wMDA2NhYPz8/lvVGo1Gv10+ZMoVlvVKpHDVqlFKpZFmvUCgSExP9/f1Z1kul0tGjRwcEBLCst9vt58+fv/fee1nWC4XC5ORklUrFsp6I3F2fm5v74IMPek4/7u7f2fpu+fi5UX9/f/aZxdSzzywi0mg0TvWjVqudqmd/1pWhUCgEAgH7evaBwvDz82OfiUw/7DOXiKRSqVMtSaXSsLAw9vVisZj931QiEolEQ4cOZV8vFAqd/ZUOHz7cqXr2f+P7V+9s/+7ux939O1vfLRzUAwBwInA4HHz34EabNm266667EhMTWdabzWaRSCQSidzalftYrVYiEou99SDDbrdbrVanDiA8jclkYn+SxAMZjUanDiA8jbP9Z2Zmvv322xxX6q3/3tzEq/8BkzcHKEMoFHr7fwKvzlAi8uoMJZ76x0E9AAAn3r3zAgCDRFNTU1tbG/M6ODiY/eUxA8DHY1Sv169evVqj0SQkJCxZsoTvdgCgn+bNm6dSqYRC4c8///zee+/ddtttXEbLy8vLyckhooqKCu69+XiMKpXKVatWsf+KCQA81ieffCKXy//yl79wHyo1NTU1NZWIMjMzuY+Gc6MAAJwgRgEAOEGMAgBwghgFAOAEMQoAwAliFACAEx5itLa2dubMmQqFIikp6eTJkywL1q5dO2bMGKFQmJ2d3VE5depUwXVz584doA0AAOiEhxhdvHhxcHBwVVXV8uXL58yZYzAY2BTExMS8+eabt956a5fiPXv2GI1Go9H46aefDtAGAAB0MtAxqtVq9+7dm5WVpVarMzIy1Gr14cOH2RTMnz9/xowZN847IJFI5HK5XC73qJvDAGDwGOgYLSkpkUgkHc8US05OLioqcqqgi8WLF8fExDzyyCNlZWVu6hkAoBcDfTOoXq/v/EgAlUql0+mcKujsxRdfjI6OFggE69atmzVrVmFhYZdp1hoaGubPny+XyzUaTUpKChFNmjTp4Ycfdtn2AID3KC8v37RpExHV1NQUFxeTi+aWHOgYVSqVWq22421ra2tcXJxTBZ098MADzIvs7OzQ0NCzZ8/edNNNnQuCgoKys7Pj4+NFIhFz1O/tM3ICQL9FRUVlZWURkd1uN5vNRLRq1Sruww70QX18fLzZbC4tLWXeFhUVjRkzxqmCbjFT1jNzv3cmEAgCAgLUanVAQABzChUxCjBoCQQCJgf8/PzUarVarXbJoy4GOkYDAgLS09OzsrJ0Ol1OTk5jY2NaWhoRvfPOO6dOneqlwGq1mkwm5iETJpPJ4XC0tLR88cUXNTU1zHf6CoWCOWwHABhIPFzwtHnz5urq6pCQkDVr1uzcuZN50uTHH39cUFDQS8GyZcsUCsWRI0eeffZZhUJx4sQJu92+fv36uLi40aNHnz17dv/+/d7+/AMA8EY8HOGGh4cfOHCgy0JmV7SXgo0bN27cuLHLwhMnTrijQwAA9nAzKAAAJ4hRAABOEKMAAJwgRgEAOEGMAgBwghgFAODEx2/pcTgcDQ0NtbW1UqlUrVbz3Q4AeASTydTa2kpEFouF+2g+HqMmk+mzzz4LDQ2NiYl5/PHH+W4HADxCSUnJ7t27iai2tpb7aD4eowqF4tlnn01MTOS7EQDwIMnJycnJyUSUmZnJfTScGwUA4AQxCgDACWIUAIATxCgAACeIUQAAThCjAACcIEYBADhBjAIAcIIYBQDgBDEKAMCJj98MqtfrV69erdFoEhISlixZwnc7AOAR8vLycnJyiKiiooL7aD4eo0qlctWqVbinHgA6S01NTU1NJdxTDwDgCRCjAACcIEYBADhBjAIAcIIYBQDgBDEKAMAJYhQAgBPEKAAAJ4hRAABOEKMAAJz4+M2gNpvtwoULZrM5ICAgJiaG73YAwCM0NzdfuXKFiAwGA/fRfDxGrVbr2bNn6+rqhg8fjhgF8CkOGwlE/ftofX19bm4uEen1eu6N+HiMymSyWbNmYWoSAF9T+Ao1/UR37O3fp+Pi4uLi4oiosLCQey84NwoAXkZIFqo6QOIA0pbw3QsRYhQAvI7KVEjhd1PY7VR3nO9eiBCjAOB1AtvOUPgdpE6mlvN890Lk8+dGAcD3BFiKKXgiOWykvch3L0SIUQDwOiKblqRBRERtDXz3QoSDegDwLv5ig0UUxHcXv4EYBQBvEi6vM0ii29+I5GQz8dkNESFGAcC7hMvrDZLrt9LIw8lUzWs7RIhRAPAuodIGgySq/Y0shEx1vLZDhBgFAO8SIms0iiPb38hDqY3/GPXxb+r1ev3q1as1Gk1CQsKSJUv4bgcAuJIKLTaBvP2NLITa6vsxSF5eXk5ODhFVVFRwb8nHY1SpVK5atQr31AP4CJvRYu+UWtJgMlztxzCpqampqalElJmZyb0pHNQDgPfQXmowa359K1WTpZm/btohRgHAe2hL6ts6x6iGzIhRAAD2dJcbsTcKANB/2ks37I028ddNO8QoAHgPfXmTRf3rW4mKLK38ddMOMQoA3sOqN9slv74ViMlu5a+bdp4So7W1tTNnzlQoFElJSSdPnmRZsHbt2jFjxgiFwuzs7IHtFwAGnN1MAk+8RtNTYnTx4sXBwcFVVVXLly+fM2fOjY/r67YgJibmzTffvPXWW/loGQAGlMB4lZQj+O6iGx4Ro1qtdu/evVlZWWq1OiMjQ61WHz58mE3B/PnzZ8yYoVAoeGocAAaOQHeZ/GP57qIbHhGjJSUlEokkNrb9F5ScnFxUVORUAQD4PIG+nPxveEy6QEQOGx/t/MojTjTo9XqVStXxVqVS6XQ6pwp6Gfndd98NDQ1VqVTx8fFEFBERkZSU5KLGAcC9Pv30U7PZTERVVVUCQzmFziY6f+DAgdLSUiKaOHHiaLE/WXUkCWQzWnNzM/N4+sbGxrKyMiK6erU/95J24RF7o0qlUqvVdrxtbW319/d3qqAXbW1tJpPJZDKZzWaz2Wy18v+9HgCw9I9//IP5l7ts2TKRsYL8Y+fMmRMdHW02m0+dOnXgwAESK8mqZzmaw+FgRuuIBbvdzr1Jj9gbjY+PN5vNpaWlI0eOJKKioqInnnjCqYKeKJXKpUuXYmoSAC+l0Wiefvrp9jeHc0gxJDV1KDOryM6dO69cueJUjGo0mvvvv7/zEt+ZmiQgICA9PT0rK0un0+Xk5DQ2NqalpRHRO++8c+rUqV4KrFYr8/eEeeFwOHjeEgBwH7uNSNB1IXNQzyvXxOjBgwdtNk5neTdv3lxdXR0SErJmzZqdO3f6+fkR0ccff1xQUNBLwbJlyxQKxZEjR5599lmFQnHixAnu2wIAnsjcRLLunmQn9mO/N+omrjmoX7lyZVVV1bx58x599NHJkycLBDf8xehLeHj4gQMHuixkdkV7Kdi4cePGjRv70TAAeBndZQqI62a52J/3GHXN3mhBQcHRo0eDg4MXLVoUExPz0ksvnTlzxiUjAwAQEelKSXnD1U5ETp0bdROXnRtNTEx87bXXLly48Mknn+zbt2/cuHFJSUlvvPGGXs/zFgKAL9CVdXPRKPlWjFoslv379y9cuDA9PT0iIuKjjz7avHlzXl7ePffc46pVAMDg1dMtTGI/snW9d3yAuebc6NNPP/3ll1/Gx8c/9thjb7zxRlhYGLN8ypQparW6988CAPRNV07K6G6WixRk6s9T7VzINTEaHR2dl5cXE9N1l1ssFjM3GwAAcGIzkNivm+Ui/vdGXXNQf++993bJ0Pz8fOZFx54pAEA/2S09TpEnUpDNOLDddOWaGJ09e3aXJczl8QAALqAvI//o7n8k9iOrl58bra6uJiKbzca8YJSWlnrI5HV6vX716tUajSYhIWHJkiV8twMA/aLteYq8fh3U5+Xl5eTkEFFFRQXH1oh7jDK3q7e2tnbcty4UCsPCwtauXcu1NVdQKpWrVq3CPfUA3k1XSv4ju/+RuD8H9ampqcyN+S65p55rjDY3NxPRM8888+6773LvBgCgG7pLFPVo9z8SeflB/dWrVzUajVKpfOWVV26cti8iIoLL4AAA7XrZG/WAr5g4xei4ceP++c9/ZmRkjBs37saf1tfzfDEXAPgIUz3JQrr/kUhONtPAdtMVpxjtCEokJgC4i8NOvcx2JJKRvW0Au+mGay540mq1bW1tRGSz2Xbt2rVnzx5M/QkArmG4Qn69PBBUQHynjWvuYrrzzjs//PDDlJSUv/71r3v37pXJZCdPnvSQL+sBwLtpL3U/RZ7HcM3eaGlpaUpKit1u//DDD//9738fPXp069atLhkZAAY7bQkFxPPdRG9cE6Nisbi1tfX7778PCwuLjo5WKpUGA8+XIACAj9BeJNUovpvojWsO6h9//PHbbruttbX1+eefJ6KffvopOjraJSMDwGDXWkL+Hn1Q75oYXb9+/fHjx6VS6aRJk4hIKpXi2R4A4BrmRpIF891Eb1wTowKBYNSoUefPn9+3b59LBnSVtra23bt35+bmDh06dPr06Xy3AwB9a2pq6rgr0m5tI2F38+N15ejmoaE9Ky0tPX78OLnoYk3XxOiGDRuysrKSkpLkcnnHwi7Pg+aFSCSKjY2Njo5WqVR89wIArNTX13/11Vd/+tOfiOium5eS/3d9fEAkI7uZhDL2qwgMDExKSiKio0ePcui0nWti9O233z537tyQIUNcMpoLicXisWPHYmoSAO8SGRmZnp5ORHR1D7X29f2SUEY2k1MxGhwcHBwcTERKpbL/XXasn/sQRBQcHOyBGQoAXq/1HAWO7qNGJCcbnzcyuSZGFyxYkJmZefHixfpOXDIyAAxqrRdJldBHDbM3yh/XHNS/+uqrRPTBBx90XqjT6VwyOAAMXvqyHids7iCS83tbvWtiFIkJAG5ht5BA1EcN35M8uew59Y2NjZ9//vmWLVuIqLq6urKy0lUjA8AgZbhKiuF9l/lGjH799ddJSUk7duxYvXo1ERUXFz/99NMuGRkABq/W832fGCUiIc9z5bkmRjMzM7/88ssvvviCuXrg1ltvzcvLc8nIADB4Nf9M6rF9l4l4/orJNTFaWVnJPB+qA+YbBQCuWotJ1dfVTsT/V0yuidGJEyd2von+rbfemjp1qktGBoDBS8viaidiLnjy/m/qN2/enJ6evnHjxqqqqjFjxohEov3797tkZAAYrBxkt5BQ2nch3wf1ronRkSNHFhYW5ubm/vLLL8OHD580aZJY7JqROdLr9atXr9ZoNAkJCUuWLOG7HQBgTV9ByihWlUKnD+rz8vJycnKIqKKioh+tdeGysGtubg4LCxs1ahRzp6qHUCqVq1atwj31AN6n+SwFsvh+iYhEUrI4d+l6amoq83VOZmZmP1rrwgXnRj/77LOEhISQkJC4uLiQkJDExMQvvviC+7AAMKg1nyV1MqtKb7/g6eDBg0899dTChQsLCwvr6uoKCwsff/zxJ5988vDhwy7pDwAGqaYC0oxnVcn31CRcD+rXr1+flZW1fPly5m1ISEhycrJMJlu3bl1aWhrn9gBgsDJcIz8WtzARkVDq3XujeXl5s2fP7rJwzpw5+fn5HEcGgEFLLrKQxJ9tNd97o1xjtLm5OSIiosvCESNGNDY2chwZAAatiIAmUqewrRZKyW52Zzt94HpQb7fbjx07JhL9ZgoWm81mt9s5jgwAg1ZUQCPbE6PE/1dMXGM0KioqIyOj2+UcRwaAQSta1UBBrGNU5OV3MZWXl7uiDQCAX4X7tVAAi9tAGXzvjbpsvlEAAJcQ2E1Wu4gErNOJ771RxCgAeBaZ4dwvuiAnPsD33qhH3PnuPkaj8f333w8NDY2Kipo/fz7f7QBA3+S6M6UtIU58wPlv6ouKipjpk6qqqpz6YLd8PEZlMtl9990XFxcnl8v57gUAWJHrz1xuCXXiA84f1MfExDz22GNEVFpa6tQHu+XjMSoUCocMGXLjla0A4LEk5sp64zgnPuD8XUx+fn5+fn5EJJWymIivz/VzHwIAwGWMVTaJk7PECUTk4PNCdcQoAHiS+lNG5QS+m3AOYhQAPEn9CVPATXw34RzEKAB4koY8k98YvptwDmIUADyGVUcCkUMo47sP5/AQo7W1tTNnzlQoFElJSSdPnmRZ0O3CqVOnCq6bO3fuwG0DALhD/SkKncJ3E07j4YKnxYsXBwcHV1VV7dy5c86cOaWlpcyVB70X9PSpPXv2TJ8+nYi6zDIFAN6n9iiF3kbOPVeJfwO9N6rVavfu3ZuVlaVWqzMyMtRqdZfHjXRb0MunJBKJXC6Xy+USiWSAtwUAXKzuuDfujQ50jJaUlEgkktjYWOZtcnJyUVFRnwW9fGrx4sUxMTGPPPJIWVnZQG0EALiBpZVISGIl3304baAP6vV6vUql6nirUql0Ol2fBT196sUXX4yOjhYIBOvWrZs1a1ZhYeGN9yTYbDar1cqcPyWijhcA4EJ1dXVmc/uN7eHh4WJxH9lSX1/f1tZ+61FYWJhEIqGabyn8TqdW2tLScu3aNSIKtznEDisJ+g60jhnlmRcOh8OpNXZroGNUqVRqtdqOt62trXFxcX0W9PSpBx54gFmSnZ0dGhp69uzZm276zRVn9fX1c+fOlcvlGo0mJSWFiG699VbMUQLgcunp6cxk7adPn965cyfzz60Xs2fPHjZsGBEVFBRs27YtNTWVar6mqEfZrzEyMvKzzz57/vnnDQbD6rTLKWlm6iu7y8rK3nrrLSKqqakpLi4mIpfMtjHQMRofH282m0tLS0eOHElERUVFTzzxRJ8FfX5KJBKJRCKr1dpldSEhIV9++WViYqJ7twpg0BOLxTt27CCi5557juVHtm/fLhAI/vznP7e/b8ijCW+yX+PEiROZNV65cqVh12SytxH59f6RmJiYDRs2dF6SmZnJfo09GehzowEBAenp6VlZWTqdLicnp7GxkXkO8zvvvHPq1KmeCrpd2NLS8sUXX9TU1FRVVS1fvlyhUPT5BxAAPJSujPwiSdDP620sNiGPT7Xj4brRzZs3V1dXh4SErFmzZufOncx1Sx9//HFBQUEvBTcutNvt69evj4uLGz169NmzZ/fv369QKAZ+cwDABSr/TcNm9vvTFruQxwnwebhuNDw8/MCBA10WMruivRTcuFCj0Zw4ccJNTQLAgKrcT5Pe7/enbQ7R4NobBQDoTGg3kM1EiqH9HsFiE/L4HBHEKADwLNDwAw2dwWUEq32QnRsFAOhMo/uGImdzGYHfc6OIUQDgk1hglVkqKWAUl0FwbhQABq/R6ooW/6kcB7HYcW4UAAar8UEXGvyncxzEioN6ABikLFqVRN8m4frsXqudz4N6H3/Asl6vX7t2bVBQ0KhRoxYvXsx3OwBe46WXXqqoqGBev/nmm0OH9v9qpN78svNsU+xwIiI6efLkO++8Q0Q6nS4kJMSpYZw9qM/Pz//kk0+I6MqVK06tqFs+HqNKpXLFihUJCQlCIfa7AZxw7Nix7du3E9GKFStaW1vdFaPl235sSGRitLy8PDk5mZkuQ6l0bro8Zw/qb7rppnHjxhHRr3f0c+DjMUpEIpGozzm7AKALgUAQGRlJzieaE7SXSKrRWX+dT0StVjMrdZazB/UCgYCJBZdMm4l9NADgh6D0fYrNcMlQ+KYeAAYdqchONUdoyD0uGQ3f1APAoDM9vs4ROY8Erokgfr+pR4wCAA/uH11DsU+5ajTcDAoAg0z1VxfrlSQNctV4VpwbBYDBpXj9jjPDXTjeoJv9HgAGtfpTJAupbHXBs+Q6WB0iHNQDwKDx8z9o7CuuHRLTNgPAoFHzLcnDKCCu70pn8Ls36uO399jt9traWpVKJZPJgoOD+W4HYBBpaGgoLy//7TLH6MqVftO2d6ksLi4WCoVlZWUajaZ/66pv1DY1GEp//FEul48ZM6bPeqPR2NTURERmswvOqPp4jLa1te3evTs0NDQ6OnrBggV8twMwiOzZs+eTTz5JSkrqWDIh+GJR89WF6TGdy+64444jR47k5+czr/uxIpVKNSIm/tqVr7Ye3vrtt9+eOXOmz49cvnx53759RFRfX9+PNXbh4zGqUCiefvrpxMREvhsBGIwWLFiQkXH9dk+rgb6eds/aEQt/W/PQQw899NBDXNaiVqv/3z/X07HZby15a+pUVjNAjx07duzYsUSUmZnJZdUMnBsFgAFRsolGPWewiNwyuEiGb+oBwKc1F1L11xSzsO/K/hGIyGFz1+B98fGDegDgn91Muc/S5G1ELpiVzgNhbxQA3OzsaxS9gPxj+q70TohRAHCnmiPUeoFG/ZHvPtwIB/UA4C5KQRMVvETTDvrq4TwDMQoAbiESWO+QZNPET0mq5rsX98JBPQC4g2OK6IMLtjtIM57vTtwOMQoAbvDTiy2OoZetE/nuYyAgRgHAxcLq3ieH7YwtfUDXKpTyNeWoj58b1ev1q1ev1mg0CQkJS5Ys4bsdAJd59NFHq6urmde7du3q96QefYqKilq0aJFUKiWil19++e677+7807vuusvhcBDRhQsXmCUzIk5/d7BsS35SVdV/XnjhBTd11Q2RnGwmlrV5eXk5OTlEVFFRwX3NPh6jSqVy1apVuKcefM8vv/xy5MgRIlqwYEFbmxvvg/yv//qvv/zlL0T07rvvdgR3h7a2tm+//ZbaH/juoDMvPzA91Truu9kkICLmWfADRCQnm5FlbWpqampqKrnonnofj1EAHyaRSKg9v9xIIBAwKxIKuzkH2PFTslvo1FPkHydIXS1xa0M9ESnY7426Fs6NAgBnxio6MpOG3ENjX+atB2cO6l0Le6MAwE3tUTr9PKVuoeBUPttw5qDetRCjANBf9rZnJlXQhbdo2kGShfDcDA7qAcDLNJ6mb2bU6aV02y7+M5T43BtFjAKAk8yNdPrPdGYV3frRlz8P5bub60R+ZNXzsmbEKACwJSQLXXyHjsyi8Ltp2n9IGcV3R52IlWQz8LNmXtYKAF7GZhorPTqa8omW0fTvSCjju6EbiJXYGwUAj9RWR+fW0ldTRGQ5TC/TqD96YoYSkdiPrPzsjSJGAaAHDT/Q90/S0QdJFkb3nDpjTrMRP1fWsyLibW/Uxw/qrVbruXPnDAaDSqWKi4vjux13KSoqKi8vZ16npKRERkZyGe2rr74ym9uneJg5c6ZI5J5HObrN8ePHW1pamNd33nmnUql0doSrV692POs8KiqKeRLvgQMHbLb2h6bdd999fQ5y6NAhi8XCvJ41a1a3twB1Kz8/v6amhnk9adKkkJD+fwlutVoPHjzIvJbJZGlpaZ1/Wl5eXlRU1OUjISEhk8YOpSs73n34DF3cRPHPUsgtHT8tKChQqzlNHtrS0rJ//34iampq4jJON8R+7M+NNjY2Mv9k9HoXJK+Px6jNZrtw4UJra+uwYcN8OEa3bNkiEAiGDh36008/paWlPfPMM1xGy8zMXLhwIRF9+OGH/Yshfi1fvpx57vmuXbtycnL6MaPC4cOH9+7de/PNN9fW1hoMhn/9619EtHTp0kWLFhHRu+++yyYWMzMzn3zySSL64IMPpk2b5ufnx3Ltr7322oQJE+Ry+TfffLNixYp77rnH2f47GAyG5cuXM23k5OScO3eu80//93//Ny8vj/kjQUQacUOcvEh+4QcyTqaoRzJ3jz385486199yyy3Nzc2FhYWdF/7hD39wqqWFCxcyIzz++OMuvo1VrGR/UN/U1MS0YTC44DyAj8eoTCabNWvWYJia5Kmnnho/fvzWrVtNJq5XIAcHBzNTUTBTTngdiUTC9H/58uV+DzJr1qzf//73RUVFb7/9NrMkPDycGfY///kPmxFCQkKYemYCEacsW7ZMo9F07MxyMXLkSKYNZh+wi3mz731oahhV/ofqfyC/4ZYhj85datj9/D4iMlnf7VI8fvz48eO5zsG8YsUKjiP0SKwkq45lbWxsbGxsLBGdPn3aBWvmPgQAeBNLK9V/Tw15D6u3yoQMjmMAAA1fSURBVCwyqnmQIh+m8f9NArG9rU1n/qjvETyTRE2WFl7WjBgF8HVWA7WcvW90LeUtIXMTGa9RyC0UdueeFkn48Ni5v5vLd38uIlGRGTEKANw57KQvo5Zz1HJOUZv/3/f8QEfvp4BRAnJQ/B8oMIkE7f/qLY4L/HbqYpIAsrbysmbEKIDXshrIVHPz0CrlL+/RtXrSV5DhFxLKKDCJAuJJPbYt+L6Vr+v3L/s3Ee17fupL6hS+O3YngYgcdl7WjBgF8Gw2E5lqyVhJusvU1kCGX8jwCxkr337wZzqaTsqoWHWTXayiETPJP6bLFCH21laHg6++BxHEKABfHGRuorZGMtWQzUimajLVk6lm8c0/+59+nGxNT0VUBjRvp2MjSR5KimGkGELKERQ6hfwiSR6W+bepx1d8TUQ73ph729D5gcFD+N4cTyAg4uHvhnfHaG1t7cKFC48ePRoTE5OdnT158mS+O4LBytxEVh3ZLWRuJpuBbKY7Yxvo0ntkbvy/Ey5Lf15JZCZzE5mbyW4maRAZKkjsT1INSYNIFkQSFclCyS+CQib+u+T7qYs3BIbFfvj31bfccss9d/b/utFBR6r2l9UP/Gq9O0YXL14cHBxcVVW1c+fOOXPmlJaWsr/IGQYjh739mhirgexmsreRVf/rQksLOexkt4wSfktiM51rCKsqnTMylwpeIoc1654LdPwRMjdtfKhQePg2EghJICBbG4kUJBCRQEhiJQmlJAkkiT/Jh2gUFhIrSJ7w47WgORGPilVDSKohSSAJ+rhu/0pLgF0aSgIvu3nMIyiGBikqB361XhyjWq127969xcXFarU6IyPj9ddfP3z48AMPPMB3X4OJRUsO669v7Zau1z8zwdTOQebm3/zU3ExCSfuzxW2m38y5a25uPzrryDgicljJor2+rrb2W1Y6LyQHmVs2PHCJvplOdvPz4y6NuFhI5XYiIrE/CSVERAIRSVRERPJhZDcSEUk1RESSABKImZhzCMQmh5zUKfrW4PNN12ZEziGhbMupb+96fi1J1S/8I/3AsmN93sX05c97l0f/HyI6W7PJrh5H+AM/AJRRwwNPDPxqvThGS0pKJBIJcysCESUnJxcVFXGNUYuWSjY5Ud/xj5klc8/3EcvDyFR7vY2WPr5ztJs7z8Lwx6QLw67kUVPAdFGVXWGnbz7v2phUTXT9xrvO2SeUkNif/KNJV97xkffnXaADN5NUvfaeAvnJdBJY24cSiEkS0LUTkZxEih77/M3Ol4CkaiIiierXXS2xP5GDhNLrJdczjok2aadnr3d8SiglsfKGhRIS+zPLlv1j6vEXvyKiN7f//oW0F/pxD1vJcanVZqVhs/RNReeb8ih4IhFVtsrJfyQRmazYT/RU6t8lhnW9+WoAeHGM6vV6lUrV8ValUul0XW8Fq6+vnz17tkwmU6vVKSkpRDR58uQFCxb0OKjhCgml5Bfxm4UdOzLsdf733zux3/Vpx3o5O349gHpwqfVAQOodAcOGfbdjx9///ne1mtmtk3cq6XyHqPi3/92NROevv5YTkUAw7tixY0S0ZvO8yl2VAoGg01A3PqTBSNTLHBM8HGHV1rb/NfL391+4cKFUKu29/kYtLS0rV64kIrFY/PXXX0+dOpWIZLL22eECAwNvv/32Pgfp2F1VqVTTp09nf//41atXmc/6+fmtXLkyKyur2zKjsf2/RUBAwOzZs7udQcZut48YMYJ57XA4mA3p0NTUtHbt2s5LBAJBSUkJU1Zfz8NJRq6G3D1iTFovPy8tLd2wYQMR1dTUXLhwgYgUip53AlgTOLz2gojTp0/fdtttHRO0zJs3b8KECcztwx02bdp01113DYZ76gGgHzIzMzumTeg3L55vND4+3mw2l5aWMm+LiorGjBnDb0sAMAh5cYwGBASkp6dnZWXpdLqcnJzGxsYu0ykCAAwAL45RItq8eXN1dXVISMiaNWt27tzJ/WonZn5Jl/TGi5aWFtfPhjuAjEZjdXU13130n91uv3LlCt9dcFJWVsZ3C5zw0r93x2h4ePiBAwdMJtP58+enTJnCfcB9+/ZdvHiR+zh8yc3NPX78ON9d9F9FRcWuXbv47qL/jEbj5s2b+e6Ck/Xr1/PdAie89O/dMQoAwDsfj9Ha2tobr4Lqvb611Ym5tsrLyxsaGtjXX7lypeNyHDauXbvm1EFuY2OjU/3U1tY6dRDa0NDQ8dAnNpqbm52agl6n0zn1+9Hr9cXFxezrTSbTjU8f6oXZbO7yzIze2Wy2jscoseFwOH788Uf29USUn5/v1nqnfv/9GN/T+ne2vls+HqPFxcVO/ZouXrxYVVXFvv7EiROXLl1iX5+bm+vUP/vTp0+fPXuWfX15eblTsVVcXJybm8u+/tKlSydOOHGXSEVFhVNPIqmtrXXq91NbW8vykR6MpqamPXv2sK/X6XROnWRgzi+xr7fb7du2bWNfT0Rbt251a71Tv/9+jO9p/Ttb3y0fj1EAAHfz4svv2di0adOpU6fYPxL24sWLQUFB7B9pW19fL5fL/f39WdY3NjZKJJKAgBtuqexBU1OTUCgMDAxkWX/16lWbzRYVFcWyXqvVWiyWoKAglvU6nc5kMrH//RgMBr1eHxoayrK+ubm5qqpq9OjRLOvb2tqampqGDGE7R5zFYqmrqxs2bBjLeqvVWl1dHRER0Xfp9frCwsIJEyawrCeiiooK9v+9iKi8vDw6Otp99bm5uRMnTvScftzdf2Njo7MHBDfy8RgFAHA3HNQDAHCCGAUA4AQxCgDACWIUAIATxCgAACeIUQAAThCj7Wpra2fOnKlQKJKSkk6ePMl3O73R6/XPPPNMZGSkv7//tGnTOu5u7HYTPHm7vv76a4FA8MYbbzBvvav/devWRUZGKhSKSZMmWSwW8qr+z5w5M3nyZKVSGRcX9/HHHzMLPbz/tWvXjhkzRigUZmdndyxk37N7N8QBDofD4Zg9e/Zjjz3W1NSUnZ0dHh6u1+v57qhHNTU1L7zwwpkzZ6qqqv74xz+OGjWKWd7tJnjsdlkslgkTJkyYMOH1119nlnhR/++///6oUaOOHTtWU1Nz8OBBq9Xq8Kr+x4wZ8+qrr5pMpm+++UahUJSVlTk8vv/t27cfOHBg8uTJ7733XsdC9j27dUMQow6Hw9Ha2ioWiy9dusS8TUhI2L17N78tscTMGNDc3NztJnjydm3YsGHlypUPP/wwE6Pe1X9MTMyhQ4c6L/Gi/u12u0gkunz5MvM2MTHx0KFD3tL/3Xff3RGj7Ht294bgoJ6oh4eM8tsSS7m5uSNGjAgMDOx2Ezx2u+rq6jZv3rxq1aqOJV7Uf319/ZUrV/Lz80NCQkaOHPmvf/2LvKp/gUAwY8aMrVu3ms3mI0eOtLS03HzzzV7Ufwf2Pbt7Q7z4yaAuxOYhox6orq5u6dKl69atox42wWO3669//esLL7zQeW4BL+q/srLSZrPl5+eXlpYWFxenpaWlpKRYrVZv6Z+INmzYcNddd7322msSieSjjz7SaDQ///yzF/XPYP//jLs3BHujRERKpVKr1Xa8bW1tZT/bCF+0Wu2sWbOeeeaZuXPnUg+b4Jnbdfr06dOnTy9atKjzQi/qn3kk78qVK1Uq1cSJEx988MGDBw96Uf8mkyktLe21115ra2v74Ycfnn/++R9++MGL+u/Avmd3bwhilMgLHzJqNBrT09PvvPPOl156iVnS7SZ45nZ99913Fy5cGDZs2JAhQ/bt25eVlbV48WIv6j86OtrPz6/jufPMCy/qv6SkpLa2dtGiRVKpdPz48VOmTPnmm2+8qP8O7Ht2+4a48DyrV5s9e/aTTz6p1Wq3bt3qOd+odstisdx3330LFiwwGAxGo9FoNNrtdkcPm+CB26XX66uuu//++1955ZXm5uaeWvXA/h0OR0ZGxrx583Q63Y8//hgYGPj99987vKf/1tZWpVL50UcfMdP6hYaG7t27t6dWPad/i8ViNBqnTZu2efPmfvw/79YNQYy2q66unjFjhkwmS0xMPH78ON/t9ObG+dXPnz/v6GETPHy7Or6pd3hV/y0tLXPnzvXz84uOjs7OzmYWelH/+/fvT0lJ8fPzi4iIyMrKYhZ6eP9Llizp/P/8sWPHnOrZrRuC+UYBADjBuVEAAE4QowAAnCBGAQA4QYwCAHCCGAUA4AQxCgDACWIUAIATxCgAACeIURgUzp07FxISUlZWRkRHjhyJiIhobm7muynwEYhRGBSSkpJWrFiRkZGh1WozMjK2bNmiVqv5bgp8BG4GhcHCZrNNnjxZp9PddNNNW7du5bsd8B3YG4XBQiQSPfnkk+fOnfvTn/7Edy/gU7A3CoNFbW1tSkrK9OnTr127xjyUlO+OwEdgbxQGi+eee+6JJ574n//5H4PBsGXLFr7bAd+BvVEYFHbs2PHqq68WFBTI5fLz58/ffvvt+fn5UVFRfPcFvgAxCgDACQ7qAQA4QYwCAHCCGAUA4AQxCgDACWIUAIATxCgAACeIUQAAThCjAACcIEYBADhBjAIAcIIYBQDgBDEKAMAJYhQAgBPEKAAAJ/8fvNkD02YXEXQAAAAASUVORK5CYII=",
"prompt_number": 20,
"text": [
"FramedPlot(...)"
]
}
],
"prompt_number": 20
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment