Skip to content

Instantly share code, notes, and snippets.

@abikoushi
Created December 26, 2018 15:07
Show Gist options
  • Save abikoushi/faeb072a817d0e70615c65b757ff2907 to your computer and use it in GitHub Desktop.
Save abikoushi/faeb072a817d0e70615c65b757ff2907 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import tensorflow as tf\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"%matplotlib inline\n",
"\n",
"def tf_diff_axis_0(a):\n",
" return a[1:]-a[:-1]"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"data = pd.read_csv('MERS2015.csv', sep=',')\n",
"d_Y = np.diff(data[\"Recovered\"])\n",
"N=len(d_Y)"
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"tf.reset_default_graph()\n",
"\n",
"Y = tf.placeholder(tf.float64, shape=[None])\n",
"S_raw = tf.Variable(5, dtype=tf.float64)\n",
"a_raw = tf.Variable(-1, dtype=tf.float64)\n",
"b_raw = tf.Variable(-1, dtype=tf.float64)\n",
"\n",
"S0 = tf.exp(S_raw)\n",
"a = tf.exp(a_raw)\n",
"b = tf.exp(b_raw)\n",
"\n",
"def SIR(state, t):\n",
" S, I, R = tf.unstack(state) \n",
" dS = -a * S * I \n",
" dI = a * S * I - b * I\n",
" dR = b * I\n",
" return tf.stack([dS, dI, dR])\n",
"\n",
"t = np.linspace(0, N, N+1)\n",
"tensor_state = tf.contrib.integrate.odeint(SIR, [S0,1,0], t)\n",
"diff_R = tf_diff_axis_0(tensor_state[:,2])\n",
"\n",
"cost_func = - tf.reduce_sum(Y * tf.log(diff_R) - diff_R)"
]
},
{
"cell_type": "code",
"execution_count": 94,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"EPOCH: 200, cost: -150.224\n",
"EPOCH: 400, cost: -151.046\n",
"EPOCH: 600, cost: -151.057\n",
"EPOCH: 800, cost: -151.057\n",
"EPOCH: 1000, cost: -151.057\n",
"EPOCH: 1200, cost: -151.057\n",
"EPOCH: 1400, cost: -151.057\n",
"EPOCH: 1600, cost: -151.018\n",
"EPOCH: 1800, cost: -151.057\n",
"EPOCH: 2000, cost: -151.057\n"
]
}
],
"source": [
"n_epochs = 2000\n",
"n_epochs_report = n_epochs // 10\n",
"\n",
"train = tf.train.AdamOptimizer(0.1).minimize(cost_func)\n",
"sess = tf.Session()\n",
"sess.run(tf.global_variables_initializer())\n",
"\n",
"for epoch in range(n_epochs):\n",
" cost_val,_ = sess.run([cost_func,train], feed_dict={Y: d_Y})\n",
" if (epoch+1) % n_epochs_report == 0:\n",
" print('EPOCH: %d, cost: %.3f' % (epoch+1, cost_val))\n",
"\n",
" \n",
"state,dR, loga, logb,logS0 = sess.run([tensor_state,diff_R,a_raw,b_raw,S_raw])\n",
"sess.close()"
]
},
{
"cell_type": "code",
"execution_count": 96,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.PathCollection at 0x15c5751d0>"
]
},
"execution_count": 96,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAFkCAYAAABxWwLDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xd4VVXWx/Hvzk0vBBJICC0UIyBNCVUEHVQEARFUxmDF\n7gg6iG3UsY29ALZRxzaOBUTAFxsWFCyhaeIIMgihJRQJcAMhjdT9/nGSkMRQQspN+X2e5zy5d59z\nz1nXQha7rG2stYiIiIhUhZenAxAREZGGRwmEiIiIVJkSCBEREakyJRAiIiJSZUogREREpMqUQIiI\niEiVKYEQERGRKlMCISIiIlWmBEJERESqTAmEiIiIVNlxJxDGmBuMMb8YY9KLj2XGmJEVrnnIGLPT\nGJNtjPnKGHNChfN+xpgXjTF7jTEZxph5xpiI441JRERE6kZ1eiC2AXcCfYFY4BtgoTGmO4Ax5k5g\nCnAdMADIAr4wxviWuccsYDRwATAMaAPMr0ZMIiIiUgdMTW6mZYxxA7dZa980xuwEnrLWziw+1wxI\nBa6w1s4tfr8HuNha+2HxNV2BdcAga+2qGgtMREREalSNzIEwxngZYy4GAoFlxphOQGvg65JrrLUH\ngJXA4OKmfoB3hWvWAyllrhEREZF6yLs6HzbG9ASWA/5ABjDeWrveGDMYsDg9DmWl4iQWAJFAXnFi\ncbhrKntmOHAOsBU4WJ34RUREmhh/oCPwhbXWXZ0bVSuBAH4D+gChwIXAf4wxw6p5z6M5B3i3lp8h\nIiLSmF0CvFedG1QrgbDWFgCbi9/+bIwZANwCPAkYnF6Gsr0QkcDPxa93Ab7GmGYVeiEii88dzlaA\nd955h+7du1cn/Hpj2rRpzJw509Nh1IjG9F1A36c+a0zfBfR96rPG9F3WrVvHpZdeCsW/S6ujuj0Q\nFXkBftbaLcaYXcCZwGoonUQ5EHix+NoEoKD4mrKTKDvgDIsczkGA7t2707dv3xoO3zNCQ0P1Xeop\nfZ/6qzF9F9D3qc8a03cpo9pTAI47gTDGPAoswpn0GILTHXI6MKL4klnAvcaYjTiZzj+A7cBCcCZV\nGmNeB2YYY/bhzKF4DojXCgwREZH6rTo9EBHAW0AUkI7T0zDCWvsNgLX2SWNMIPAK0Bz4Hhhlrc0r\nc49pQCEwD/ADPgduqkZMIiIiUgeOO4Gw1l5zDNc8ADxwhPO5wNTiQ0RERBoI7YVRD8TFxXk6hBrT\nmL4L6PvUZ43pu4C+T33WmL5LTarRSpR1wRjTF0hISEhojJNaREREak1iYiKxsbEAsdbaxOrcSz0Q\nIiIiUmVKIERERKTKlECIiIhIlSmBEBERkSpTAiEiIiJVpgRCREREqkwJhIiIiFSZEggRERGpMiUQ\nIiIiUmVKIERERKTKlECIiIhIlSmBEBERkSpTAiEiIiJVpgRCREREqkwJhIiIiFSZEggRERGpMiUQ\nIiIiUmVKIERERKTKlECIiIhIlSmBEBERkSpTAiEiIiJVpgRCREREqkwJhIiIiFSZEggRERGpMiUQ\nIiIiUmVKIERERKTKlECIiIhIlSmBEBERkSpTAiEiIiJVpgRCREREqkwJhIiIiFTZcScQxpi/GWNW\nGWMOGGNSjTEfGmNOrHDNm8aYogrHZxWu8TPGvGiM2WuMyTDGzDPGRBxvXCIiIlL7qtMDMRR4HhgI\nnAX4AF8aYwIqXLcIiARaFx9xFc7PAkYDFwDDgDbA/GrEJSIiIrXM+3g/aK09t+x7Y8yVwG4gFvih\nzKlca+2eyu5hjGkGXAVcbK39trhtMrDOGDPAWrvqeOMTERGR2lOTcyCaAxZIq9B+RvEQx2/GmH8a\nY8LKnIvFSWK+Lmmw1q4HUoDBNRibiIiI1KDj7oEoyxhjcIYifrDW/q/MqUU4wxFbgC7AY8BnxpjB\n1lqLM6SRZ609UOGWqcXnREREpB6qkQQC+CdwEjCkbKO1dm6Zt2uNMWuATcAZwJIaeraIiIjUsWon\nEMaYF4BzgaHW2t+PdK21dosxZi9wAk4CsQvwNcY0q9ALEVl87rCmTZtGaGhouba4uDji4irO0RQR\nEWl6Zs+ezezZs8u1paen19j9jTOScJwfdpKHccDp1trNx3B9OyAZGGet/aR4EuUenEmUHxZf0xVY\nBwyqbBKlMaYvkJCQkEDfvn2PO3YREZGmJjExkdjYWIBYa21ide513D0Qxph/4izJPA/IMsZEFp9K\nt9YeNMYEAffjzIHYhdPr8ASwAfgCwFp7wBjzOjDDGLMPyACeA+K1AkNERKT+qs4Qxg04qy6WVmif\nDPwHKAR6A5fjrNDYiZM43GetzS9z/bTia+cBfsDnwE3ViEtERERqWXXqQBxxCai19iAw8hjukwtM\nLT5ERESkAdBeGCIiIlJlSiA8qKAAPvwQfvoJUlOhqMjTEYmIiBybmqoDIcdh506YMOHQe19faNcO\n2reHLl3gxBOha1fnZ5cu4OfnuVhFRETKUgLhQe3bOz0P27aVP1JSYPVq+OADyMhwrvXygpgY6Nv3\n0HHKKdCihWe/g4iINE1KIDzIGIiIcA5nWW551sLu3bB+vXOsWQOJifDRR5CV5VwTEwN/+tOhIzLy\nj/cRERGpaUog6jFjnIQgMhKGDTvUXlgISUmQkADLlsGSJfCvfznnuneHs86CCy6A004Dl8szsYuI\nSOOmBKIBcrmgWzfnuOQSp+3332HpUieZ+PBDeP55aN3amWNx0UUwdKiSCRERqTlahdFIREVBXJzT\nE5Gc7PRMxMXBxx87Qxvt2sH99zsTN0VERKpLCUQj5OUFgwfDjBmwdSssX+70RDzzDERHO4nFsmXO\nHAsREZHjoQSikfPygkGD4MUXYccOePpp+PFHGDIE+veHTz5RIiEiIlWnBKIJCQ2FW26BDRvg008h\nOBjGjnWGOH780dPRiYhIQ6IEogny8oJzz3UmXH76KezdCwMGwMUXw6ZNno5OREQaAiUQTZgxTiLx\nyy/w+uvw/ffOMtD77oP8/KN/XkREmi4lEILLBVdd5dSWuOsuePRRZxLmunWejkxEROorJRBSKjAQ\nHnoIVqxwKl327QvPPqtNvkRE5I+UQMgf9OvnlMy+7jr4619hxAjYvt3TUYmISH2iBEIqFRDg9D58\n9ZWzD0e/fs624yIiIqAEQo7irLOcPTc6dXL241i40NMRiYhIfaAEQo4qIgK++QZGj4bx42HWLBWf\nEhFp6pRAyDEJCID334fbb4dp0+Dmm6GgwNNRiYiIp2g3TjlmXl7wxBPQpQv85S/OxMq5c8HHx9OR\niYhIXVMPhFTZddc5cyE+/RQuvxwKCz0dkYiI1DUlEHJcRo+GOXOcHogbbtCcCBGRpkYJhBy3CRPg\nzTfhtdfg1luVRIiINCWaAyHVcvnlkJEBU6Y4u30+8ICnIxIRkbqgBEKq7aabnCTib3+DkBCYPt3T\nEYmISG1TAiE14q67ID0dbrsNunaFMWM8HZGIiNQmzYGQGvPIIzBuHFx6qbOzp4iINF5KIKTGeHnB\nW29BZKQzwTIz09MRiYhIbVECITUqNBQ+/BC2bIFrrtHKDBGRxkoJhNS4k05ylne+/z7MnOnpaERE\npDYogZBacdFFzr4Zd9wBS5d6OhoREalpSiCk1jz6KJx+OkycCLt2eToaERGpSUogpNZ4ezvlro2B\n66/XfAgRkcbkuBMIY8zfjDGrjDEHjDGpxpgPjTEnVnLdQ8aYncaYbGPMV8aYEyqc9zPGvGiM2WuM\nyTDGzDPGRBxvXFK/tGoFL78MH30E77zj6WhERKSmVKcHYijwPDAQOAvwAb40xgSUXGCMuROYAlwH\nDACygC+MMb5l7jMLGA1cAAwD2gDzqxGX1DPjx8Mll8DNN8OOHZ6ORkREasJxJxDW2nOttW9ba9dZ\na9cAVwIdgNgyl90C/MNa+4m19lfgcpwE4XwAY0wz4CpgmrX2W2vtz8BkYIgxZsDxxib1z3PPQUAA\nXHuthjJERBqDmpwD0RywQBqAMaYT0Br4uuQCa+0BYCUwuLipH0457bLXrAdSylwjjUBYGLz6Kixa\n5CzxFBGRhq1GEghjjMEZivjBWvu/4ubWOAlFaoXLU4vPAUQCecWJxeGukUZi9GiYPBn++ldISfF0\nNCIiUh01tZnWP4GTgCE1dL+jmjZtGqGhoeXa4uLiiIuLq6sQ5DjMnAlffQVXXw1ffums0BARkZo3\ne/ZsZs+eXa4tPT29xu5vbDUHpI0xLwBjgaHW2pQy7Z2ATcDJ1trVZdqXAj9ba6cZY/4ELAZalO2F\nMMZsBWZaa5+t5Hl9gYSEhAT69u1brdjFM778Es45B954w+mREBGRupGYmEhsbCxArLU2sTr3qtYQ\nRnHyMA74U9nkAcBauwXYBZxZ5vpmOKs2lhU3JQAFFa7pijMZc3l1YpP6a8QImDTp0BbgIiLS8FSn\nDsQ/gUuASUCWMSay+PAvc9ks4F5jzFhjTC/gP8B2YCGUTqp8HZhhjDnDGBMLvAHEW2tXHW9sUv89\n9RRkZ8MDD3g6EhEROR7V6YG4AWgGLAV2ljkmllxgrX0Sp1bEKzirLwKAUdbavDL3mQZ8Aswrc68L\nqhGXNABt2sDf/w7PPw+//urpaEREpKqqPQeirmkOROORlwe9ekHbtvD115pQKSJS2+rNHAiR6vD1\ndQpMLVkCH3zg6WhERKQqamoZp8hxOeccOP98mD7dqRMRFOTpiESkXrL2UBnbij+r0lbxXMXXXl7g\n7//H6+UPlECIx82YASed5Gz//cgjno5G5BC3201aWhphYWGEh4cfe5u1kJtL2o4d7N+1ixYBAbQI\nDISDB0lPTSVj716a+fvTzM8P8vLI2LuXrH37CPbzI9jXF/Lzydq3j5wDBwj09SXQxwfy88nJyOBg\nZiYB3t74e3tDQQEHs7LIy87Gz+XCz+WCwkIoLCQvJ4f83Fx8vLzw9fKCoiLyc3MpyMvD2xh8itsK\n8vIozMvD5eWFtzFQVERhfj6FBQW4jMFV0lZQQFFhIS7Aq7itqKiIooICvIxxurOtLW23RUV4Aaa4\n3RYfBjDFyYAt8wu+tK3sew842KsXWUuWlP67lcNTAiEe16kT3HknPPaYUxfihBOO/hmRGlVURNqW\nLaQnJ9PCGALy8ljx2Wek/fILfpmZ+BlDq5Yt8cnNJXf7dnxzctgH5AcE4J2bC+npROTn41NYSIG1\nuPLzMdYSBoRVeFRo8VFWSPEBUOhyUeRy4QUEGEOBy0Wmvz9FLhe5+flYazno7Y13cDDW5SInKwtr\nLUXe3viHhdGydWt2791Lxr592KIiinx8CAwPx3p5kZGWhi0owPr6EhQRgfXy4sCePc4vdy8vQlq3\nxnp5kZ6aCnl5FPn50axNGzCG/b//js3Px/r6EtquHRjDvh07IC8P6+dHeHQ0XXv04Ld169ibnIzN\nzQU/P5p36ABeXqSlpEBuLtbPj7DoaKyXF2llrgvr2BEAd3IylLR16uS0bd1a+tnwTp3AGPZu2eJc\n5+9P5Akn0Ovkk1n9yy+kJiWVfj68c2fn2s2b4eBB8PevtK1lSdumTeQUFbHxb38jZsgQJkycSEBA\nAFI5TaKUeiEnB7p1g/79Yd48T0cjDV5REezZA7//7hy7dzvHnj2wezd5O3ZQuGcPvhkZeKWnY9PS\n8DrMn4XW5aLQxwd3URH5LhctQkPxDQggG/glPR3j5UW3tm0JCQpiX2Eha7Oz2R0aSrOcHHqFhdGy\nWTNSc3P5z4YNBBjD+N69iWrenO3Z2Ty7ejUhxnBF//60b9GC5IwMHvvxR8KA6wcMIDo0lOT0dB5b\nteqY2halppISHk4Ht5tRkZFV/nx12hrbsxelptJu3DguueKK2v1vtY7V5CRK9UBIvRAQAP/4B1xx\nBaxcCQMHejoiqbesBbcbtmxxNlUpe2zbBjt3QmoqFBSU+1hRSAhEROA2Bnd+PtkuFznBwewNC8O3\nTRtOat2a8LAwNmRl8WpiIud27MjYP/0J4+1NenY2D3/8MScDcWPH4hMYSEF2Nh8Ut8WOGIF/YCBR\nQFJyMrOXL+feQYOILv5bdWh2Nu6NGzkZaNexI4GBgURkZ+OTkEA3oHNYGIGBgbQ1hpbZ2ZwMnBAS\nQqC/P22Lio6prbe/P/sPHuSH77/n8kGD6B0ZCXDMn69OW2N7du/iORAL4uNxjxmj4YzDUAIh9cYl\nl8DTTzsVKr/5Rss6m6Jycwl8feG338hctYq8NWsI2rULv23bYNOm8iVMAwOhQwfyoqLIiY7G5/TT\nCezShdywML5Zu5Z1SUlkFRZSGBrKHmtpm5rKmLZt6RIayi+pqby5ZAmXdupE9KmnAhDiduP722+Y\nrCyy8/II9PYmLSeHEGOIBHJycggMDKy0DSDYxwffnBwifH1LQzzWz1enTc+uuWcDRIeGQkoKaWlp\nSiAOQwmE1BsulzMPYswY+OILGDnS0xFJbSpNFkJDCdqxg4RXXsEsW0bEnj24MjMhMxOAYCA9MJCd\noaHkdelCxplnsjk9nX3e3uyPiCDitNPAy4vkVaucz6SnExMQQN6BA+xes6a0S/uX1FRmLFnCnzp1\nKv1balRwMNE+Pvjv20d2djaBgYGEBQQQFBhIalZW6S+UsIAAMqwlFUrHxCtrA8jMzycvIIDdeXl0\nKm471s9Xp03PrrlnAySnp0NwMGFhFWexSAklEFKvnHsuDBvm9EKMGOGsqJKGrzRZaNGCoL17SZg1\nC+/vv6fN7t2E7N+Pb0EBQ4C8kBBMZCT7oqJ4MS2NjIAARg8ZQoeWLckoGaveurX8+PUrrxwav+7Q\ngeT0dObNncumrCzu7N69yslCeGAgIeHhfH/gAD0zMvAOCmJHRgZ7AwNJBPpnZBDt5VVpW3J6OsvT\n0+kwdChL3W4CUlOJDg095s9Xp03Prrlnl8yBiBk3Tr0PR6BJlFLvrFgBgwfDO+84wxrScFRczpiT\nnc03M2bg+vRTOuzYQXu3m5DsbAByQ0MxbduyLTCQf6ek0C8mhnFnneXcJzubB8rMOQgMDDzmNoCP\n16/n36tW8cKYMUS1aFF6z2e++IIu2dmcP3Jk6S+Gx777jvVbtnDj6afTrXVrktPT+WjnTn6PiKCV\nMU6vRnAw0f37gzGHejoO0xYzZAijxo5l0ccfkxQfX+XPV6dNz665ZzfWVRg1OYlSCYTUSxMmwM8/\nw2+/gZ+fp6ORo8nJyWHB3LkkxcfTbPduOu/fT5/cXFqvXUtARgZFLhdFUVHsbN6ct1NT6RUTw3ln\nnw1AktvNC59/ztlBQQwfMYLAwECS3G7e/PprTgUGn3km4eHhx9wGsGXfPm7+5BNuHzCAYV27lsZ5\nrMlCyS+P7Ozs46sDUey460hUs03Prtl7NiZahSGN3iOPQM+e8MorcPPNno5GjshavnnmGTq89x4T\n9uwhYO9eLLAzNJSvmjen5+mn07lPH7y8vcl1u9n9+ed4ZWXVyJyDw41fZ+Tl4de+PUsPHKB5cZd2\ncno6hSEhBJ11FouMYVFKipMsTJjA9MMkCwEBAX/4RRIeHn5MbVW5tqbb9OyavadUTgmE1EvduztF\npf7xD7jySmjWzNMRSTnWkr54MUXvvUfw4sWM3r6dQl9fXCeeCKefjuncmcRt23hr1Spe6NwZvJ0/\namp6zsGRxq9HX3IJvr6+LIiPd5Z4HkeyICKHpwRC6q0HHoB334VnnoEHH/R0NAJASgr5b75J9ssv\nE7prFxkBAaxo04YVXbsyeeRIWjZvXnppz4gI8oCk3btL5yFUliwcrmdgxA03gDEsWLWqNAE41raY\nceNKx6/dY8YoWRCpBUogpN5q1w6mTIGZM+GWW0CrqTwkL88pD/raa7B0KcbHh22RkU6Bne7dsXv2\nEL9kCZ1++YULTz+99GM1NYzgHj/+uNtA3dIitUWTKKVe273b2Stj2jR4+GFPR9PEpKY6k1Beegl2\n7SJ74EAOjB/PG+vWMSYgoHR5JFQ+OXFRaioRo0bh6+tbbmb8kSYoikjt0iRKaTIiIpxeiGefdZII\n/Z6pA4mJMGsWzJmD9fEhaeBAvjrtNPb4+OBOTCQjOZk2Q4aU+8jkU05hSnY2c3JyCNUwgkiToARC\n6r3bb4cXX3TmQjz6qKejacQSEpzJJh9/TH67duTdey+ftWzJ5sWLGdWy5aFqjrt2sTQxsdxwxe7s\nbHr27s0lt98OoGEEkSZACYTUey1bwtSp8PzzcOutznupQWUShwORkSw+4wx+jYwkc+NGdnzzDbfG\nxJQOVwyLjia+c2c+2bKF9tHR5YYrYsaNIyYmxsNfRkTqigoFS4Nw223Oz6ef9mwcjUpSEowfD/36\nwW+/EX/ddbw8aBAndOvGLR07MjQ/n/yNGzm4c2e5j00+5RQyIyOZk5PDzJQUFhQW0q54uEJEmg71\nQEiDEB7uFJR69lmYPh1atfJ0RA1PaZU9Hx/CX3oJZs6kMDKSPU8+ScaYMXw1cyYToqJKextObd+e\nL0NC+G3zZk7p0aO0TPTRhitEpGlQAiENxvTpzjDGU0/Bk096Opr66w/7URSXmd74ww/0XrOGs1av\npqCggNWjRvFZeDgFP/+M+6ef/jA5MjwwkN6dO/P96tXEbNtGn+hoDVeISCklENJghIU59SBmzHCG\nNCIiPB1R/VJ2P4qySybz8vMpmjOH29asIWj3bvZ36cJfAwKw27dzfZs2REdEHHZyZJ82bfgiK4vP\nfHz4psLqChFp2pRASINy663w3HNOD4TmQ5S3YO5cti9cyITIyNJtrRfMmUPX9ev5c3IyXmFhcOWV\nFLZqhfn4Y04pKOCEkBAC/f0POzlyidvN+MmTGVnJMkwRadqUQEiD0qKF0wvx9NNwxx3qhSgZrgBI\nio9nQmRk6RyG3pmZtP/xR5q53WTHxhI8ciR4e5PmdhNiDJFQuh8FHL2WgxIHESlLCYQ0OH/9q1Pe\nesYMePxxT0fjGRWHK9yFhYfmMFgLK1bAN98QFBzMnW3acF5MDMPKbGhV2Q6WmhwpIlWhBEIanLAw\npzrlCy84Raaa4u+4isMVJXMYlq1cyXmpqbBxIwwaxPqePdmalFRuP4oj7WCpyZEicqyUQEiDVDIX\nYtYsZ8vvpsTtdv9huGJYdDSbWrXi1BUryPf2Jm/iRDaFhR12W+sj7WApInIslEBIg9SqFdx4o5NE\nTJ8OZXaRbvTS0tIgM5PoDh2cBmshMZErf/6Zjf7+/POMMyhyuaCw8Kj7URxuB0sRkaNRAiEN1m23\nOXtkPPcc3Hefp6OpO2FhYRAcTHJ6Or1dLvjsM/jvf3GfdBJzBg0i7q67Sq872n4U2qNCRI6XSllL\ng9W6NVx3nTOMceCAp6OpO+Hh4cQMGcKSlBSyXn8d++uvpAwfzutdu9J52DBiYmKIiYlRYiAitUoJ\nhDRod9wBWVnwz396OpLa53a7SUpKwu12c8GAAVyZmAjp6bw6fDhvdOmi/ShEpE5pCEMatLZt4aqr\nnK2+p06FoCBPR1TzKi7ZjDp4kMu//pqQVq1I/+gj/tS8ueYwiEidO+4eCGPMUGPMR8aYHcaYImPM\neRXOv1ncXvb4rMI1fsaYF40xe40xGcaYecaYJl4aSKrqrrtg/354+WVPR1I7SpdsulzcXlDANZ98\nwj6XiwW33EKL/v01XCEiHlGdIYwg4L/AXwB7mGsWAZFA6+IjrsL5WcBo4AJgGNAGmF+NmKQJio6G\nK65wNtnKyfF0NDWjZLgiKSmJpPh4RkVG0js1laD583F17Ih7wgR+XbMGt9vt6VBFpIk67iEMa+3n\nwOcAxhhzmMtyrbV7KjthjGkGXAVcbK39trhtMrDOGDPAWrvqeGOTpufuu+Hf/4Z//cspdd1QHa7C\nZKcWLeDzz+Hkk2HsWDrk5UFKCmlpaep9EBGPqO1JlGcYY1KNMb8ZY/5pjAkrcy4WJ4H5uqTBWrse\nSAEG13Jc0sh07gyXXgpPPAEHD3o6muNXdrhiWocOXBQYSO9NmwguSR7OOw+KK0cSHOws6RQR8YDa\nTCAWAZcDw4E7gNOBz8r0VrQG8qy1FRfgpRafE6mSu++G1FR4/XVPR3J8SipMjiquMBnq78+wjAz+\nuncvPwQHs7J3b9Jzc1mdmuqUnR4yRL0PIuIxtbYKw1o7t8zbtcaYNcAm4AxgSXXvP23aNEJDQ8u1\nxcXFERdXcZqFNBUnnghxcc4GW9dcA35+no6oav5QYXLtWliwgNzu3Xk+JIS2Bw/+YZdMEZHDmT17\nNrNnzy7Xlp6eXmP3r7NlnNbaLcaYvcAJOAnELsDXGNOsQi9EZPG5I5o5cyZ9+/atnWClwbr3Xnjv\nPWc+xPXXezqaY1N2S+7SCpObN8P8+dCzJ0mDB3OStdolU0SqpLK/VCcmJhIbG1sj96+zBMIY0w4I\nB34vbkoACoAzgQ+Lr+kKdACW11Vc0rh06wZ//jM8+ihMngy+vp6O6PAqTpgkOJg91rL6l1/ouXIl\nBd268dvgwSzas0e7ZIpIvXPcCYQxJginN6FkTkNnY0wfIK34uB9nSeau4uueADYAXwBYaw8YY14H\nZhhj9gEZwHNAvFZgSHXcey/07Alvvw1XX+3paA6v4pbcyenpLFu/ngkJCWyJiuK9Xr0otFbDFSJS\nL1WnB6IfzlCELT6eKW5/C6c2RG+cSZTNgZ04icN91tr8MveYBhQC8wA/nGWhN1UjJhF69IALL4RH\nHoHLLwcfH09H9EeVbcndOz+fbj/+yL7gYPjoIy4OCdFwhYjUW9WpA/EtR17FMfIY7pELTC0+RGrM\nvfc6qx7ffReuvNLT0fzRHyZM5ubCe+/hDbw7dChjQ0I0ZCEi9Zo205JGqU8fOP98pxeioMDT0fxR\n2S25KSyEefNg3z42jhrFgYgI1XcQkXpPCYQ0WvfdBxs3Or0Q9UVJiWqAmCFDWLRrF+7587GbN7P5\n7LP5sKBA9R1EpEHQbpzSaJ1yCowfDw8+CJMmeXYuRGUrLqL792doUBDh69bxfwMH8t+oKGKGDNGE\nSRFpEJRASKP24IPOcMabb8J113kujspWXPzyzjsMjI8n59pr6XH77QzVhEkRaUA0hCGNWq9eMHEi\nPPywM08MIaYaAAAgAElEQVTREyorUd07KIiLf/qJHS1bkv3gg9qSW0QaHCUQ0ug98ADs2AGvvuqZ\n55euuCgpvV48adJlDHMHDyYtM9MzgYmIVIMSCGn0unVzdup89FHIyan755dbcQGweDFs28bms84i\ns1UrrbgQkQZJCYQ0CffdB7t3w0sv1f2zw8PDnRUXqakkx8fDihXsGDyY+V5eWnEhIg2WEghpErp0\ncfbGePxxZxFEXZswcSLdBwwgaulSVkdH82r37rRTiWoRacC0CkOajHvvhbfeghdegLvuqttnB7hc\nnPfBBxR27EjgvHlMbddOPQ8i0qCpB0KajOhouPZaePJJKJmOUGf+/ndYswbX++9zQp8+Sh5EpMFT\nAiFNyt13OxMpn3qqDh/6zTfOAx95BPr2rcMHi4jUHiUQ0qS0bQvTpsGMGc7SztpSUrI6beNGZ0vQ\nP/0Jpk+vvQeKiNQxzYGQJufOO52aEPffD6+9VrP3LleyOiODiStXErR/P0Uvv0yAl/J1EWk89Cea\nNDmhoc6yzjffhF9/rdl7l5asdrm4Iz2dk7Zs4bMePViwbFnNPkhExMOUQEiTdP310KmT0xtRXSXD\nFUlJSYdKVvv4ELh4MZxyCl369iUpPh632139h4mI1BMawpAmydcXHnvM2Sfjm29g+PCq36PiDpvu\nwkIykpNpc+qpsGABBAXByJFEFxVBSgppaWlafSEijYZ6IKTJuvBCGDgQbr8dioqq/vmywxXTOnTg\nosBA9u/aRconn0ByMowdC76+Tgnr4GCVrBaRRkUJhDRZxjirKxMTYc6cqn22sh02h0VHc1abNnTf\nuJHdMTGkt2nD6tRUFqWmqmS1iDQ6SiCkSRs6FMaNc+pDVGW77z/ssAlgLdf+/jvZ3t7M7NGDmSkp\nLCgsVMlqEWmUNAdCmrzHH4eePeHZZ+GOO47tM2V32Ozt7+80rl6N39atzB8xgquefLL0OvU8iEhj\npB4IafK6dYObboKHHoLt24/tM2V32FydmsoBt5uCRYtY264ddtIkYmJiiImJUfIgIo2WEggR4MEH\nITgYbrvt2D8zYeJE2o0bx4LCQrYtXEiuMay7804NV4hIk6AEQgRo3tyZUPn++/D118f2mYCAAC65\n4gqmn3YaPbZto+j557lwyhQCAgJqN1gRkXpACYRIsUsvhdNOgylTIC/vGD+UmUnIPffA6NGEXH11\nrcYnIlKfKIEQKWYMvPgiJCXBrFnH+KGHHoK9e+H5550biIg0EUogRMro3RumTnXygm3bjnLx2rUw\ncybce69TF1tEpAlRAiFSwQMPQEjIUXbfthb+8hfo3LlqMy9FRBoJ1YEQqSA0FJ5+2pkT8dVXcPbZ\nh8653W7S0tJovXgxId9951zg5+e5YEVEPEQJhEglJk2C115zdu1cvRpcrkMbZ/mnpTHls8/YOmAA\nkUOGoDUXItIUaQhDpBLGOAlEaqqz2Va5jbM2byagoIBPw8JYMHeup0MVEfEIJRAih9Gli1Mb4uWX\n4fP305yNswoL8fv5Z1zDhzO0Y0eS4uNxu92eDlVEpM4pgRA5ghtugMGDs/no28to7tcKFi2CiAgY\nMMDZSCsz09lYS0SkiTnuBMIYM9QY85ExZocxpsgYc14l1zxkjNlpjMk2xnxljDmhwnk/Y8yLxpi9\nxpgMY8w8Y0zE8cYkUtO8vOCf/8zlYH4QC+bibJYxahR4eZGcng7Bwc7GWiIiTUx1eiCCgP8CfwFs\nxZPGmDuBKcB1wAAgC/jCGONb5rJZwGjgAmAY0AaYX42YRGrcySe34Jo/x3PR76+zsfUppLduzerU\nVBalphIzZIg2zBKRJum4V2FYaz8HPgcwptISfLcA/7DWflJ8zeVAKnA+MNcY0wy4CrjYWvtt8TWT\ngXXGmAHW2lXHG5tITZvVdjFFxk3f/YsYk/Q8AeGGmHHjtHGWiDRZtbKM0xjTCWgNlG5LZK09YIxZ\nCQwG5gL9ip9f9pr1xpiU4muUQEj9sHUrPrNmcWDqbex4qydrc2fw70fzaNlSPQ8i0nTV1iTK1jjD\nGqkV2lOLzwFEAnnW2gNHuEbE826/HcLDafbIXbz8sheffhrCvHlKHkSkaVMhKZEj+fZbmDcP3n4b\ngoO5+GL44Qe45Rbo1885RESaotpKIHYBBqeXoWwvRCTwc5lrfI0xzSr0QkQWnzuiadOmERoaWq4t\nLi6OuLi46sQtUsq9ezfBN9yAV79++EyaVNr+zDPw449w4YWQkACaQyki9dHs2bOZPXt2ubb09PQa\nu7+x9g8LKKp+E2OKgPOttR+VadsJPGWtnVn8vhlOMnG5tfaD4vd7cCZRflh8TVdgHTDocJMojTF9\ngYSEhAT69u1b7dhFKsrJccpWN3vtNcb+8AOvjhtH4PjxTJg4kYAAp3B1Sgr07QsDBsAnnzjLPUVE\n6rvExERiY2MBYq21idW5V3XqQAQZY/oYY04ubupc/L598ftZwL3GmLHGmF7Af4DtwEJwJlUCrwMz\njDFnGGNigTeAeK3AEE9aMHcuu+fNY9RPP5HXowcDo6LYvnBhubLVHTrAu+/C55/DI494MFgREQ+p\nzt+b+uEMRyTgTJh8BkgEHgSw1j4JPA+8AqwEAoBR1tq8MveYBnwCzAOWAjtxakKIeITb7SYpPp5L\nduzAu6AA3xEj6B0ZyajIyD+UrT7nHLj/fuf46isPBi0i4gHVqQPxLUdJQKy1DwAPHOF8LjC1+BDx\nuLS0NEJTU2m1ejUMGQLNmgE4ZatTUkhLSytXOOrvf4flyyEuDpYtgxNP9FTkIiJ1SyO3ImWEhYVx\n5tq1FPj5OQlEscOVrfbygvfec7bHGDnS2b1TRKQpUAIhUkb4pk303rSJxV27snrfPtIPHjxq2eqw\nMGePrYMHYfRoyMz0QOAiInVMCYQIxXMfNmwg/+abKerZk31Tp7KgsJCZKSksKCyk3VHKVkdHO0nE\nhg1w0UWQn1+HwYuIeIAKSUmTVrJkMyk+nu7r1hGzciVf33Yb4y++mOzzziMtLY2wsLBj2jCrTx/4\n8ENns87rroM33oBKd4kREWkE1AMhTdqCuXPZvnAhFwAX/vorB9q356dNm1gwdy7h4eHExMRUabfN\nM8+EN9+Ef//bmWApItJYKYGQJqtkyeaoyEh6JSfjSk+n2ZgxlS7ZrIpLLoEnn3TqQzz8cA0HLSJS\nT2gIQ5qstLQ0yMykY6tW8N13TmnJiAiiDx6sdMlmVdx+O+Tlwb33Qna2k0xoOENEGhMlENJkhYWF\nQXAwuUuXQlERnHEGcPglm1V1zz0QEADTpztJxMyZSiJEpPFQAiFNVnh4OKfExBD22mvs6tePAG9v\nkkuWbI4bd9y9D2XdequTRPzlL5CTAy+9pH0zRKRxUAIhTdrolSvJbd6c13v2JD8lBYKDiTnKks2q\nuvFGJ4m4+moniXjjDfDW/3ki0sDpjzFpupYvx3v+fLzffJMbxo6t0pLNqrrySieJuPRS2LUL5sxx\nClCJiDRUSiCkabIWbrvNKd5w2WWEu1y1kjiU9ec/Q6tWTqGpAQPgo4/gpJNq9ZEiIrVGo7HSNC1Y\n4Ox+9fTT4HLV2WOHD4cff3R6IwYNgo8/rrNHi4jUKCUQ0vTk5cGddzolI886q84f37mzk7uceSaM\nGwePPup0iIiINCRKIKTpeekl2LLFqfbkISEhMH++U63ynnvgvPO0k6eINCxKIKRp2bcPHnrIWRLR\ns6dHQ/HyggcfdOZCrFrlhPPhhx4NSUTkmCmBkKbl4YchN9f5zV1PjB0La9bAkCEwYQJMngwHDng6\nKhGRI1MCIU3Hxo3w/PPwt79BVJSnoyknIsLpfXjjDWdoo3dvWLzY01GJiByeEghpOu68E1q3dspD\n1kPGOL0Pq1dDx45w9tlwwQXOdA0RkfpGCYQ0Dd995yzdfOwxZw1lPdaxIyxZAu++CytXQvfuzmTL\nrCxPRyYicogSCGn8ioqcXof+/SEuztPRHBNjYNIkWL/eqXf11FPQrRu88w4UFno6OhERJRDSFLzz\nDiQkONthNrCdrIKCnHmf//ufk/9cdpnTI/Hvf0N+vqejE5GmrGH9aSpSVVlZcPfdTv3oIUMAcLvd\nJCUl4Xa7PRzcsevc2RmB+ekn6NHDmSvRtSu8+qpTF0tEpK5pLwxp3J55BvbsgccfJycnhwVz55IU\nHw+Zmc7Om0OGMGHiRALq+byIErGxzmqN1audnonrr3dWpF5/PVxzTb1bXCIijZh6IKTx2rkTnngC\nbrkFOndmwdy5bF+4kAkuF9M6dGCCy8X2hQtZMHeupyOtst69Ye5c+PVXpyL3449Dhw4wcaIzAVOl\nsUWktimBkMbrjjucSQR33+0MW8THMyoykt6RkYT6+9M7MpJRkZEkxcc3qOGMsk46yRnG2LEDZsxw\nClINH+60P/ywU/pCRKQ2KIGQxik+3lkH+dhj0Lw5aWlpkJlJdGhoucuiQ0MhM9M534A1bw5TpzqT\nLZcscYY6Hn8cYmKcyZfPPAPbtnk6ShFpTJRASONTWAhTpji/OSdPBiAsLAyCg0lOTy93aXJ6OgQH\nO+cbAWPgjDOchSe7dzvDHB06OBt2degA/frBfffBihVaDioi1aMEQhqf116D//7XKVtdvGwzPDyc\nmCFDWJSayurUVNIPHmR1aiqLUlOJGTKE8PBwDwdd8wIDncUn8+c7ycTbb8OJJ8ILL8DgwU5Rzssu\nc5aEbt6seRMiUjVahSGNS1qa89ftK6+EgQPLnZowcSILgAXx8ZCS4qzCGDeOCRMneiTUutSsGVx6\nqXMUFDgVLj/7zDnefddJHtq1g2HDnOPUU516E976E0JEDsPYBvbXDmNMXyAhISGBvn37ejocqW+m\nTHH+qr1hA0RGVnqJ2+0mLS2NsLCwRtnzUFVpac6Uke++c46EBGd4IzAQTjnFGfbo1w/69nXmVPj4\neDpiETleiYmJxMbGAsRaaxOrcy/9/UIaj19+gZdecuo+HyZ5AGc4Q4nDIWFhzpbiY8c67zMznSTi\np5+c45NP4NlnnXO+vk4Bq549naNHDyep6NwZ/P099x1EpO4pgZDGwVq4+Wbnt9vUqZ6OpkELDobT\nT3eOEmlpTvGqX3+FtWudn4sWwf79znljoH17J5mIiYHo6PJHVFSDqyIuIkdRqwmEMeZ+4P4Kzb9Z\na08qc81DwDVAcyAeuNFaq9XrUjVz5jj9719+qT72WhAW5qzuOOOMQ23Wwq5dTq2JpKRDx/Llzr+O\nkuQCnH8lUVHQpk35IzISIiLKH4GBdf3tROR41EUPxK/AmYApfl9QcsIYcycwBbgc2Ao8DHxhjOlu\nrVWFfzk2+/bBX/8KF14IZ5/t6WiaDGOcpCAqCoYO/eP5AwcgOdmZr5qc7BQGLTm+/dYpflVZ+Y2A\nACdhCQ93fpYcoaFOvYuyP0NCDh3Bwc7PoCBwuWr/+4s0dXWRQBRYa/cc5twtwD+stZ8AGGMuB1KB\n84GGV19YPOOOOyA3F557rlyzJkt6VrNm0KuXcxxOfj7s3essMy053G7nSEs79HrbNkhPd3o19u8/\n+gZifn5OIlFyBAQ4PRsBAeUPf3/nWn//Q699fZ2fJa9LDh+fQz/LHt7e5X+6XM7rksPlOtRW8trl\nchIwkYasLhKIGGPMDuAgsBz4m7V2mzGmE9Aa+LrkQmvtAWPMSmAwSiDkWHz3nVP34eWXS3eSagyb\nZjUVJUMbVd0E7OBBJ6HIzISMjPJHVpZzZGcfep2T47zPyXEOt9vJOQ8eLH/k5jrJSW6uc9TmIjVj\nnETCy+tQUuHldeh9yevK3htz+J8VXx+prbKjJLbK2sqeq9h2uHOVXXe0toqva6utMl26OFVc5ehq\nO4FYAVwJrAeigAeA74wxPXGSB4vT41BWavE5kSPLzYXrrnO26b722tLm0k2zIiOJ7tCB5PR0Fi1c\nyALgkiuu8Fy8UmNKegyOsNimRhQUOAlFfv4ff+bnO+dLXufnO8tfCwoOHSVtJUdBQfn3RUXlX5e8\nL/va2kPvK7aVPVe2reLrI7WVJElHel+xreR1xZ8V2yq77nCfKetIbZXd52iO9RngJKVybGo1gbDW\nflHm7a/GmFVAMjAR+K02ny1NwGOPOSUU588vneJfsmnWhOJNswB6F68vXBAfj3vMGA1nyDErGYYQ\nkT+q0/81rLXpxpgNwAnAUpyJlZGU74WIBH4+2r2mTZtGaIWNkeLi4oiLi6uxeKUeW7cOHn0U7rrL\nKUZQrHTTrA4dyl0eHRoKKSmkpaUpgRCRJmH27NnMnj27XFt6hf2AqqNOEwhjTDBO8vCWtXaLMWYX\nzgqN1cXnmwEDgRePdq+ZM2eqEmVTVVQE118PnTrB3XeXO1V206zeZSobNbZNs0REjqayv1SXqURZ\nbbVdB+Ip4GOcYYu2wINAPjCn+JJZwL3GmI04yzj/AWwHFtZmXNLAvfYafP+9s291hfKHpZtmLXT+\nE4oODXXmQKSmEjNunHofRERqSG33QLQD3gPCgT3AD8Aga60bwFr7pDEmEHgFp5DU98Ao1YCQw9q0\nCW69Fa65pnxVozKa8qZZIiJ1RZtpScNRUODUV961y9muOySk9FRlNR9UB0JEpDxtpiVN0xNPwIoV\nzvBFcfJwpJoP2jRLRKT2aHsbaRh++gkeeAD+9jc49dTS5tKaDy4X0zp0YILLxfaFC1kwV3XIRERq\nkxIIqf+ys+Gyy6B3b7jvvtLmkpoPo4prPoT6+9M7MpJRkZEkxcfjdrs9GLSISOOmIQyp/+68E7Zu\nhcRE8PUtnduwb98+1XwQEfEQJRBSv33xBbzwAjz3HDkdO7LgrbdK5ztke3uzY+dONgYFEdu2belH\nVPNBRKT2KYGQ+mv7drj8chgxAm66iQVvv/2HPS4ey8jgjcREfLy9VfNBRKQOKYGQ+ik3Fy680Nk/\n+e23ce/bV+keF9P69uXJpCTey8rCv7jn4Wg1HwqKCsgtyCW/KJ/8wnzyi/IpKCogPCCcIN+guvqG\nIiINmhIIqZ9uvhl+/hl++AEiIkhLSqp0vkPX8HB6ZmVx7pQptGjR4g81H5L3J/NL6i+sSV3D6t2r\nWZO6hg3uDRTawkofGxEUQcfmHenUvBMdm3ekT2Qfzuh4BlEhVdxvWkSkkVMCIfXPa6/Bv/7l/Ozf\nHzj6HhddunQhPDwcay1rUtcw73/z+OB/H7Bu7zoAWvi3oFdkL4Z3Gs7UAVNp5tcMH5cPvi5ffLx8\ncHm52J21m637t7J1/1a27N/Ciu0reCL+CQC6hnfljI5ncEbHMxjeaTgRQRF1/89FRKQeUQIh9cvK\nlXDTTc5mWVdfXa6a5JH2uNhv9jPj6xnMWzePDe4NhPqFMq7bOB4Z/gj92/anbUhbjDFVDic1M5Vv\nk79l6dalLN26lFcSXsHLeDG803Au7nExE7pPoEVAi5r+pyAiUu+plLXUH6mpEBsL7duTs2gRCxYu\nLFdhMrp/fzCG5FWrStu8+kXx32ZrWLhhIaF+oZzf7XwuPOlCzup8Fr4u3xoPcVfmLj5e/zFz1s5h\nyZYleHt5M/KEkUzqNYnx3cbj5+1X488UEakpNVnKWgmE1A8HD8I558D69ZCQwLuLF7N94UJGRUaW\n621oN24cI0afy4e/fshbSW+xbOcyurTowm2n3sYVfa4gwCegzkL+PeN3PvjfB8z5dQ7Lty+nZWBL\nrjr5Kq6LvY4uYV3qLA4RkWOlvTCkcSkogEmTYNUqWLwYt79/pSsuAJ75cS5PZs1g9Z7VDGw7kHkX\nzeP8bufj8nLVedhRIVHcPPBmbh54M7/t/Y1XfnqFfyX+iyeXPck5Xc7hhn43MPbEsR6JTUSktqmU\ntXiWtXDDDfDRRxx47TWSIiLYtGmTs+IiNLT0srXs5vZWX/KfVp/hVeTFkiuWsPzq5Vxw0gX14hd0\nt5bdmDlyJjtu3cGb495k38F9jH9/PCe+cCLPrXyOjNwMT4coIlKjNIQhnnXXXfDEEyy79lq+hEMV\nJlNSuDUmhjZtm3EfS3iD/9K2MISB6f156b45tGzZ0tORH9WPO35k5oqZzF07l2DfYK7tey1TB06l\nQ2iHo39YRKQW1OQQhnogxHOeegqeeIKf4uL4fu/e0h01Lw0KoiDzAFPTF3GCfZ75rOO2jMFc/78/\ncX7sZQ0ieQDo37Y/713wHltu2cL1sdfzauKrdH62M5PmTyLx92r9fysi4nFKIMQz3ngD7riD7Ftv\n5ZPg4HI7arojs/lx5BaWn7SD7gc6c93G0QRk9KDjeeOPWGGyvmof2p4nzn6C7bduZ+Y5M1m+fTmx\n/4rlzP+cyaKkRTS0XkAREdAkSvGE115z6jxcfz07rr8eHniA6A4d2MEBpvMl75u1DPJqy4jkQVx1\n44OVVphsiIJ9g5k6cCo39r+RBesW8NSypzj3vXPp0aoH0wdPZ1KvSVoGKiINhnogpO5YC48+Ctde\n60ycfPFFwsLDKQwO5KH8b+nGiyxhK2/Z83l5zxha+7SnS5cuxMTENPjkoSxvL28m9pjIqmtW8e2V\n39KpRSeu+ugqomdF8/B3D+POdns6RBGRo1ICIXWjqAhuuQXuuQcefBD3gw+StHkzy3Ys498dP2dm\nyArGZp/IqoPXcHJqaz5P3U3MkCGNKnGoyBjDsOhhfBz3Mb/d9BvndzufR75/hPYz23PjJzeyfu96\nT4coInJYWoUhtS8vD664At5/n7xZs/ggNJSfl3/FV4HxrA7dShefzlza/BL4767SCpMxQ4YwYeJE\nAgLqrjBUfbAnaw8v//QyL/z4AruzdnNOl3OYOmAqo2JG4WWU74tI9aiQlDQc6ekwcSIsXQoffMDs\njP18EP8030VtwMe4uD/9dPy2tqDDuC6MfOyW0n0vGnPPw5G0CmrF30//O3cMuYP3177P86ueZ8zs\nMXRp0YWb+t/E5FMm09y/uafDFBHREIbUooQECk8+mcLly0l//33ei85l+sbpLGqzlstMbzYwhQea\nncHoyNbOnhfQ6OY7HC8/bz8u73M5q65ZxfKrlzOw3UDuWHwHbZ5pw+SFk1m+bblWb4iIRymBkJpn\nLXkzZlA4aBC7MjN58OxYYldO4ZJPLyEk148fcifzIqMJJxBwdtYkM5O0tDQPB17/GGMY1G4Q7054\nl5S/pnDP0HtYsmUJp75xKr1f7s1zK59jX84+T4cpIk2QEgipMW63m00JCeSedx6+06ezolMbpl8e\nzqO9viXPO4OJG2IZ/GM3gvaX3yUzOT0dgoMJCwvzUOQNQ1RIFPcMu4fNt2zmi0u/oGt4V6Z/OZ2o\nZ6K46IOL+Gj9R+QV5nk6TBFpIjQHQo6L2+0una8QGBjIgrlzyZ03j/FLllCYl8eUM8N49bQdNMOP\nJzmbm1z9Wd/MzUM71vF/KSkA5XbZjBk3TkMXx8jLeDGiywhGdBnBrsxdvLP6Hd5e/Tbj5owjPCCc\ni3tezGW9L2NA2wEYYzwdrog0UlqFIVWSk5PDgrlznTkLxSsmcrOymPTdd/Tavp2tbYMZOSGHTSFF\nXLn7BGa0vYAQnOJI6QcP8uTmzYT160fGxo1NfsVFTftl1y+8vfpt3lvzHr9n/k50aDQXdL+AC0+6\nkIHtBmoVh4hoFYZ4zoK5c9m+cCETIiOJbteOPd9+S9RPP1LgY7hxrIv/nJLDNQWxZH+yj0F5Llxj\nCyme6kByejo+zZtz5dVXAzT5FRc1rU/rPvRp3YcnznqCpVuXMn/dfN5d8y4zVsygbUhbLuh+AWO7\njmVY9DB8Xb5Hv6GIyBEogZBj5na7SYqPZ0JEBL2zs8n6cC6dU9N4NRYeHWq4NHAgW8ypRPgE8Wq7\nH/l+9Wpitm2jT3R0pUMVShxqh8vLxZmdz+TMzmfy/Kjnid8Wz7z/zWP+uvk8t+o5gn2DObvz2YyO\nGc25MecSFRLl6ZBFpAFSAiFHVTLfYV9aGh03rKXVxrWwI521beDRq4M5uDeUOxf58eez+xIeHgRA\nnzZt+CIri898fPgmJcUZqhg3rkFuhtWQubxcDIsexrDoYTw78ll+Sf2FTzd8ymcbP+O6T66jyBbR\nK6IXwzsNZ3in4QyLHqY6EyJyTJRAyGGVznf4/jtabVrJmWuTuHJPHvHt4eVJbTjlhNOYb7ry5OYf\nWJW7hdiMDLyDgkhOT2eJ2834yZMZOWaMhirqCWMMJ7c+mZNbn8w9w+5hb/Zevtz0JV9v/pr/++3/\neHbls3gZL/pG9eX06NM5tf2pDG43WD0UIlIpTaKUckp6G5o1b8Z//vUALT56n7M37Cd6v+W7Ti7e\nim5OQW5rpvcZWrqK4qOdO/k9IoJWxmhiZAO2Zd8WvtnyDd9s/YYfUn4gJd1ZLdOxeUdObX8qA9sO\nJDYqlj6t+xDsG+zhaEXkeGgSpdSIiksx3579Jt//OI/oLWs4c72b27daMn1h9UktSes3kKFtY/Hf\nsZMnk5J4LysL/+L6DTETJjB94kSys7PV29CAdWrRiatbXM3VfZ1JrjsO7GD59uUs27aMZduWMe9/\n88grzMNg6NqyK7FRsZzS+hR6RvSkR0QP2oa01bJRkSZECUQ9VPYXe8kv4ppsKxmaWLNsCclFW8jy\n2k1MyjZO35zFy1sgKB9+bRvIE+1dDO0+iFMHn1oaW9fwcHpmZXHulCm0aNGi3HMCAgKUODQibZu1\n5cKTLuT/27v34LjK847j32d1Wd1sSZZsyWD5BhiwCTb3S7gkgRacDGGgjMHxNE2ZDtCQTuopQyBp\nmhQmQGlCCkkoDG25BZzSNhRMQg0U2lIINcSJC8bYpmBA2JItW5YsrS57efvHe1ZerWVbx155dcTv\nM3Pm7L7n3fXzaNe7zzl73vdcMf8KAJLpJOu2r2PN1jVDy5PvPEkimQCgNl7LgmkLWDB1AfMa5nHM\nlGOY1zCPufVziZfGi5mKiIyBcVFAmNn1wA1AM7AW+BPn3OvFjerwGGlCptw5FmaddhqY8cHq1YfU\nNswSnHEAABAFSURBVNjbydb6fpLzJrF950amvL2OhVv6WfouLGyHlMH7jZUMnn8y1fNPYXo8TuvK\nlaz/8CMWLUxQVeXHYmZnjTzqqKNULHzClJWUDZ1DcfVJVwOQcRk279rMum3reGvbW6zbvo7Xt7zO\n428+Tm+yF/ATX82sncmcujnMqZvD7LrZQ0tLbQtHTDpCw0pFIqjoBYSZXQn8ALgGWA0sB1aZ2Tzn\nXEdRgxtDI03ItN05jmxv5/Ijj2TWzJl80NXF7ffdxxTg2tNPP2BbHY7Lz1xAcrbxWv8G7nj5KVxF\nH3MyaWb0JPj0ejj7H2F+8Fftq46TOmo2WxfO4u5Nmzhn8mSOPOlsqKqiAThx7twDDsWUT7aYxZhb\nP5e59XO55NhLhtqdc7T1tLFxx0Y27dzEph2b2Ny1mTe3vcnKjSvZntg+7HmmVU9jxuQZzJg8g+k1\n02mqbqK5ppmmGr+eVj2NxqpGauO1+plEZJwoegGBLxjud849AmBm1wFfAK4G7ixmYGNp2IRMM2ey\ntr2du156ic/OmcOJTU0AHJnJ0JhIsBBH06Q4HRUJWuli15Q2YhVpHp+ylvayPt4t28GWM7bS3J+i\nq+03zFsHi7bB77fB3F1Q4sAB/U0NtNXVcY91c9XixUybPRvMGEwkiLW10d7bS19f39DRBg3FlINl\nZkyfNJ3pk6Zz/uzz99reO9jLB10f0NrdSmt3Kx93f8zHuz+mtbuV17e8TntPO+297aQyqWGPK42V\n0lDZQGNVIw1VDdRX1FNfWe/XFfXUVdRRW1FLbbyWyfHJ1Fb49aTySdSU11BVVqUCRKRAilpAmFkZ\ncApwW7bNOefM7AXgrKIFNsayEzJd1FzPy9M+4BcMsHnaLt47exc/qVzLg+mP6C5JsrW8m56Lu/m3\nVIb7dm2guRWO2A1nJGHGDpi9oYPZXcaRnWkqk3tG06SqK+mvn8JaEiSa4sw+5xxqjz6ayngcOjt5\n/plnWDQ4yLTgg7ShqopJDQ283N3NCRqKKYdBdXk186fOZ/7U+fvsk3EZdvbtHComdiR20JHoYEef\nX3ckOujs72RDxwY6+zvp7Ouks79zvxcUM4ya8hpqymuoLq+mqqyK6jK/zi6VZZVUlFRQWVZJZWkl\nFaUVxEvjfl0SH7pdXlJOeUk58ZL40O3yknLKSsooi5XttS6NlVJW4telsVJKrETFjERasY9ANAIl\nQHteeztw7OEP5/DYuXMn9PQwbVoNfc8/y/SBUhb0G0t609T3O6akeqjtd0xKpCnNDB9m68zYVRIj\nVVpG/YwZlM6so3d+DQ9s2kRDWRkXX3opVXV1DCQS/GzlShYBx8yZA3F/EtvuwUHiLS38R3c3de3t\nQ0Mx05MmUX3hhTxrxrN5Rxt0cqQUQ8xiNFY10ljVyAIWjPpxA6kBuge66Rro8uv+LnoGe4aW3YO7\n6RnsoXewl0QyQW9yz7oj0UFfqo/+VD99yb6h2wOpAfpT/fSn+nEUbuh7iZXsKShiJZRYyT7XMYvt\ndTtmsaElW5DktsUshuHbstty7+ffDrsGhrWBL9Ly7+f2y8rfntuW/7jRtOUaqTAbbb+WyS0sP2v5\nXu2yt2IXEAdt+fLl1NbWDmtbunQpS5cuLVJEozdlyhSoqaGva5Ab1tdDRQXE42zo6WF7Xx9HtLRQ\n2dBAuxn/0toKpaVcvGgRzU1NbE4muf2NN/acAxEUAC92dzMFOHpggFn9/Xy8ezcdVVWsAU7bvZtZ\nsdjQOQxfWLaM8vJyfv7KK5AtFjQUUyaIeGmcqaVTmVo9teDP7ZwjlUnRn+onmUkykBpgMD3IYHqQ\ngfQAyXSSZCbJYHqQZNqvU5kUqUyKZCbp1+kkaZceak9n0iQzSdKZNGmXJp3x27K30y5NxmWGbXc4\nMi4zrN05R4aMX7vMnragn8MNPTbb7nAkM0mcc0Pto10D+72d/Xvtb1tuMZY7J1H+43Pb8l+P/T1m\ntI/NdcK0E1jOxCggVqxYwYoVK4a1dXV1Fez5izqRVPATRgL4Pefc0zntDwG1zrnLRnjMhJhI6rGH\nH6b1qadY3NS03wmZCjUKY6QJnkYa8ikiIhPXhJlIyjmXNLNfAxcATwOYP6Z0AXBPMWMba5cvWcLP\nYdRHAXZcdllB28BfzEqFg4iIHIyiT2VtZkuAh4Dr2DOM8wrgOOfc9hH6T4gjEFk6CiAiIofLhDkC\nAeCce8LMGoFbgCbgt8BFIxUPE5GOAoiISBQVvYAAcM7dC9xb7DhERERkdGLFDkBERESiRwWEiIiI\nhKYCQkREREJTASEiIiKhqYAQERGR0FRAiIiISGgqIERERCQ0FRAiIiISmgoIERERCU0FhIiIiISm\nAkJERERCUwEhIiIioamAEBERkdBUQIiIiEhoKiBEREQkNBUQIiIiEpoKCBEREQlNBYSIiIiEpgJC\nREREQlMBISIiIqGpgBAREZHQVECIiIhIaCogREREJDQVECIiIhKaCggREREJTQWEiIiIhKYCQkRE\nREJTASEiIiKhqYAQERGR0FRAiIiISGgqIERERCQ0FRAiIiISmgqIcWDFihXFDqFgJlIuoHzGs4mU\nCyif8Wwi5VJIY1ZAmNlmM8vkLGkzuzGvT4uZ/cLMes2szczuNLNPXFEzkd6cEykXUD7j2UTKBZTP\neDaRcimk0jF8bgf8OfAAYEHb7uzGoFD4JbAFOBM4AngUGAweJyIiIuPUWO/t9zjntjvntgVLX862\ni4DjgGXOuTedc6uAbwPXm9lYFjYiIiJyiMa6gLjJzDrMbI2Z3WBmJTnbzgTedM515LStAmqBBWMc\nl4iIiByCsdzTvxtYA+wEzgbuAJqBG4LtzUB73mPac7at3cfzVgCsX7++kLEWVVdXF2vWrCl2GAUx\nkXIB5TOeTaRcQPmMZxMpl5zvzopDfS5zzo2+s9ntwDf208UBxzvnNo7w2K8A9wM1zrmkmd0PzHTO\nLc7pUwn0AouDnzRGiuFLwGOjDlpERETyLXPOPX4oTxD2CMT3gQcP0Oe9fbSvDv692cAmoA04La9P\nU7Bu28/zrwKWAZuB/gPEIiIiIntU4L+HR9xJDyNUAeGc2wHsOMh/6yQgA2wL7v8K+KaZNeacB/G7\nQBfw9gFiOKSqSURE5BPs1UI8yZicA2FmZwJnAC/hh26eDdwFPOqc6wq6PYcvFB41s28A04FbgR87\n55JjEZeIiIgURqhzIEb9pGYnAfcCxwJx4H3gEeCHucWBmbUAfwt8Bn/uw0PAzc65TMGDEhERkYIZ\nkwJCREREJrZP3LTRIiIicuhUQIiIiEhokSogzOybZvZKcPGtnfvoE5kLdJnZ9Wb2vpn1mdlrZpY/\nrHVcMrNzzexpM/s4uFDaF0foc4uZbTGzhJk9b2ZHFyPWAzGzm81stZl1m1m7mT1pZvNG6BeVfK4z\ns7Vm1hUsr5rZxXl9IpFLPjO7KXi/3ZXXHol8zOw7eRcYzJjZ23l9IpFLlpkdYWaPBjMOJ4L33sl5\nfSKRU/BZnP/6ZMzsRzl9opJLzMxuNbP3gljfNbO9rjF1qPmMyy/W/SgDnsCfeLmXnAt0leKnyv4D\n4CvALYcpvlEzsyuBHwDfwQ9xXQusMrPGogY2OtXAb4Gv4icPGyYYVfM14BrgdPwJsqvMrPxwBjlK\n5wI/wo8auhD/HnsumNQMiFw+H+EnezsZOAV4EXjKzI6HyOUyJCiuryFvhtoI5vMWfr6b5mA5J7sh\narmYWR3wCjCAv7bR8cCfAZ05faKU06nseV2agd/Bf749AZHL5SbgWvxn9HHAjcCNZva1bIeC5OOc\ni9yCLwx2jtC+GEgCjTlt1+Lf0KXFjjsv1teAu3PuG9AK3Fjs2ELmkQG+mNe2BViec38y0AcsKXa8\no8inMcjpnImQTxDvDuAPo5oLUANsAD6HHxp+VxRfG/zOwpr9bI9MLkF8dwD/eYA+kcopL/a/ATZG\nMRdgJfBAXts/A48UMp+oHYE4kEhcoMvMyvB7h/+ebXP+FXwBOKtYcRWCmc3BV++5uXUD/0M0cqvD\n73XshGjnExzGvAqoAl6NcC4/AVY6517MbYxoPscEP/39n5n91PxQ9qjmcgnwhpk9Efz8t8bM/ii7\nMaI5AUOf0cuAvw/uRy2XV4ELzOwYADNbCHwaf4S+YPlMtMtmH+wFug63RqCEkWM99vCHU1DN+C/g\nkXJrPvzhjJ6ZGX6v47+dc9nfpiOXj5mdgJ/ptQI/kdtlzrkNZnYW0cvlKmAR/vByvqi9Nq/hf1Ld\ngJ8477vAfwWvV9RyAZgL/DH+p9jv4Q+D32NmA865R4lmTlmX4Xc8Hw7uRy2XO/BHFN4xszT+dIVv\nOed+FmwvSD5FLyDsEC7QJVJg9wLz8ZV6lL0DLMR/AF4BPGJm5xU3pPDMbAa+oLvQTYDZad3wCwS+\nZWargQ+AJfjXLGpiwGrn3LeD+2uDYug64NHihVUQVwPPOuf2d12m8exK4EvAVfgZnxcBd5vZlqC4\nK4jx8BPG9/EneexrOZ59X6ArXxt7LsiVNZoLdB1uHUCakWMdT3EejDb8+RyRys3Mfgx8HviMc25r\nzqbI5eOcSznn3nPO/cY59y38kbevE71cTgGmAmvMLGlmSeB84OtmNojfW4pSPsM4P63/RuBoovfa\nAGwF1ue1rQdmBrejmBNmNhN/QvUDOc1Ry+VO4A7n3D8559Y55x4DfgjcHGwvSD5FLyCcczuccxsP\nsKRG+XS/Aj6VN5LhgBfoOtyCvalfAxdk24LD5xdQoIucFItz7n38GzA3t8n4UQ7jMregeLgU+Kxz\n7sPcbVHMZwQxIB7BXF4APoXfe1oYLG8APwUWOufeI1r5DGNmNfjiYUsEXxvwIzDyf3I9Fn9UJcr/\nd67GF6e/zDZEMJcq/E5qrgzBd37B8in22aIhzyxtwX+I/AW+KMh+qFQH22P4va1ngRPxQ4vagVuL\nHfsIuSwBEsCX8Uda7sefLT+12LGNIvbq4O++KHhT/mlwvyXYfmOQyyX4L4B/xV/CvbzYsY+Qy734\nUTrn4qvv7FKR0ydK+dwW5DILOAG4HUgBn4taLvvIL38URmTyAf4aOC94bc4Gng8+nxqilksQ76n4\nIZw3A0fhD5nvBq6K4usTxGvAZuB7I2yLTC7Ag8CH+KOqs/DndGwDbitkPkVP9CD+KOkRlvNy+rQA\nzwA9wX/OvwJixY59H/l8NXiz9uGPnpxa7JhGGff5+MIh/3X4h5w+38UPE0rgR8IcXey495HLSHmk\ngS/n9YtKPn+H/8mvD7+H8RxB8RC1XPaR34vkFBBRygdYgR+q3Rd8uD8OzIliLjnxfh743yDedcDV\nI/SJTE74uR/S+4oxKrngd/Luwl/IsjcoDP6SvOkMDjUfXUxLREREQiv6ORAiIiISPSogREREJDQV\nECIiIhKaCggREREJTQWEiIiIhKYCQkREREJTASEiIiKhqYAQERGR0FRAiIiISGgqIERERCQ0FRAi\nIiIS2v8D9z0ZIsM6g7cAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x158ef3cc0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.plot(t,state[:,0])\n",
"plt.plot(t,state[:,1])\n",
"plt.plot(t,state[:,2])\n",
"plt.scatter(t,data[\"Recovered\"],alpha=0.5,c=\"red\")"
]
},
{
"cell_type": "code",
"execution_count": 99,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.002378053866825905, 0.2950335336536682, 241.3428258354608]\n"
]
}
],
"source": [
"print([np.exp(loga),np.exp(logb),np.exp(logS0)])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [default]",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment