diff --git a/examples/toy/100 b/examples/toy/100 new file mode 100644 index 000000000..16eeae493 --- /dev/null +++ b/examples/toy/100 @@ -0,0 +1,15 @@ +Timer unit: 1e-06 s + +Total time: 5.4e-05 s +File: +Function: _propensities at line 1 + +Line # Hits Time Per Hit % Time Line Contents +============================================================== + 1 def _propensities(xs, ks): + 2 return [ + 3 1 37.0 37.0 68.5 1e5 * xs[0] ** 2 * ks[0], + 4 1 10.0 10.0 18.5 xs[0] ** 3 * ks[1], + 5 1 1.0 1.0 1.9 2e5 * ks[2], + 6 1 6.0 6.0 11.1 xs[0] * ks[3] + 7 ] \ No newline at end of file diff --git a/examples/toy/model-logistic.ipynb b/examples/toy/model-logistic.ipynb index 0ed958e51..3a3540924 100644 --- a/examples/toy/model-logistic.ipynb +++ b/examples/toy/model-logistic.ipynb @@ -17,7 +17,10 @@ "cell_type": "code", "execution_count": 1, "metadata": { - "collapsed": true + "collapsed": true, + "jupyter": { + "outputs_hidden": true + } }, "outputs": [], "source": [ @@ -200,9 +203,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.7" + "version": "3.7.4" } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/examples/toy/model-michaelis-menten.ipynb b/examples/toy/model-michaelis-menten.ipynb new file mode 100644 index 000000000..0c9ed5cc0 --- /dev/null +++ b/examples/toy/model-michaelis-menten.ipynb @@ -0,0 +1,287 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Michaelis-Menten\n", + "\n", + "- X1+X2 -> X3 \n", + "- X3 -> X1+X2 \n", + "- X3 -> X2+X4" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import pints\n", + "import pints.toy.stochastic\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import math" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import line_profiler" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Specify initial concentration, time points at which to record concentration values, and rate constant value (k)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "x_0 = [1e4, 2e3, 2e4, 0]\n", + "model = pints.toy.stochastic.MichaelisMentenModel(x_0)\n", + "\n", + "times = np.linspace(0, 24, 100)\n", + "k = [1e-5, 0.2, 0.2]\n", + "\n", + "values = model.simulate(k, times)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEGCAYAAACtqQjWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de7QU5Znv8e8DottrBIxyv+hgkEsGw444R5NxjTEimqBnibeMisma7TIyMRn1QExm6SRqMI5GczERj3LxiKgxGRmjYzgsjcYLujEcEZGEUYEtBAhgwBhE5Tl/VDUUm+7d1d1V3V3dv89ae7H329Vdb2VHfrzv89Zb5u6IiIgkoVutOyAiIo1DoSIiIolRqIiISGIUKiIikhiFioiIJGafWneg2g477DAfMmRIrbshIpIpixcv/pO7f7zYcU0XKkOGDKG9vb3W3RARyRQzWxXnOE1/iYhIYhQqIiKSGIWKiIgkRqEiIiKJSS1UzGygmT1pZsvNbJmZXRG2X2dmb5vZkvBrQuQ93zSzlWa2wsxOjbSPD9tWmtm0SPtQM1tkZn8wswfMbN+0rkdERIpLc6TyIXClux8DHA9cbmYjwtd+4O5jwq/HAMLXzgNGAuOBO8ysu5l1B34CnAaMAM6PfM5N4WcNA7YAX0nxekREpIjUQsXd17n7y+H324DlQP8u3jIRmOfu77v7m8BK4Ljwa6W7v+HuO4B5wEQzM+AfgJ+H758NnJnO1YiISBxVuU/FzIYAxwKLgBOAKWZ2EdBOMJrZQhA4L0Te1sHuEFrTqX0c0Bt4x90/zHN85/O3AW0AgwYNqvyCRETq3JYHHmTro4/u+nm/Y4bT55prUj9v6qFiZgcBDwNfd/etZvZT4LuAh3/eAnwZsDxvd/KPpryL4/dudJ8BzABobW3VA2REpCFFg+S9l14C4IBPf7qqfUg1VMysB0Gg3OfuvwBw9/WR1+8CclHaAQyMvH0AsDb8Pl/7n4BDzWyfcLQSPV5EpCkUCpIDPv1pDjnjDHqee05V+5NaqIQ1j7uB5e5+a6S9r7uvC388C3g1/H4+MNfMbgX6AcOAFwlGJMPMbCjwNkEx/wJ3dzN7EjiboM5yMfBIWtcjIlJLnaezcuohSKLSHKmcAFwILDWzJWHbNQSrt8YQTFW9BVwK4O7LzOxB4DWClWOXu/tHAGY2BXgC6A7c4+7Lws+bCswzs+uB3xGEmIhIQ4gznVUPQRJlzfaM+tbWVteGkiJSr7oKklqGh5ktdvfWYsc13S7FIiL1ICvTWaVSqIiIVEkWp7NKpVAREUlJ59FI1kchcShUREQS1NVopFGDJEqhIiJSoXq7V6SWFCoiImVQkOSnUBERiUlBUpxCRUSkCwqS0ihUREQ6UZCUT6EiIk2rUW9ArCWFiog0lWa4AbGWFCoi0vA0nVU9ChURaUgKktpQqIhIw1CQ1J5CRUQyqxn31qp3ChURyZRm31ur3ilURKTuaVorOxQqIlKXFCTZpFARkbqhIMk+hYqI1I2tjz7K9tdfp2X4cAVJRilURKSmoqOTXKAMvndOjXsl5VKoiEjVFZrmahk+nEPOOKOWXZMKKVREpCpUL2kOChURSY2CpPkoVEQkUQqS5qZQEZGKKUgkR6EiIhXTUmDJUaiISFm0FFjyUaiISGxaCizFKFREpEuql0gpFCoisgc9o0Qq0S2tDzazgWb2pJktN7NlZnZF2N7LzBaY2R/CP3uG7WZmPzSzlWb2ipl9KvJZF4fH/8HMLo60jzWzpeF7fmhmltb1iDSLXNE954BPf5o+//ZvDL53DoPvnaNAkS6lOVL5ELjS3V82s4OBxWa2AJgMLHT36WY2DZgGTAVOA4aFX+OAnwLjzKwXcC3QCnj4OfPdfUt4TBvwAvAYMB54PMVrEmlIKro3iPaZsPTn+V/rMxpOm556F1ILFXdfB6wLv99mZsuB/sBE4KTwsNnAUwShMhGY4+4OvGBmh5pZ3/DYBe6+GSAMpvFm9hRwiLs/H7bPAc4kpVC56cWbAJh63NQ0Pl6k6lR0z7BC4bHqt8Gfg0+sbn8iqlJTMbMhwLHAIuCIMHBw93Vmdnh4WH9gTeRtHWFbV+0dedrznb+NYETDoEGDyrqG1ze/XvwgkTqnonsd6WpUUUyh8Bh8Iow+G1ovqaxvFUg9VMzsIOBh4OvuvrWLske+F7yM9r0b3WcAMwBaW1vzHiPSDHSTYpV1FRyVjCrqIDwKSTVUzKwHQaDc5+6/CJvXm1nfcJTSF9gQtncAAyNvHwCsDdtP6tT+VNg+IM/xIhKheklK4ow0ugqOOg6GSqQWKuFKrLuB5e5+a+Sl+cDFwPTwz0ci7VPMbB5Bof7PYfA8AdyYWyUGfB74prtvNrNtZnY8wbTaRcCP0roekSxRvaQMpU5HxRlpNGhwdCXNkcoJwIXAUjNbErZdQxAmD5rZV4DVwKTwtceACcBK4D3gEoAwPL4LvBQe951c0R64DJgF7E9QoNfKL2laqpeQTp2ikCYMjDjSXP31W/LXPQBOznO8A5cX+Kx7gHvytLcDoyropkjDaOh6SdywaNA6RZbojnqRDGu4ekmlS2UVDDWnUBHJsOjopO7rJZUUthUWmaFQEcmYuh+dVDLaUHhknkJFJAPqZjWXRhtShEJFJAOqWoSv9IY9hUdTU6iI1KlUprl0w56kTKEiUkcSm+ZSXUNqpGiomFkLcAbwGaAf8FfgVeBX7r4s3e7VlxWbV3DJfwX/oU04cgKTjp5U5B0ipSl5mqvU8FBgSMq6DBUzuw74AsFeW4sI9ulqAY4GpoeBc6W7v5JuN2tvwpETdn2/YvMKAIWKJCLWNJfCQzKi2EjlJXe/rsBrt4bb1pe3l3zGTDp60q4QyY1WRMrR1eN695jmigaJwkMyostQcfdfFXl9A7t3GRaRQiIBsXXuWrZv2EHL4fsCcMDAFg4ZcRA9x2wKjn1vNsycvWeQKDwkI+LUVAYA5wMn0qmmAjzu7jtT7aFIVkWCZMvCl9m6an9o+diuQBl8Qb+u368gkQwqVlOZSfA0xUeBm9izpjIe+JaZTXP3p9PuqEjdilHv2Lq+D9v/4rQMGU1LH4IprkbZ7FEkothI5RZ3fzVP+6vAL8xsX5qkpiJSarF8y5ZPBqOT/+7N9nc20jKyzrZTEUlBsZrKqwBmdoW73x59LdK2MsX+iVRfQiuttl54Edvffp2W4X3qf7NHkYTEvfnxYuD2Tm2T87SJZFNCK63qfrNHkZQVq6mcD1wADDWz+ZGXDgY2pdkxkcSUujVJiQXyutnsUaQOFBupPAesAw4Dbom0bwMa/oZHybA4I4+oClZaNfQTF0VKVCxUVrv7KuDvCh1gZhY+ClikuuLuppvC0lxNc4nkVyxUnjSzh4FH3H11rjFc9XUiQa3lSWBWaj2sU9oHrIrK2Rwx5Xs8MvXERZEqKhYq44EvA/eb2VDgHWB/oBvwa+AH7r4k3S7WH+0DlpI6399KoxOR4ootKd4O3AHcYWY9CGorf3X3d6rRuXqlfcDKkNEnBqoIL1Ka2M9TcfcPgHVmdqCZfQm4wN1PT69rknlVLJanRUV4kdLECpWwhjKBYHnxeOBh4Gcp9kuyqlCQ1GFgFKJpLinXBx98QEdHB9u3b691V8rW0tLCgAED6NGjR1nvL3afyikEm0meSlCQvxc4zt3r/28GqZ4GCJIoFeGlXB0dHRx88MEMGTIEM6t1d0rm7mzatImOjg6GDh1a1mcUG6k8ATwDnOjubwKYme6il4YLEo1OJAnbt2/PbKAAmBm9e/dm48aNZX9GsVAZC5wH/F8zewOYB3Qv+2ySXZ0L7Q0QJFEanUhSshooOZX2v9jqr98BvwOmmtkJBFNh+5rZ48Av3X1GRWeX+tZVob0BgkSjE2lEa9as4bOf/SyLFy+mV69ebNmyhU996lM89dRTXHrppbzwwguceOKJPBp5+miSusU90N2fdfcpBM9XuY0u7rKXBrH05/DHpcH3g0+EM26DS361+yvDgQK7RyeARifSMAYOHMhll13GtGnTAJg2bRptbW0MHjyYq6++mnvvvTfV8xcr1A9x97eibeGTHp8AnrBgnNTf3TvS66JUVXR08sel0Gd0ECANQqMTaQbf+MY3GDt2LLfddhu//e1v+dGPfgTAySefzFNPPZXquYvVVG42s27AI8BiYCPBkx//BjgJ+BxwLbBXqJjZPcAZwAZ3HxW2XQf8U/g5ANe4+2Pha98EvgJ8BHzN3Z8I28cTbLHfHfjf7j49bB9KUOPpBbwMXOjuO0r+XyABmd+ypdA0V5/RwRRXA1HtRKrl3/5zGa+t3ZroZ47odwjXfmFk0eN69OjBzTffzPjx4/n1r3/Nvvvum2g/ulKspjLJzEYAXyLYrqUv8B6wHHgMuDG86z6fWcCPgc7/DPyBu/97tCE8x3nASKAfwcKAo8OXfwKcQhBcL5nZfHd/jeDxxj9w93lm9jOCQPpp8UtOVkNs2ZKb5uozuiFqJZ1pdCLN6PHHH6dv3768+uqrnHLKKVU7b9GbH8O/wL9V6ge7+9NmNiTm4ROBee7+PvCmma0EjgtfW+nubwCY2TxgopktB/6B4GZMgNnAddQgVDK7ZUuDT3NFaXQitRBnRJGWJUuWsGDBgl1F+fPOO4++fftW5dyxt2lJ0BQzuwhoB6509y0Exf8XIsd0hG0Aazq1jwN6A++4+4d5jpc4oqOTBpzm0uhEmpW7c9lll3HbbbcxaNAgrr76aq666iruu+++qpw/9uqvhPwUOAoYQ/Dwr9yDv/ItjPYy2vMyszYzazez9kpu6sm89pkw8/TgKzo6aYCVXJ1pZZc0q7vuuotBgwbtmvL66le/yuuvv85vfvMbPvOZzzBp0iQWLlzIgAEDeOKJJxI/f1VHKu6+Pve9md0F5BZKdwADI4cOANaG3+dr/xNwqJntE45WosfnO+8MYAZAa2tr8z5QrIFHJ9GRCWh0Is2rra2Ntra2XT93796dxYsXA/DMM8+kfv7YoWJm/YHB0fe4+9OlnMzM+rr7uvDHs4BXw+/nA3PN7FaCQv0w4EWCEcmwcKXX2wTF/Avc3c3sSeBsghVgFxOsUJPOGrh2UmhbetDoRKRW4u5SfBNwLvAawZJfCKabCoaKmd1PsOz4MDPrIFh6fJKZjQnf+xZwKYC7LzOzB8PP/xC43N0/Cj9nCsF9Md2Be9x9WXiKqcA8M7ue4K7/u+NdcpNp4NGJtqUXqT9xRypnAp8IV2fF4u7n52ku+Be/u98A3JCn/TGC5cud299g9woxiWrg0UlnmuISqS9xQ+UNoAcQO1SaVV3cCNnAo5N8q7pEpH7EDZX3gCVmtpBIsLj711LpVUbV9EbIBh6d6JG+ItkRN1Tmh1/ShZreCNnAoxPVTkSyI1aouPvs8JHCua1TVoTPrJcqm7toNY8seXvXzxPH9OeCcYOCHxpodNKZaici8RTa+n7WrFlMmzaNrVu30r17d771rW9x7rnnJn7+uKu/TiLYCuUtgmW+A83s4lKXFEvXOgdGPove3AzAuKG9OHL1Qxz59nMsW7gPQz54g3d7HsMR1ehoFah2IlKe6Nb3M2bM2LX1fd++fZkzZw7Dhg1j7dq1jB07llNPPZVDDz000fPHnf66Bfi8u68ACDd7vJ/gyZBSps4hEg2MQsYN7bVrdLL+h9/moC2reYsjWfrRIB5Z/7e8cefze71nj9FMRmi/LpHy5dv6PrpTcb9+/Tj88MPZuHFjzUKlRy5QANz992bWI9GeNKFHlrzNa+u2MqLvIcCegRHHEQe3wMHHMvKSXzF30WreyDPKWfTmZha9uXlXeNVzwGi/Lmkoj0/b/ZC7pPQZDadNL3pYsa3vX3zxRXbs2MFRRx2VbP+IHyrtZnY3kHtk2JcInq8iJYqOTnKB8sClJTxEM98qL+CCcYPyhkX0fPUeMBqdiCSn0Nb369at48ILL2T27Nl065b89o9xQ+Uy4HLgawQ1laeBOxLvTYPJ3bOyYev7fLD1b+n50Wf3mOIa0fcQJo4pcXPlEld5RcOmHgNGoxNpWDFGFGkptPX91q1bOf3007n++us5/vjjUzl33NVf7wO3hl8SQ/SelY73VrLzw7/Q0z5b8hRXXmWu8qrHgNHoRCRZhba+nzlzJmeddRYXXXQRkyald/9csWfUP+ju55jZUvJsLe/un0ytZxn3wZZxvLdqAAA7fToH7LcPD1xSwjRXVIEpr0rUMmA0OhFJT76t72fNmsX3vvc9nn76aTZt2sSsWbMAmDVrFmPGjEn0/OZeeCf43K7CZjY43+vuvirR3lRBa2urt7e3p36ec+98flfN5K19/53DDtyPX50zt7wPiz7/BFJ93G/ngIHdq9GSCphVF160xzJh3cwojWL58uUcc8wxte5GxfJdh5ktdvfWYu8t9oz63Db1X3X3qZ1OcBPBTsFSQK4If8l/HVL5h1XpxsZCI5jX1m3d9Xo5NDoRaQ5xC/WnsHeAnJanranlW9mVZdGAyY28zo3cB1PKyEW1E5HmUKymchnwVeBIM3sl8tLBwLNpdiyLovedlLWyq451vpY4tReNTkSaT7GRylzgceB7wLRI+zZ335xarzKs5PtOCkmhOF+JzvfBxCnua3Qi0nyK1VT+DPwZOB/AzA4HWoCDzOwgd1+dfhcbQ8nPWanzXYfj1l40OhFpLnE3lPwCwT0q/YANBM+qXw6MTK9r2RCnjlL2c1YysutwNGC+f/n3OHLhszz+0D702biaD4YOI+/SQRFpSHHv0b8eOB74vbsPBU5GNRVgdx0FKFhHmXT0JGaOn8nM8TP5RK9PVLuLVXXS20s4autaAH5/UF9mthzNuXc+z7l3Ps/cRRrYiqRtzZo1DB06lM2bgwrFli1bGDp0KKtWrWLs2LGMGTOGkSNH8rOf/SyV88dd/fWBu28ys25m1s3dnwyXFAuNW0eJK1qQ/9jbb9LyyVEce+8c5i5azeaEliSLSDxdbX3/3HPPsd9++/Huu+8yatQovvjFL9KvX79Ezx83VN4xs4MI9vy6z8w2AB8m2hOp+zpKIYUK8l0tSa63zSxFGkmxre/ff/99du7cmcq544bKROCvwDcIdij+GPCdVHqUAanej5KROkpnxQry0WnBetnMUiRNN714E69vfj3RzxzeazhTjyt+e2Chre/XrFnD6aefzsqVK7n55psTH6VAjFAxs+7AI+7+OWAnwRMgm1ql96OUvBKsDpX6ZMa07tQXkfzybX0/cOBAXnnlFdauXcuZZ57J2WefzRFHJPu82KKh4u4fmdl7ZvaxcImxUH4dpeyVYHWmkntQNC0mzSDOiCIthba+z+nXrx8jR47kmWee4eyzk51mjzv9tR1YamYLgL/kGt39a4n2pglMOnrSrhC55KHTghrKzNODF+u4OB8dmUByd8hHR3katYhUrtDW9zfddBO9e/dm//33Z8uWLTz77LP8y7/8S+Lnjxsqvwq/ogpvbyzx/GUj7PgL5B7MXMfF+ejIBEjsDnmNWkSSVWjr+7vvvpuHH34YM8Pdueqqqxg9Ovl/xMYNlUPd/fZog5ldkXhvmtG+B8Lk+izMV3vvLhXzRSrX1tZGW1vbrp+7d+/O4sXB09+vvfba1M8fN1QuBm7v1DY5T1vDSmvF1wp27CraQ30V7qu9d5eK+SLZV2yX4vOBC4ChZjY/8tLBwKY0O1Zv0tiBeIIfCLb753os3Ndq7y5Ni4lkU7GRynPAOuAw4JZI+zbglbzvaGCJ3DkfuWt+0h/fZFKf0TB+JsAeI5ZaKXWpcDVoWkwkO4rtUrwKWAUksAeJAHV/13w9blevaTGR7Ii7S/H/BG4CDieYsDHA3b1gYcHM7gHOADa4+6iwrRfwADAEeAs4x923mJkR1GcmAO8Bk9395fA9FwPfDj/2enefHbaPBWYB+wOPAVe4ezZWpNX5XfP1vF29psVE6lvcQv33gS+4+/ISPnsW8GMg+rfTNGChu083s2nhz1MJHk08LPwaB/wUGBeG0LVAK8ES5sVmNt/dt4THtAEvEITKeIIHimVaLe62r8cprzi6mhbLva6QEamuuFvfry8xUHD3p4HOT4ecyO5tXmYDZ0ba53jgBeBQM+sLnAoscPfNYZAsAMaHrx3i7s+Ho5M5kc/KrAlHTti1Nf6KzSt47I3HqnLe3JQXJHf/STVcMG4QD1z6dzxw6d9x41mjGTe0167XXlu3dY+AEWkWXW19D7B161b69+/PlClTUjl/3JFKu5k9APwH8H6u0d1/UeL5jnD3deF714VPkgToD6yJHNcRtnXV3pGnPS8zayMY1TBoUHn/ch3RL8FNIwvY4277lIv2jfb8+M6PO9bUmDSrQlvfDx4cPC7vX//1X/n7v//71M4fN1QOIah1fD7S5kCpoVKI5WnzMtrzcvcZwAyA1tbWsuou136hgodc1uFzUuqxIJ8krRiTZpZv63uAxYsXs379esaPH097e3sq544VKu6e1D+b15tZ33CU0pfg0cQQjDQGRo4bAKwN20/q1P5U2D4gz/H1qcwVX2nXV7I+OulKoRVjChiplj/eeCPvL0926/v9jhlOn2uuKXpcvq3vd+7cyZVXXsm9997LwoULE+1XVNzVX0cTFMaPcPdRZvZJ4Ivufn2J55tPcHf+9PDPRyLtU8xsHkGh/s9h8DwB3GhmPcPjPg980903m9k2MzseWARcBPyoxL5UV4krvtLYzTirBflKaUmyNKPOW9/fcccdTJgwgYEDBxZ/cwXiTn/dBVwN3Ang7q+Y2VyCZ9fnZWb3E4wyDjOzDoJVXNOBB83sK8BqIPe35GMEy4lXEkyzXRKeZ7OZfRd4KTzuO+6eK/5fxu4lxY/TACu/otKorzT6lFccWpIs1RJnRJGWfFvfP//88zzzzDPccccdvPvuu+zYsYODDjqI6dOnJ3ruuKFygLu/GNxOskuXjxN29/MLvHRynmMduLzA59wD3JOnvR0Y1VUfZG+NPOVVKtVdpBEV2vr+vvvu23XMrFmzaG9vTzxQIH6o/MnMjiIshpvZ2QTbt0iVlFtfadYprzhUd5FGVGjr+9/85jeprvrKiRsqlxOsnhpuZm8DbwL/mFqvZA+V1Fc05RWP6i7SKLra+j5n8uTJTJ48OZXzx1399QbwOTM7EOjm7ttS6U0jSXAZcaX1FU15lUZ1F5HyxV39dSPwfXd/J/y5J3Clu3+763c2sRpuHKkpr+So7iJSmrjTX6e5+66lDOEmkBPYvdGj5JPSxpHF6iua8kqO6i4ipYkbKt3NbD93fx/AzPYH9kuvW1JI3PqKprySp4CRONydTitlM6XSzd7jhsr/ARaa2UyCFWBfZvfGkFJFheormvKqLgWM5NPS0sKmTZvo3bt3JoPF3dm0aRMtLS1lf0bcQv33zWwpwT0mBnzX3Z8o+6ySmNxU2HlzltH/jx9w8IjRmvKqMgWM5AwYMICOjg42btxY666UraWlhQEDBhQ/sIC4IxXcveHuWs+66FTYex/+lbf77M+pmvKqKS1Nbm49evRg6NChte5GTaX25EdJ3+d+5xz36EcAbF7vbCi4+b/UgpYmSzNK88mPkrLoKq8N/Q/gyeEfMa/KT42UePSUSmkWcUOl5Cc/NqUaPDclt8rrxd8/xPrwSZFJ7WosySk0LQaaGpPGYnGWj5nZ7UAfKn/yY821trZ6Wg+nYebpe4bJ6LOhNb0nOK668CKAvZYOX/Jfl7Bi84pdjybWqKW+5abGRvQNZpM1apF6ZGaL3b212HH18uTHxpHSDY85cZYOp/EsFkmP7tqXRhJrpNJIUh+pQKqhsurCi/YIk0POOIOe555T8HiNWrKl85JkgHFDe+11nMJGqi3RkYqZDSB4suIJBCOU3wJXuHtHRb2UspRyt7xGLdnSVe0lR6MZqWdxayoLgLnAvWHTPwJfcvdTUuxbKrI4Usk35VXOFiwatTSGrkYzChhJS9I1lY+7+8zIz7PM7OvldU1KldQGkRq1NAbdwS/1rJQnP/4jcH/48/nApnS6JPkksUFk533Dyn2apNQPBYzUm7ih8mXgx8APCGoqz4VtklEatTQeBYzUA63+qlS+Gx4TqKkkVUeJQ7WWxqYajCQh6dVfswlWe0Wf/HiLu2u0ktITHqv5oK3oqKV9fTvt69t5LLw7XwGTfRrBSDXFXf31O3c/tlhbFiQ+UklpxVehu+XT9tDvH9oVKLnRy8zxM4u8S7JIIxgpRdKrv7qZWU933xJ+eK8S3isZomJ+84g7ggGFjMQXNxhuAZ4zs58TFOrPAW5IrVdNqt6e3tjVtFjudYVMY+jqpktNk0kpYhfqzWwE8A8Ez1JZ6O6vpdmxtNTz9FepW7BUU3RaDIKQAWg9Yu/RsMKmsWiaTCD56S/CEMlkkGRJmqu8KhGdFoO9QyZHy5Mbjwr9UgrVRaQsnUMmR3WYxlZKHSZHYdNcFCo1FK2hQH3UUSql5cnNQ5tfSj66+bFSFdRUOtdQoL7qKJWKTpF1VYOJUvA0FtVjGkfcmopCpVIVhgpU/16UWihUg4lS8b+x6Vkx2VbXoWJmbwHbgI+AD929Nbz35QFgCPAWcI67bzEzA24HJhA8fXKyu78cfs7FwLfDj73e3WcXO3cioZLQ1izNFCpxFAoehU3j6Wq6DDSaqUdZCJVWd/9TpO37wGZ3n25m04Ce7j7VzCYA/0wQKuOA2919XBhC7UArwb0zi4GxuRs0C0kkVCp4Fn019/RqFOWEDShwskijmfqVxVBZAZzk7uvMrC/wlLt/wszuDL+/P3pc7svdLw3b9ziukMRCBRKpozRSDaXauppS6xw4Cpjs0WimviR+n0rCHPi1mTlwp7vPAI5w93UAYbAcHh7bH1gTeW9H2FaofS9m1ga0AQwaVPv/82l0koxCy5ph70UCnXcDyEfBU1+iq8uitKVMfatVqJzg7mvD4FhgZq93cazlafMu2vduDEJrBgQjlVI7K9kTDZy4iwS0/DkbtKVMfatJqLj72vDPDWb2S+A4YL2Z9Y1Mf20ID+8ABkbePgBYG7af1Kn9qZS7LhnU1YgmJ+7IRmFTXzqPZnRDZu1VvaZiZgcC3dx9W/j9AuA7wG11AUcAAAfESURBVMnApkihvpe7/y8zOx2Ywu5C/Q/d/biwUL8Y+FT40S8TFOo3d3X+WtRUVJzPFq1CawxxazJRCpvC6rZQb2ZHAr8Mf9wHmOvuN5hZb+BBYBCwGpjk7pvDJcU/BsYTLCm+xN3bw8/6MnBN+Fk3uHvRB3/UIlRUnG8MCpvGoAUA5anbUKm1WoUK6H6URhU3bBQw9U3LmbtW76u/RBpGoZpNOSvQOlMQVU85e5lFNWvYdKZQEUlJqSvQOtOKtNqJs5w5SqvOdtP0VzliTH+pOC+Virshp8Km9uJOnUF2A0c1lQLKDpUS9/tScV6SpEUC2VFoNAPZrtUoVAooO1RK3O9LxXmphnL3RctR8FRXlpc5K1QKqChUQCu+JBP0qIFsyULYaPWXSBMrdReBqGIr1RQ4yat0YUBUrUc2GqnEpeK8NIlSdn8uRMGTrnJGNiP6HcK1XxhZ9jk1UqmBrY8+uitMWoYP55Azzqh1l0RKFnf350K0d1r6Sh3ZVJNGKnHFGKmojiKixQONSiMVEamJODsMFKJRTvYpVESkKtJcPFCIgqj6FCoVylecF5HyVDLK6aycIFIIVU6hUiEV50XSF2eU01mpQaQQSoZCJQFaOixSf0oNorRCqNmCR6EiIkI6IVTK6KdRwkehIiJShkoWHnTWSFNvChURkZTEHf000tSbQiWu3O7EIiIJq8bU2/Bew5l63NSK+hmHQiWu06bv+lbLiEWklpKcekuaQqUMWkYsIvWunGXYSVColEnLiEVE9tat1h0QEZHGoVAREZHEKFRERCQxChUREUmMCvUxaRmxiEhxGqnElFtGDGgZsYhIARqplEDLiEVEuqaRioiIJEahIiIiicl8qJjZeDNbYWYrzWxarfsjItLMMh0qZtYd+AlwGjACON/MRtS2VyIizSvToQIcB6x09zfcfQcwD5hY4z6JiDStrK/+6g+sifzcAYzrfJCZtQFtAIMGDSrrRPsdo/tSRESKyXqoWJ4236vBfQYwA6C1tXWv1+Poc8015bxNRKSpZH36qwMYGPl5ALC2Rn0REWl6WQ+Vl4BhZjbUzPYFzgPm17hPIiJNK9PTX+7+oZlNAZ4AugP3uPuyGndLRKRpZTpUANz9MaD6D2IWEZG9ZH36S0RE6ohCRUREEqNQERGRxChUREQkMeZe1r2AmWVmG4FVZb79MOBPCXYnS5r52qG5r7+Zrx2a+/qj1z7Y3T9e7A1NFyqVMLN2d2+tdT9qoZmvHZr7+pv52qG5r7+ca9f0l4iIJEahIiIiiVGolGZGrTtQQ8187dDc19/M1w7Nff0lX7tqKiIikhiNVEREJDEKFRERSYxCJQYzG29mK8xspZlNq3V/qs3M3jKzpWa2xMzaa92ftJnZPWa2wcxejbT1MrMFZvaH8M+etexjWgpc+3Vm9nb4+19iZhNq2ce0mNlAM3vSzJab2TIzuyJsb/jffRfXXvLvXjWVIsysO/B74BSCh4K9BJzv7q/VtGNVZGZvAa3u3hQ3gJnZZ4F3gTnuPips+z6w2d2nh/+w6OnuU2vZzzQUuPbrgHfd/d9r2be0mVlfoK+7v2xmBwOLgTOByTT4776Laz+HEn/3GqkUdxyw0t3fcPcdwDxgYo37JCly96eBzZ2aJwKzw+9nE/wH13AKXHtTcPd17v5y+P02YDnQnyb43Xdx7SVTqBTXH1gT+bmDMv/HzjAHfm1mi82srdadqZEj3H0dBP8BAofXuD/VNsXMXgmnxxpu+qczMxsCHAssosl+952uHUr83StUirM8bc02Z3iCu38KOA24PJwikebxU+AoYAywDriltt1Jl5kdBDwMfN3dt9a6P9WU59pL/t0rVIrrAAZGfh4ArK1RX2rC3deGf24AfkkwJdhs1ofzzrn55w017k/VuPt6d//I3XcCd9HAv38z60Hwl+p97v6LsLkpfvf5rr2c371CpbiXgGFmNtTM9gXOA+bXuE9VY2YHhoU7zOxA4PPAq12/qyHNBy4Ov78YeKSGfamq3F+oobNo0N+/mRlwN7Dc3W+NvNTwv/tC117O716rv2IIl9HdBnQH7nH3G2rcpaoxsyMJRicA+wBzG/36zex+4CSCbb/XA9cC/wE8CAwCVgOT3L3hCtoFrv0kgukPB94CLs3VGBqJmZ0IPAMsBXaGzdcQ1BYa+nffxbWfT4m/e4WKiIgkRtNfIiKSGIWKiIgkRqEiIiKJUaiIiEhiFCoiIpIYhYpIiszsUDP7avh9PzP7ea37JJImLSkWSVG4j9KjuR1/RRrdPrXugEiDmw4cZWZLgD8Ax7j7KDObTLDbbXdgFMGeSvsCFwLvAxPcfbOZHQX8BPg48B7wT+7+evUvQyQeTX+JpGsa8N/uPga4utNro4ALCPZTugF4z92PBZ4HLgqPmQH8s7uPBa4C7qhKr0XKpJGKSO08GT67YpuZ/Rn4z7B9KfDJcMfY/wE8FGzNBMB+1e+mSHwKFZHaeT/y/c7IzzsJ/tvsBrwTjnJEMkHTXyLp2gYcXM4bw+dZvGlmkyDYSdbM/jbJzokkTaEikiJ33wQ8a2avAjeX8RFfAr5iZv8PWIYeZS11TkuKRUQkMRqpiIhIYhQqIiKSGIWKiIgkRqEiIiKJUaiIiEhiFCoiIpIYhYqIiCTm/wMjmuN6nRDGcwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.step(times, values[:,0], label = 'X1')\n", + "plt.step(times, values[:,1], label = 'X2')\n", + "plt.step(times, values[:,2], label = 'X3')\n", + "plt.step(times, values[:,3], label = 'X4')\n", + "plt.legend()\n", + "plt.xlabel('time')\n", + "plt.ylabel('concentration (A(t))')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can also evaluate this model using tau-leaping for more efficient but approximate simulations" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "values_exact = model.simulate(k, times)\n", + "values_approx = model.simulate(k, times, approx=True, approx_tau=0.0125)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "def plot_output(vs, suffix=''):\n", + " plt.step(times, vs[:,0], label = 'X1'+suffix)\n", + " plt.step(times, vs[:,1], label = 'X2'+suffix)\n", + " plt.step(times, vs[:,2], label = 'X3'+suffix)\n", + " plt.step(times, vs[:,3], label = 'X4'+suffix)\n", + " plt.legend()\n", + " plt.xlabel('time')\n", + " plt.ylabel('Molecule Count')\n", + " plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEGCAYAAACtqQjWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de7xd853/8ddHJE5dMoSSk8tJwiMVETNap9KLR8eMUREqdKjLDEJ/c1LFtJ3ykxr96UU76Zgq1TKiJOGBoNpSgzJ+4tK6nWNSRGjzQ+RIGiRpj0hd8/n9sdZOlpN9WXvvtfbea6/38/HI4+z93Wuv9V09Td6+12XujoiISBK2anYFRESkfShUREQkMQoVERFJjEJFREQSo1AREZHEbN3sCjTaLrvs4uPHj292NUREMqWvr+81d/9gpeNyFyrjx4+nt7e32dUQEckUM1se5zh1f4mISGIUKiIikhiFioiIJEahIiIiiUktVMxsrJndZ2ZLzWyJmX0pLP+Gmb1sZovDP9Mj3/mamS0zs+fM7JBI+bSwbJmZzY6UTzCzR83s92Z2o5kNS+t+RESksjRbKu8CX3X3vYCPAaeb2eTwsx+4+77hnzsAws+OA/YGpgGXmdkQMxsC/Bg4FJgMHB85z/fCc00E1gGfT/F+RESkgtRCxd1XufsT4evXgaXA6DJfmQEsdPe33P0FYBmwf/hnmbs/7+5vAwuBGWZmwN8CPw2/vwA4Mp27ERGROBqyTsXMxgMfBh4FPgmcYWYnAb0ErZl1BIHzSORr/WwOoRWDyqcCOwN/dPd3ixwvIpJr6268iYHbb9/0fpu9JjHy3HNTv27qoWJm2wO3AF929wEzuxz4NuDhz+8DpwJW5OtO8daUlzm+WB16gB6Arq6uam9BRCQTokGy4fHHAdj2ox9taB1SDRUzG0oQKNe5+88A3H115PMrgUKU9gNjI18fA6wMXxcrfw3Y0cy2Dlsr0ePfx93nAnMBuru79VQyEcm0wa2QgmiQbPvRjzL88MPZ6djPNbRuqYVKOOZxFbDU3S+KlHe6+6rw7VHA0+Hr24DrzewiYBQwEXiMoEUy0cwmAC8TDOaf4O5uZvcBRxOMs5wM3JrW/YiINFOcVkizgiQqzZbKJ4ETgafMbHFYdi7B7K19CbqqXgRmAbj7EjO7CXiGYObY6e7+HoCZnQH8ChgCXO3uS8LznQMsNLMLgP8hCDERkbZQKkhaITxKsbw9o767u9u1oaSItKpyLZJmBomZ9bl7d6XjcrdLsYhIK2jlcZF6KFRERBokK+Mi9VCoiIikZHBrJOutkDgUKiIiCSrXGmnXIIlSqIiI1CmLs7TSolAREamBgqQ4hYqISEwKksoUKiIiZShIqqNQEREZREFSO4WKiAgKkqQoVEQktxQkyVOoiEiuKEjSpVARkbanIGkchYqItCUFSXMoVESkbShImk+hIiKZlccNG1udQkVEMiXvGza2OoWKiGTKwO238+azz9IxaZJCpAUpVESk5UVbJ4VAGXftNU2ulRSjUBGRllSqm6tj0iSGH354M6smZShURKRlaPZW9ilURKSpFCTtRaEiIg2nIGlfChURaTjN4GpfChURaQjN4MoHhYqIpEYzuPJHoSIiidJ4Sb4pVESkbgoSKVCoiEjdNPAuBQoVEamJBt6lGIWKiMSmgXepZKu0TmxmY83sPjNbamZLzOxLYfkIM7vHzH4f/twpLDcz+6GZLTOzJ83sI5FznRwe/3szOzlSvp+ZPRV+54dmZmndj4hs7uaCIExGfvObjLv2GsZde426uwRIt6XyLvBVd3/CzHYA+szsHmAmcK+7zzGz2cBs4BzgUGBi+GcqcDkw1cxGAOcD3YCH57nN3deFx/QAjwB3ANOAO1O8J5HcUTeXVCO1UHH3VcCq8PXrZrYUGA3MAA4MD1sALCIIlRnANe7uwCNmtqOZdYbH3uPuawHCYJpmZouA4e7+cFh+DXAkKYXK9x77HgDn7H9OGqcXaRnlnqaobq4W0TsPnvppdd8ZuQ8cOied+kQ0ZEzFzMYDHwYeBXYLAwd3X2Vmu4aHjQZWRL7WH5aVK+8vUl7s+j0ELRq6urpquodn1z5b0/dEskBPU2whcQJj+UPBz3EHpF+fKqUeKma2PXAL8GV3Hygz7FHsA6+hfMtC97nAXIDu7u6ix4jkmaYEJ6iWVkRUnMAYdwDsczR0n1L7dVKSaqiY2VCCQLnO3X8WFq82s86wldIJvBKW9wNjI18fA6wMyw8cVL4oLB9T5HgRiUFjJXUqFR71tiJaODDiSC1UwplYVwFL3f2iyEe3AScDc8Kft0bKzzCzhQQD9X8Kg+dXwHcLs8SATwNfc/e1Zva6mX2MoFvtJODStO5HpB1oSnBM9XRBZTwU6pVmS+WTwInAU2a2OCw7lyBMbjKzzwMvAceEn90BTAeWARuAUwDC8Pg28Hh43LcKg/bAacB84AMEA/Sa+SVShrq5BqmntZHz8CglzdlfD1F83APgoCLHO3B6iXNdDVxdpLwXmFJHNUXanrq5qD48FBg104p6kTbU1t1ctQyEKzwaRqEi0oYy282V1nRahUfDKFRE2kSmurk0ltG2FCoiGdaS3VyaOZVrChWRjGmZB2KptSFFKFREMib18ZK4A+FqbUgRChWRDEhlvKTeFeEKDylCoVKF59Y+xyl3BX+Bpu8+nWM+dEyFb4gkI9o6qXq8RGs0pIEUKjFN3336ptfPrX0OQKEiqaq6daLwkBagUInpmA8dsylECq0VkaSVnM01ekeGj3gB5h1W+ssKD2kBChWRFjJww09484WX6dh1GNuO7WD45O3Zad81sPzR8AjNqJLWplARabRB3VTrFg8w8Mx6AN78wwY6dnyHcSeMf/93FBiSEQoVkSTFmI677t4nGFj+Aej4CwA2rHgTgG3HdtAxcluGH/wpOOWS1KsqkgaFikgt6piOO7B6JG++4XSM3weAbUeSnb25RCpQqIiUk9CMqvfN5Prjq3Ts3cL7conUQaEiMlg0SBKaUVXXOhORDFGoSL5Uu9lhHQPkmdo1WCQhFUPFzO5194MqlYm0lCZtdtiSuwaLNFDJUDGzDmBbYBcz24nNjwYeDoxqQN1EyivX6mjSQsDMPhxLJCHlWiqzgC8TBEgfm0NlAPhxyvVqedoHrEnijHcUyhq0rkPdXCKblQwVd78EuMTMznT3SxtYp5anfcAaIE73VRMXBKqbS6S4imMq7n6pmX0CGB893t1z+59i2gcsQRndBFHdXCLFxRmovxbYA1gMvBcWO5DbUJE6pTBltxHUzSVSWZwpxd3AZHf3tCsjbaxUkLRgeJSitSZSyTvvvEN/fz9vvvlms6tSs46ODsaMGcPQoUNr+n6cUHkaGAmsqukKki8tPhZSLbVOpBr9/f3ssMMOjB8/HjOr/IUW4+6sWbOG/v5+JkyYUNM54oTKLsAzZvYY8Fbk4kfUdEVpPxntzipFg/BSqzfffDOzgQJgZuy88868+uqrNZ8jTqh8o+azS/tqg+6sUjQIL/XIaqAU1Fv/OLO/7q/rCtIeBndrtVmQqJtL2sWKFSv41Kc+RV9fHyNGjGDdunV85CMfYdGiRcyaNYtHHnmEAw44gNvD/78nLc7sr9cJZnsBDAOGAm+4+/BUaiSto1y3VhsESZQG4aVdjB07ltNOO43Zs2czd+5cZs+eTU9PD+PGjePss89mw4YNXHHFFaldP05LZYfoezM7Etg/tRpJc7Vxt9Zgap1Iu/rKV77Cfvvtx8UXX8xDDz3EpZcG69cPOuggFi1alOq1q96l2N1/YWaz06hMVmV+y5YcBUmUWieSpm/+cgnPrBxI9JyTRw3n/M/sXfG4oUOHcuGFFzJt2jTuvvtuhg0blmg9yonT/fXZyNutCNatVFyzYmZXA4cDr7j7lLDsG8A/AYWpBee6+x3hZ18DPk+wwPKf3f1XYfk04BJgCPATd58Tlk8AFgIjgCeAE9397Ur1SlpbbNny1E/hD0/ByH3aPkjUOpG8uPPOO+ns7OTpp5/m4IMPbth147RUPhN5/S7wIjAjxvfmAz9iy5X3P3D3/4gWmNlk4Dhgb4INLP/bzD4Ufvxj4GCgH3jczG5z92eA74XnWmhm/0kQSJfHqFeiMrtlS7R1UgiUU/6ruXVKiaYISzPEaVGkZfHixdxzzz2bBuWPO+44Ojs7G3LtOGMqNf1L6e4PmNn4mIfPABa6+1vAC2a2jM3jNsvc/XkAM1sIzDCzpcDfAieExywgmPrc8FDJrGjrZOQ+QeukTWmKsOSJu3Paaadx8cUX09XVxdlnn81ZZ53Fdddd15Drx+n+GgNcCnySoNvrIeBL7t5f4zXPMLOTgF7gq+6+DhgNPBI5pj8sA1gxqHwqsDPwR3d/t8jxxe6hB+gB6OrqqrHabSCnrRN1c0meXHnllXR1dW3q8vriF7/I/Pnzuf/++znvvPN49tlnWb9+PWPGjOGqq67ikEMOSfT6cbq/5gHXA4WBgn8My2rppLsc+DZBOH0b+D5wKpuf1RLlBGM4xcpLHV+Uu88F5gJ0d3fnaw+zUoPwOWqdqJtL8qSnp4eenp5N74cMGUJfXx8ADz74YOrXjxMqH3T3eZH3883sy7VczN1XF16b2ZVAYfVNPzA2cugYYGX4ulj5a8COZrZ12FqJHi9RORmEj7ZMQK0TkWaJEyqvmdk/AjeE748H1tRyMTPrdPfCxpRHEWxWCXAbcL2ZXUQwUD8ReIygRTIxnOn1MsFg/gnu7mZ2H3A0wQywk4Fba6lTW8pRN1dBtGUCqHUi0iRxQuVUgllcPyDoYvpNWFaWmd0AHEjwjPt+4HzgQDPbNzzPiwSPLMbdl5jZTcAzBDPMTnf398LznAH8imBK8dXuviS8xDnAQjO7APgf4KoY95IPORmE17iJSOuJM/vrJaDqHYnd/fgixSX/4Xf37wDfKVJ+B3BHkfLn0cr+0nLWOlHLRKQ1lAwVM/t34Hl3/89B5V8BRrr7OWlXLouatrq+WJdXG1LrRKS1lWupHA5MKVJ+CfAkQfeTRKS1uv76R1/i1sUvF/1sxr6jOWFqV266vNQ6EWlt5ULF3X1jkcKNlvUHBqSk3tX1pcLj0RfWAjB1woj3le/+0s3s/vJvWHLv1ox/53leHLo733r7PABmvDd608rQrFPrRCS+Ulvfz58/n9mzZzMwMMCQIUP413/9V4499tjEr18uVDaY2UR3/3200MwmAn9OvCY5NDhESoXH1AkjNrdIIlb/8Dy2X/cSL7I7Lw7dnV9/4G82nefRF9ZuOnex72aJWici8ZXa+r6zs5NrrrmGiRMnsnLlSvbbbz8OOeQQdtxxx0SvXy5U/g9wZzi7qi8s6wa+BtS0TkXeHySDQ6RUeJSy2w4dsMOH2TsckN+bYNuAwdfIYsCodSJSu2Jb30d3Kh41ahS77rorr776auNCxd3vDJ+dcjZwZlj8NPD37v5UorXIkVsXv8wzqwaY3Dm86hABYg3InzC1a9M5sxowap1I5t05O/g7mqSR+8ChcyoeVmnr+8cee4y3336bPfbYI9n6UWFKsbs/TbCwUOoQ/Ye9ECg3zvp4bSerckC+VMA8s2pg0+etQq0TkeSU2vp+1apVnHjiiSxYsICttiq2E1Z9qn5Il1Qv2jqZ3DmcGfuW3PsynhrXoEQD5tgrHuaZVQMce8XDWxzXrBaMWifSVmK0KNJSauv7gYEBDjvsMC644AI+9rGPpXJthUpKXhl4i/4Ny5g67+/Z4O/SOeYT3Pi/WmcWdqlga3QXmVonIskqtfX9vHnzOOqoozjppJM45pj01s8pVFLyzsBfsfHdN2Ab2KpjFUO3/W3tJ0thYWO01RLV6C4ytU5EklVq6/t/+7d/44EHHmDNmjXMnz8fgPnz57Pvvvsmen1zL78TfPgExsuB3dx9ipn9JXCEu1+QaE0apLu723t7e1O/TqFb6cZZH9+0ZmXetHnlvlLavMPeHyYN2mm40EU2uXM4kE6rZfmJJwGodSJtYenSpey1117Nrkbdit2HmfW5e3el78ZpqVxJMAPsCgB3f9LMrgcyGSqZ1YS9vKJdZIO7xQqf1xIyxbq8RKQ9xAmVbd39sUGL6N8tdXCeFZvllWWlZo5BfV1j6vISaV9xn6eyB+GTFc3saGBV+a/kU6KzvFpsg8jBYzCDZ49VarVoQF4kH+KEyukEj+KdZGYvAy8QPFJYiqhrDUpUi28QWa5rrFjAqHUikg9xnqfyPPB3ZrYdsJW7v55+tQRo6WeixFlUqdaJSP6Ue57Kv5QoB8DdL0qpTpkSdxylac9ZaYBSiypn3nw9Y9b2M3zvyWqdiOREuZbKDg2rRYbFGUeJ/ZyVFhtHqcWp637LxofuAmCXPyznmb8Yxc8OOA2AGePbZzt+kVZVauv7RYsW8dnPfpb33nuPd955hzPPPJMvfOELiV+/3IaS30z8am2q0jhK7OestPg4ShyTlvyGN/+0ko5Jk1i9/YdYPjpYWNWKe42JtKNyW9//5je/YZtttmH9+vVMmTKFI444glGjRiV6/YpjKmY2j3DmV5S7n5poTSTQwuMopZQaOxkH7B8eU+1sMRGpXaWt79966y02btziGYyJiDP76/bI6w7gKGBlKrXJiHZbj1KvODO7qp0tJpJ133vsezy79tlEzzlpxCTO2b/yHoKltr5fsWIFhx12GMuWLePCCy9MvJUC8WZ/3RJ9b2Y3AP+deE0yJPFdh9tApZld7fCMF5EsKbb1/dixY3nyySdZuXIlRx55JEcffTS77bZbotetZUPJiUDu/9bXsx5l00yw1//A9PUbOIbtgw8yNDhfz1YrWXrGi0it4rQo0lJq6/uCUaNGsffee/Pggw9y9NHJjt1WfEKLmb1uZgOFP8AvgdbZwz1jpu8+nT1H7AnAc2+8zB3vvLL5wwwNzhe6vIC6pgufMLWLG2d9nBtnfZzJncM3jbsce8XDXP/oS0lWWSQXSm1939/fz5///GcA1q1bx69//Wv23HPPxK8fp/tLU4sT9L6ZYPO7YdgwmNn6A/PRlgmks5hR4y4i9Su19f1VV13FLbfcgpnh7px11lnss0/yPSNxZn8dBfxfd/9T+H5H4EB3/0XitZGWFR2Mh/paJ6WoW0ykfj09PfT09Gx6P2TIEPr6+gA4//zzU79+nDGV893954U37v5HMzsfyFWoaMZX5cH4JJV79LFaLSKtK06oFBt3yd0TIxOb8RVdNf/OGzw3bOj7FkS20hYurfLcE3WLiWRHnHDoNbOLgB8TLII8E+hLtVYtKpEdiCOr5qcP3RW223bTR2W3cGmCVtlZWN1iItkRJ1TOBL4O3Bi+vxs4L7Ua5UG4av4YIBofZbdwaZBW31lY3WIirS3O7K83gNlmtr27r497YjO7GjgceMXdp4RlIwjCaTzwIvA5d19nwdbHlwDTgQ3ATHd/IvzOyWwOsQvcfUFYvh8wH/gAcAfwJXffYjsZqU6rtE7iSOtxxyJSuzizvz4B/ATYHugys78CZrn7Fyt8dT7wIyD6n7mzgXvdfY6ZzQ7fnwMcSrCociIwFbgcmBqG0PlAN0HXW5+Z3ebu68JjeoBHCEJlGnBnnJuW92v11kkp5R53rLEXkeaouPgR+AFwCLAGwN1/C3yq0pfc/QFg7aDiGcCC8PUC4MhI+TUeeATY0cw6w+ve4+5rwyC5B5gWfjbc3R8OWyfXRM6VaYXV9qfcdQo3/+7mhlwzqYWMzRRdRHnjrI/z3aP2YeqEEUAw9hINHJF2tmLFCiZMmMDatcE/v+vWrWPChAksX74cgIGBAUaPHs0ZZ5yRyvVjzeJy9xWFh3OF3qvxeru5+6rwnKvMbNewfDSwInJcf1hWrry/SHlrivmclNjPXUlBVloncWnsRfKq1Nb348aNA+DrX/86f/3Xf53a9eOEyoqwC8zNbBjwz8DShOthRcq8hvLiJzfrIegqo6urtn9MJo+qY11KzOekxH7uSgJaZbpwI2hKsuRNsa3vAfr6+li9ejXTpk2jt7c3lWvHCZUvEAyiF1oHdwOn13i91WbWGbZSOoHCxlf9wNjIcWMIttfvBw4cVL4oLB9T5Pii3H0uMBegu7u7psH88z+zdy1f26zFnpOSpQH5emmHZGmGP3z3u7y1NNmt77fZaxIjzz234nHFtr7fuHEjX/3qV7n22mu59957E61XVJzZX68B/5DQ9W4DTgbmhD9vjZSfYWYLCQbq/xQGz6+A75rZTuFxnwa+5u5rw40uPwY8CpwEXJpQHVtG2s+1b7curzi05kXyYvDW95dddhnTp09n7Nixlb9ch5KhYmaXUqZLyd3/udyJw+euHAjsYmb9BLO45gA3mdnngZfYvEzjDoLpxMsIphSfEl5jrZl9G3g8PO5b7l4Y/D+NzVOK76TNZn6lMb6Spy6vODTuImmK06JIS7Gt7x9++GEefPBBLrvsMtavX8/bb7/N9ttvz5w5cxK9drmWSl0dbu5+fImPDipyrFOiS83drwauLlLeC0ypp46tLI3xlTx1eVVL4y7SLkptfX/ddddtOmb+/Pn09vYmHihQJlQKiwwLzGyHoDj+AkhpPXns8opD4y7SLkptfX///fenOuurIM7ixynAtcCI4K29Cpzk7kvSrpzUT11e1VPASJaV2/q+YObMmcycOTOV68eZ/TUX+Bd3vw/AzA4ErgQ+kUqN2kXMtSlx1Tpory6v+mhgX6Q6cUJlu0KgALj7IjPbLsU6tYeYa1PiqHfQXl1eySg3sB+lFozkWZxQed7Mvk7QBQbwj8AL6VWpjSS0NqXaQXt1eaWv1PN01IKRvIsTKqcC3wR+RrCS/QHCKb/SHJW6wtTllb5oqyVKU5PF3Rm0rVWm1LvZe5zFj+sItmaRFhC3K0xdXs2hqcn51tHRwZo1a9h5550zGSzuzpo1a+jo6Kj5HOUWP95W4eJH1HxVqVmprjB1ebUGzRzLtzFjxtDf38+rr77a7KrUrKOjgzFjxlQ+sIRyLZWPE+wQfAPBVijZi90cKHSFHXfNEkb/4R12mLyPurxaRNyAAYVMuxg6dCgTJkxodjWaqlyojAQOBo4HTgD+C7hB61NaR7QrbMO7f+blkR/gEHV5tSQ9UEzyotyK+veAu4C7zGwbgnBZZGbfcve227wxEQmvTank7/7H2f/24NE2a1c7r7TuE2UkYvAgv9a/SDspO1AfhslhBIEyHvghwSwwKSbBtSlxRGd5vTJ6W+6b9B4LU9zVWNKhjS2lnZQbqF9AsGHjncA33f3phtUqyxr83JTCLK/Hfnczq5+/A2j8UyMlOZo9JllnpeYkm9lG4I3wbfQgI9hYso5HITZPd3e3p/XEM+YdFvxMMVSKzfIaPHX4lLtO4bm1z7HniD0BtVqyavDgPsDUCSMABYw0npn1uXt3pePKjalslWyVJAlxFjY281n3khxNT5YsKtlSaVdZb6ksP/EkgNgLG9VqaT9qwUgz1N1SkfagVkv7UQtGWplaKklKqaUSZxwlDrVa2ptaMJImtVTaSFIbREZbLb2re+ld3csd4YwxBUz2qQUjrUChkhFJbBAZ3Tfs5t/dvClQ1C3WfhQw0izq/qpXsVX0CXd/VTs4Xy11i+WHusikVur+apSUVtE3ctdhdYvlh1owkja1VOqV0uD88hNPel+YDD/8cHY69nOJXqOYaLdY7+rgf6fu3YL/OFHAtK9yLRhQyEj8lopCpV4phgqk1+UVR7mAAYVMuyq2izKomyzv1P0ldSs1sA8a3G9n5XZRVjeZVKKWSr0SbKkktR6lEQYP7kepBdO+NNCfX2qpZFBS61EaITq4H6WB/vamgX6pRC2VeiXYUmmFcZR6VRqHKVDYtBe1YNqfWirSFOXGYQrUmmk/asFIgVoq9aqjpRIdQ4HWH0dJiloz+VFpqnKBwqb1aUpxCa0UKoPXokDj1qO0inKtGdAamXYyeKpygbrLsqGlQ8XMXgReB94D3nX3bjMbAdwIjAdeBD7n7uvMzIBLgOnABmCmuz8Rnudk4LzwtBe4+4JK104kVBLamqUdxlDSotZMfmg8JhuyECrd7v5apOzfgbXuPsfMZgM7ufs5ZjYdOJMgVKYCl7j71DCEeoFugscd9wH7ufu6ctdOJFTmHbY5TCDYmqX7lKpPo1CJR62Z/FDAtK4shspzwIHuvsrMOoFF7r6nmV0Rvr4helzhj7vPCsvfd1wpiYUK1D2OkpcxlLTEbc2AAieLNB7TWlp99pcDd5uZA1e4+1xgN3dfBRAGy67hsaOBFZHv9odlpcq3YGY9QA9AV1dz/8+XpbUorS7OTDPYcrZZlMKmdZWaURY1eHYZKGSarVmh8kl3XxkGxz1m9myZY61ImZcp37IwCK25ELRUqq1s0tQ6SV40YAaLO7W5FAVP8w3eOqag2D5lmsLcXE0JFXdfGf58xcx+DuwPrDazzkj31yvh4f3A2MjXxwArw/IDB5UvSrnqkkGlAqdc66ZArZzWpn3KWk/Dx1TMbDtgK3d/PXx9D/At4CBgTWSgfoS7/28zOww4g80D9T909/3Dgfo+4CPhqZ8gGKhfW+76zRhT0ThKdmmSQHZpTCZZLTtQb2a7Az8P324NXO/u3zGznYGbgC7gJeAYd18bTin+ETCNYErxKe7eG57rVODc8Fzfcfd5la7fjFBp1rNRJD2a8pwtcdfIRCls3q9lQ6XZmhUqoKnD7Spua6YchU9zaEFmfK0++0ukbdQzZgOaMNBMcSYAFJthVqCw2ZJCRSQl5WakRdUyYUABk65apzMX5Dls1P1VixjdXxqclyRVs9CzGIVQ8moZp4HsBo7GVEqoOVSq3O9Lg/OSlrjdagWaqdZYpcIGsj1Wo1ApoeZQqXK/Lw3OS6vQTLXWkeVpzgqVEuoKFdCML8k0zVRrHVmb5qzZXyKyBc1Uax1xt54pyMrEALVU4tLgvEjsmWpQudWj4KlOs1s26v4qIc1Q0eC8SG3Bo4CpXdywmTxqOOd/Zu+ar6PuryZR60TyLs76nMGTB9SlVrtqu9HSplARkYaL+yycgrjBM1ieg6hU2KRNoSIiTVVtyyYuPbagORQqItLy4m55E1XLw2iQQZUAAAYZSURBVNkUNvVTqNSp2IwvEWm+aqdPa2wnGQqVOumZ8yLZktaTQKPyHDwKlQRoxpdI9iU1tlPNpIJ2DB+FiohITElOKqhlRlsWQkihIiKSoCSfoxOVlRBSqMRV2J1YRCQB1c5oqzeEJo2YxDn7n1N1PaulUInr0DmbXmrGl4g0WtohlBSFSg0040tEWl0ta3uSoFCpkWZ8iYhsaatmV0BERNqHQkVERBKjUBERkcQoVEREJDEaqI9J04hFRCpTSyWmwjRiQNOIRURKUEulCppGLCJSnloqIiKSGIWKiIgkJvOhYmbTzOw5M1tmZrObXR8RkTzLdKiY2RDgx8ChwGTgeDOb3NxaiYjkV6ZDBdgfWObuz7v728BCYEaT6yQikltZn/01GlgRed8PTB18kJn1AD0AXV1dNV1om720LkVEpJKsh4oVKfMtCtznAnMBuru7t/g8jpHnnlvL10REciXr3V/9wNjI+zHAyibVRUQk97IeKo8DE81sgpkNA44DbmtynUREcivT3V/u/q6ZnQH8ChgCXO3uS5pcLRGR3Mp0qAC4+x1A4x/ELCIiW8h695eIiLQQhYqIiCRGoSIiIolRqIiISGLMvaa1gJllZq8Cy2v8+i7AawlWJ0vyfO+Q7/vP871Dvu8/eu/j3P2Dlb6Qu1Cph5n1unt3s+vRDHm+d8j3/ef53iHf91/Lvav7S0REEqNQERGRxChUqjO32RVoojzfO+T7/vN875Dv+6/63jWmIiIiiVFLRUREEqNQERGRxChUYjCzaWb2nJktM7PZza5Po5nZi2b2lJktNrPeZtcnbWZ2tZm9YmZPR8pGmNk9Zvb78OdOzaxjWkrc+zfM7OXw97/YzKY3s45pMbOxZnafmS01syVm9qWwvO1/92XuverfvcZUKjCzIcDvgIMJHgr2OHC8uz/T1Io1kJm9CHS7ey4WgJnZp4D1wDXuPiUs+3dgrbvPCf/DYid3P6eZ9UxDiXv/BrDe3f+jmXVLm5l1Ap3u/oSZ7QD0AUcCM2nz332Ze/8cVf7u1VKpbH9gmbs/7+5vAwuBGU2uk6TI3R8A1g4qngEsCF8vIPgL13ZK3HsuuPsqd38ifP06sBQYTQ5+92XuvWoKlcpGAysi7/up8X/sDHPgbjPrM7OeZlemSXZz91UQ/AUEdm1yfRrtDDN7Muwea7vun8HMbDzwYeBRcva7H3TvUOXvXqFSmRUpy1uf4Sfd/SPAocDpYReJ5MflwB7AvsAq4PvNrU66zGx74Bbgy+4+0Oz6NFKRe6/6d69QqawfGBt5PwZY2aS6NIW7rwx/vgL8nKBLMG9Wh/3Ohf7nV5pcn4Zx99Xu/p67bwSupI1//2Y2lOAf1evc/WdhcS5+98XuvZbfvUKlsseBiWY2wcyGAccBtzW5Tg1jZtuFA3eY2XbAp4Gny3+rLd0GnBy+Phm4tYl1aajCP6iho2jT37+ZGXAVsNTdL4p81Pa/+1L3XsvvXrO/Ygin0V0MDAGudvfvNLlKDWNmuxO0TgC2Bq5v9/s3sxuAAwm2/V4NnA/8ArgJ6AJeAo5x97Yb0C5x7wcSdH848CIwqzDG0E7M7ADgQeApYGNYfC7B2EJb/+7L3PvxVPm7V6iIiEhi1P0lIiKJUaiIiEhiFCoiIpIYhYqIiCRGoSIiIolRqIikyMx2NLMvhq9HmdlPm10nkTRpSrFIisJ9lG4v7Pgr0u62bnYFRNrcHGAPM1sM/B7Yy92nmNlMgt1uhwBTCPZUGgacCLwFTHf3tWa2B/Bj4IPABuCf3P3Zxt+GSDzq/hJJ12zg/7n7vsDZgz6bApxAsJ/Sd4AN7v5h4GHgpPCYucCZ7r4fcBZwWUNqLVIjtVREmue+8NkVr5vZn4BfhuVPAX8Z7hj7CeDmYGsmALZpfDVF4lOoiDTPW5HXGyPvNxL83dwK+GPYyhHJBHV/iaTrdWCHWr4YPs/iBTM7BoKdZM3sr5KsnEjSFCoiKXL3NcCvzexp4MIaTvEPwOfN7LfAEvQoa2lxmlIsIiKJUUtFREQSo1AREZHEKFRERCQxChUREUmMQkVERBKjUBERkcQoVEREJDH/H9utT3p1QfKvAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plot_output(values_exact)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def mc_estimate_approx(n,tau):\n", + " return sum([model.simulate(k, times, approx=True, approx_tau=tau) for i in range(n)])/n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "def mse(r1,r2):\n", + " return np.square(r1 - r2).mean()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "ename": "KeyboardInterrupt", + "evalue": "", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_line_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'timeit'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'model.simulate(k, times)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;32m/opt/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py\u001b[0m in \u001b[0;36mrun_line_magic\u001b[0;34m(self, magic_name, line, _stack_depth)\u001b[0m\n\u001b[1;32m 2312\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'local_ns'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getframe\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf_locals\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2313\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2314\u001b[0;31m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2315\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2316\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m\u001b[0m in \u001b[0;36mtimeit\u001b[0;34m(self, line, cell, local_ns)\u001b[0m\n", + "\u001b[0;32m/opt/anaconda3/lib/python3.7/site-packages/IPython/core/magic.py\u001b[0m in \u001b[0;36m\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m 185\u001b[0m \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 186\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 187\u001b[0;31m \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 188\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 189\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/opt/anaconda3/lib/python3.7/site-packages/IPython/core/magics/execution.py\u001b[0m in \u001b[0;36mtimeit\u001b[0;34m(self, line, cell, local_ns)\u001b[0m\n\u001b[1;32m 1160\u001b[0m \u001b[0;32mbreak\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1161\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1162\u001b[0;31m \u001b[0mall_runs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtimer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrepeat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepeat\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnumber\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1163\u001b[0m \u001b[0mbest\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mall_runs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mnumber\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1164\u001b[0m \u001b[0mworst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mall_runs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mnumber\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/opt/anaconda3/lib/python3.7/timeit.py\u001b[0m in \u001b[0;36mrepeat\u001b[0;34m(self, repeat, number)\u001b[0m\n\u001b[1;32m 202\u001b[0m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 203\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrepeat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 204\u001b[0;31m \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtimeit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnumber\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 205\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 206\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/opt/anaconda3/lib/python3.7/site-packages/IPython/core/magics/execution.py\u001b[0m in \u001b[0;36mtimeit\u001b[0;34m(self, number)\u001b[0m\n\u001b[1;32m 167\u001b[0m \u001b[0mgc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdisable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 168\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 169\u001b[0;31m \u001b[0mtiming\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minner\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtimer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 170\u001b[0m \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 171\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mgcold\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m\u001b[0m in \u001b[0;36minner\u001b[0;34m(_it, _timer)\u001b[0m\n", + "\u001b[0;32m/opt/anaconda3/lib/python3.7/site-packages/pints/toy/stochastic/_markov_jump_model.py\u001b[0m in \u001b[0;36msimulate\u001b[0;34m(self, parameters, times, approx, approx_tau)\u001b[0m\n\u001b[1;32m 175\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mapprox\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 176\u001b[0m \u001b[0;31m# run Gillespie\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 177\u001b[0;31m \u001b[0mtime\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmol_count\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msimulate_raw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparameters\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 178\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 179\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;32mnot\u001b[0m \u001b[0mapprox_tau\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mapprox_tau\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/opt/anaconda3/lib/python3.7/site-packages/pints/toy/stochastic/_markov_jump_model.py\u001b[0m in \u001b[0;36msimulate_raw\u001b[0;34m(self, rates, max_time)\u001b[0m\n\u001b[1;32m 107\u001b[0m \u001b[0mtime\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 108\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0ma_0\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m0\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0mmax_time\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 109\u001b[0;31m \u001b[0mr_1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr_2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0muniform\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 110\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr_1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0ma_0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 111\u001b[0m \u001b[0ms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mKeyboardInterrupt\u001b[0m: " + ] + } + ], + "source": [ + "%timeit model.simulate(k, times)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%timeit model.simulate(k, times, approx=True, approx_tau=0.0125)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "empirical_mean = sum([model.simulate(k, times) for i in range(30)])/30" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_output(empirical_mean)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "exact_mse = 0\n", + "for i in range(25):\n", + " exact = model.simulate(k, times)\n", + " exact_mse += mse(exact, empirical_mean)/25" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "taus = [0.0125, 0.025, 0.05, 0.1, 0.25, 0.5, 1]\n", + "approx_mses = []\n", + "for tau in taus:\n", + " amse=0\n", + " print(\"Running for tau = \" + str(tau))\n", + " for i in range(1000):\n", + " sim = model.simulate(k, times, approx=True, approx_tau=tau)\n", + " amse += mse(empirical_mean, sim)/1000\n", + " approx_mses.append(amse)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot([0]+taus[:4], [exact_mse]+approx_mses[:4], 'bo')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.4" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/toy/model-schlogl.ipynb b/examples/toy/model-schlogl.ipynb new file mode 100644 index 000000000..7d86c1baa --- /dev/null +++ b/examples/toy/model-schlogl.ipynb @@ -0,0 +1,1168 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Schlögl" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- A + 2X -> 3X, rate k1\n", + "- 3X -> A+2X, rate k2\n", + "- B -> X, rate k3\n", + "- X -> B, rate k4" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import pints\n", + "import pints.toy.stochastic\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import math\n", + "from tqdm import tqdm" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext line_profiler" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n", + "\n", + " 0%| | 0/1000 [00:00" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "for value in valuess[:10]:\n", + " plt.step(times, value)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "376 ms ± 43.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" + ] + } + ], + "source": [ + "%timeit model.simulate(k, times)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "31.7 ms ± 1.31 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" + ] + } + ], + "source": [ + "%timeit model.simulate(k, times, approx_tau=0.0125)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "valuess = np.array(valuess)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "No handles with labels found to put in legend.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAATVElEQVR4nO3df/BldX3f8edL2bhJoALLQggL+W4M1mDSANlYZmyM1kYFomisDjaN1NCsmcEp1LTjaiaNncbpmghJ7Q/sGhhXowIdtG4C+YEM0XEa1F2k/BApi27gK5vdlaZCalHAd/+45/vxsnv3+727+73fc7/7fT5m7pxzP+ece9575uy+9nzOuZ+bqkKSJIBn9V2AJGl6GAqSpMZQkCQ1hoIkqTEUJEnNMX0XcCROOumkmpmZ6bsMSVpWduzY8Y2qWjtq2bIOhZmZGbZv3953GZK0rCT5q4Mts/tIktRMLBSSnJ7ktiT3Jbk3yeVd+7uTfD3Jnd3rgqFt3plkZ5L7k7xyUrVJkkabZPfRU8CvV9UdSY4DdiS5pVv2e1X1vuGVk5wFXAy8EPhh4NNJnl9VT0+wRknSkImFQlXtBnZ3848nuQ84bZ5NLgKuq6pvA19LshN4EfCXk6pRkpa7J598ktnZWZ544okDlq1evZp169axatWqsT9vSW40J5kBzgE+D7wYeFuSNwPbGVxN/A2DwLh9aLNZRoRIko3ARoAzzjhjonVL0rSbnZ3luOOOY2ZmhiStvap49NFHmZ2dZf369WN/3sRvNCc5FrgRuKKqHgOuBp4HnM3gSuLKuVVHbH7AaH1VtaWqNlTVhrVrRz5RJUkrxhNPPMGaNWueEQgASVizZs3IK4j5TDQUkqxiEAgfrapPAFTVnqp6uqq+C3yQQRcRDK4MTh/afB3wyCTrk6Sjwf6BsFD7fCb59FGAa4D7quqqofZTh1Z7HXBPN78NuDjJc5KsB84EvjCp+iRJB5rkPYUXA78M3J3kzq7tXcCbkpzNoGtoF/BWgKq6N8kNwJcZPLl0mU8eSdLSmuTTR59j9H2Cm+fZ5j3AeyZVk7SUZjbd1Mt+d22+sJf9qj9VNbKr6HB+RM1vNEvSMrZ69WoeffTRAwJg7umj1atXH9LnLeuxjyRppVu3bh2zs7Ps27fvgGVz31M4FIaCJC1jq1atOqTvISzEUNBRra9+fWm58p6CJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUTCwUkpye5LYk9yW5N8nlXfuJSW5J8kA3PaFrT5L3J9mZ5K4k506qNknSaJO8UngK+PWq+nHgPOCyJGcBm4Bbq+pM4NbuPcD5wJndayNw9QRrkySNMLFQqKrdVXVHN/84cB9wGnARsLVbbSvw2m7+IuDDNXA7cHySUydVnyTpQEtyTyHJDHAO8HnglKraDYPgAE7uVjsNeHhos9mubf/P2phke5Lt+/btm2TZkrTiTDwUkhwL3AhcUVWPzbfqiLY6oKFqS1VtqKoNa9euXawyJUlMOBSSrGIQCB+tqk90zXvmuoW66d6ufRY4fWjzdcAjk6xPkvRMk3z6KMA1wH1VddXQom3AJd38JcCnhtrf3D2FdB7wzbluJknS0jhmgp/9YuCXgbuT3Nm1vQvYDNyQ5FLgIeAN3bKbgQuAncC3gLdMsDZJ0ggTC4Wq+hyj7xMAvHzE+gVcNql6JEkL8xvNkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJzTF9F6CVYWbTTX2XsGL0dax3bb6wl/1qcXmlIElqxgqFJDcmuTCJISJJR7Fx/5G/GvgnwANJNid5wQRrkiT1ZKxQqKpPV9UvAecCu4BbkvyPJG9JsmqSBUqSls7Y3UFJ1gD/DPjnwJeA/8AgJG45yPrXJtmb5J6htncn+XqSO7vXBUPL3plkZ5L7k7zyMP88kqQjMNbTR0k+AbwA+Ajw6qra3S26Psn2g2z2IeA/AR/er/33qup9+33+WcDFwAuBHwY+neT5VfX0WH8KSdKiGPeR1D+oqpuHG5I8p6q+XVUbRm1QVZ9NMjPm518EXFdV3wa+lmQn8CLgL8fcXpK0CMbtPvrtEW2H+w/225Lc1XUvndC1nQY8PLTObNd2gCQbk2xPsn3fvn2HWYIkaZR5QyHJDyX5aeD7k5yT5Nzu9VLgBw5jf1cDzwPOBnYDV87tasS6NeoDqmpLVW2oqg1r1649jBIkSQezUPfRKxncXF4HXDXU/jjwrkPdWVXtmZtP8kHgj7u3s8DpQ6uuAx451M+XJB2ZeUOhqrYCW5O8vqpuPNKdJTl16Cb164C5J5O2AR9LchWDG81nAl840v1Jkg7NvKGQ5J9W1R8CM0nevv/yqrpqxGZz234ceClwUpJZ4LeAlyY5m0HX0C7grd3n3JvkBuDLwFPAZT55JElLb6Huox/spsce6gdX1ZtGNF8zz/rvAd5zqPuRJC2ehbqP/ms3/bdLU44kqU8LdR+9f77lVfUvFrccSVKfFuo+2rEkVUiSpsI4Tx9JklaIhbqPfr+qrkjyR4z4MllVvWZilUmSltxC3Ucf6abvm3ctSdJRYaHuox3d9DNJvo/BSKkF3F9V31mC+iRJS2jcobMvBD4APMhgnKL1Sd5aVX8yyeIkSUtr3KGzrwReVlU7AZI8D7gJMBQk6Sgy7tDZe+cCofNVYO8E6pEk9Wihp49+sZu9N8nNwA0M7im8AfjihGuTJC2xhbqPXj00vwf4uW5+H3DCgatLkpazhZ4+estSFSJJ6t+4Tx+tBi4FXgisnmuvql+ZUF2SpB6Me6P5I8APMfglts8w+GW0xydVlCSpH+OGwo9V1W8C/7cbD+lC4CcnV5YkqQ/jhsKT3fT/JPkJ4LnAzEQqkiT1Ztwvr21JcgLwmwx+T/nYbl6SdBQZKxSq6g+62c8APzq5ciRJfRqr+yjJmiT/MckdSXYk+f0kayZdnCRpaY17T+E6BsNavB74x8A3gOsnVZQkqR/j3lM4sar+3dD7307y2kkUJEnqz7hXCrcluTjJs7rXGxmMkipJOoosNCDe4wwGwAvwduAPu0XPAv4W+K2JVidJWlILjX103FIVIknq37j3FEjyGuAl3du/qKo/nkxJkqS+jPtI6mbgcuDL3evyrk2SdBQZ90rhAuDsqvouQJKtwJeATZMqTJK09MZ9+gjg+KH55y52IZKk/o17pfDvgS8luY3Bk0gvAd45saokSb1YMBSSBPgccB7wMwxC4R1V9dcTrk2StMQWDIWqqiT/vap+msEIqZKko9S49xRuT/IzE61EktS7cUPhZQyC4cEkdyW5O8ld822Q5Noke5PcM9R2YpJbkjzQTU/o2pPk/Ul2dp9/7uH/kSRJh2vcUDifwe8o/EPg1cAvdNP5fAh41X5tm4Bbq+pM4Fa+90jr+cCZ3WsjcPWYdUmSFtFCYx+tBn4N+DHgbuCaqnpqnA+uqs8mmdmv+SLgpd38VuAvgHd07R+uqmJwRXJ8klOravd4fwxJ0mJY6EphK7CBQSCcD1x5hPs7Ze4f+m56ctd+GvDw0HqzXdsBkmxMsj3J9n379h1hOZKkYQs9fXRWVf0kQJJrgC9MqI6MaKtRK1bVFmALwIYNG0auI0k6PAtdKTw5NzNut9EC9iQ5FaCb7u3aZ4HTh9ZbBzyyCPuTJB2ChULhp5I81r0eB/7e3HySxw5jf9uAS7r5S4BPDbW/uXsK6Tzgm95PkKSlt9DvKTz7cD84yccZ3FQ+Kcksgx/k2QzckORS4CHgDd3qNzMYdG8n8C3gLYe7X0nS4Rv79xQOVVW96SCLXj5i3QIum1QtkqTxHMooqZKko5yhIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJzcS+0SxpZZnZdFNv+961+cLe9n208UpBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGYS5WkD6HIZC0PHilIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLU9DIgXpJdwOPA08BTVbUhyYnA9cAMsAt4Y1X9TR/1SdJK1eeVwsuq6uyq2tC93wTcWlVnArd27yVJS2iauo8uArZ281uB1/ZYiyStSH2FQgF/nmRHko1d2ylVtRugm57cU22StGL19SM7L66qR5KcDNyS5CvjbtiFyEaAM844Y1L1SdKK1MuVQlU90k33Ap8EXgTsSXIqQDfde5Btt1TVhqrasHbt2qUqWZJWhCUPhSQ/mOS4uXngFcA9wDbgkm61S4BPLXVtkrTS9dF9dArwySRz+/9YVf1pki8CNyS5FHgIeEMPtUnSirbkoVBVXwV+akT7o8DLl7oeSdL3TNMjqZKknhkKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktT08XOckrSoZjbd1Mt+d22+sJf9TpJXCpKkxlCQJDWGgiSpMRQkSY2hIElqfPqoB309KSFJC/FKQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktT4jWZJOkx9jk4wqd9ymLorhSSvSnJ/kp1JNvVdjyStJFN1pZDk2cB/Bn4emAW+mGRbVX15sffl+EOSdKBpu1J4EbCzqr5aVd8BrgMu6rkmSVoxpupKATgNeHjo/Szw94dXSLIR2Ni9/dsk9y9RbcNOAr7Rw36P1HKseznWDMuz7uVYMyzPuo+45rz3iPb/IwdbMG2hkBFt9Yw3VVuALUtTzmhJtlfVhj5rOBzLse7lWDMsz7qXY82wPOue5pqnrftoFjh96P064JGeapGkFWfaQuGLwJlJ1if5PuBiYFvPNUnSijFV3UdV9VSStwF/BjwbuLaq7u25rFF67b46Asux7uVYMyzPupdjzbA8657amlNVC68lSVoRpq37SJLUI0NBktQYCvNI8neT3Dn0eizJFUneneTrQ+0XTEGt1ybZm+SeobYTk9yS5IFuekLXniTv74YSuSvJuVNW9+8m+UpX2yeTHN+1zyT5f0PH/QNTVPNBz4kk7+yO9f1JXtlHzV0do+q+fqjmXUnu7Nqn5VifnuS2JPcluTfJ5V37VJ/b89Q91ec2AFXla4wXgxvff83gSx/vBv5V3zXtV99LgHOBe4bafgfY1M1vAt7bzV8A/AmD74WcB3x+yup+BXBMN//eobpnhtebsppHnhPAWcD/BJ4DrAceBJ49LXXvt/xK4N9M2bE+FTi3mz8O+F/dMZ3qc3ueuqf63K4qrxQOwcuBB6vqr/ouZJSq+izwv/drvgjY2s1vBV471P7hGrgdOD7JqUtT6TONqruq/ryqnure3s7g+ypT4yDH+mAuAq6rqm9X1deAnQyGc1ly89WdJMAbgY8vaVELqKrdVXVHN/84cB+DkQ+m+tw+WN3Tfm6D3UeH4mKe+Rfmbd0l4LVzl65T6JSq2g2DkxQ4uWsfNZzIaUtc27h+hcH//OasT/KlJJ9J8rN9FXUQo86J5XKsfxbYU1UPDLVN1bFOMgOcA3yeZXRu71f3sKk8tw2FMWTwRbrXAP+ta7oaeB5wNrCbwWX3crLgcCLTIMlvAE8BH+2adgNnVNU5wNuBjyX5O33Vt5+DnRPL4lgDb+KZ/+mZqmOd5FjgRuCKqnpsvlVHtPV2vA9W9zSf24bCeM4H7qiqPQBVtaeqnq6q7wIfpKfugDHsmbt07qZ7u/apH04kySXALwC/VF2na9cF82g3v4NB//zz+6vye+Y5J5bDsT4G+EXg+rm2aTrWSVYx+If1o1X1ia556s/tg9Q99ee2oTCeZ/wvar8+ytcB9xywxXTYBlzSzV8CfGqo/c3dkxrnAd+cuxSfBkleBbwDeE1VfWuofW0Gv7lBkh8FzgS+2k+VzzTPObENuDjJc5KsZ1DzF5a6vgX8I+ArVTU71zAtx7q713ENcF9VXTW0aKrP7YPVvSzO7b7vdE/7C/gB4FHguUNtHwHuBu5icBKeOgV1fpzBJeiTDP63dCmwBrgVeKCbntitGwY/ZvRg9+fYMGV172TQL3xn9/pAt+7rgXsZPM1zB/DqKar5oOcE8Bvdsb4fOH+ajnXX/iHg1/Zbd1qO9T9g0P1z19D5cMG0n9vz1D3V53ZVOcyFJOl77D6SJDWGgiSpMRQkSY2hIElqDAVJUmMoSIukGxnza0lO7N6f0L3/kb5rk8ZlKEiLpKoeZjDcxeauaTOwpaZ0EEVpFL+nIC2ibmiDHcC1wK8C51TVd/qtShrfMX0XIB1NqurJJP8a+FPgFQaClhu7j6TFdz6D4SR+ou9CpENlKEiLKMnZwM8z+NWvf9nXjxdJh8tQkBZJNzLm1QzGzn8I+F3gff1WJR0aQ0FaPL8KPFRVt3Tv/wvwgiQ/12NN0iHx6SNJUuOVgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTm/wMVD8zJmmkRYwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.hist(valuess[:,-1])\n", + "plt.legend()\n", + "plt.xlabel('X')\n", + "plt.ylabel('Probability')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Given the stochastic nature of this model, every iteration returns a different result. However, averaging the concentration values at each time step, produces a reproducible result which tends towards a deterministic function as the the number of iterations tends to infinity (Erban et al., 2007): $ n_0e^{-kt} $\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "N = 100\n", + "average = np.zeros(len(times))\n", + "for i in range(N):\n", + " values = np.asarray(model.simulate(k, times))[0,:]\n", + " average += (values/N)\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "average.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "plt.plot(times, average, label = 'deterministic mean of A(t)')\n", + "plt.title('stochastic degradation across different iterations')\n", + "plt.xlabel('time')\n", + "plt.ylabel('concentration (A(t))')\n", + "plt.legend(loc = 'upper right')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The deterministic mean (from above) is plotted with the deterministic standard deviation. \n", + "The deterministic variance of this model is given by: $e^{-2kt}(-1 + e^{kt})n_0$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mean = model.mean(k, times)\n", + "variance = model.variance(k, times)\n", + "std_dev = np.sqrt(variance)\n", + "\n", + "plt.plot(times, mean, '-', label = 'deterministic mean of A(t)')\n", + "plt.plot(times, mean + std_dev, '--', label = 'standard deviation upper bound')\n", + "plt.plot(times, mean - std_dev, '--', label = 'standard deviation lower bound')\n", + "plt.legend(loc = 'upper right')\n", + "plt.xlabel('time')\n", + "plt.ylabel('concentration (A(t))')\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.4" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/toy/model-sir.ipynb b/examples/toy/model-sir.ipynb index 26ec07ee0..50b446ad4 100644 --- a/examples/toy/model-sir.ipynb +++ b/examples/toy/model-sir.ipynb @@ -285,9 +285,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.7" + "version": "3.7.4" } }, "nbformat": 4, - "nbformat_minor": 1 + "nbformat_minor": 4 } diff --git a/pints/tests/test_toy_stochastic_degradation_model.py b/pints/tests/test_toy_stochastic_degradation_model.py index a7262c90b..3c48feb67 100755 --- a/pints/tests/test_toy_stochastic_degradation_model.py +++ b/pints/tests/test_toy_stochastic_degradation_model.py @@ -1,6 +1,6 @@ #!/usr/bin/env python3 # -# Tests if the stochastic degradation (toy) model works. +# Tests if the degradation (toy) model works. # # This file is part of PINTS (https://github.com/pints-team/pints/) which is # released under the BSD 3-clause license. See accompanying LICENSE.md for @@ -10,16 +10,16 @@ import numpy as np import pints import pints.toy -from pints.toy import StochasticDegradationModel +from pints.toy.stochastic import DegradationModel -class TestStochasticDegradation(unittest.TestCase): +class TestDegradation(unittest.TestCase): """ - Tests if the stochastic degradation (toy) model works. + Tests if the degradation (toy) model works. """ def test_start_with_zero(self): # Test the special case where the initial molecule count is zero - model = StochasticDegradationModel(0) + model = DegradationModel(0) times = [0, 1, 2, 100, 1000] parameters = [0.1] values = model.simulate(parameters, times) @@ -28,7 +28,7 @@ def test_start_with_zero(self): def test_start_with_twenty(self): # Run small simulation - model = pints.toy.StochasticDegradationModel(20) + model = pints.toy.DegradationModel(20) times = [0, 1, 2, 100, 1000] parameters = [0.1] values = model.simulate(parameters, times) @@ -38,7 +38,7 @@ def test_start_with_twenty(self): self.assertTrue(np.all(values[1:] <= values[:-1])) def test_suggested(self): - model = pints.toy.StochasticDegradationModel(20) + model = pints.toy.DegradationModel(20) times = model.suggested_times() parameters = model.suggested_parameters() self.assertTrue(len(times) == 101) @@ -46,7 +46,7 @@ def test_suggested(self): def test_simulate(self): times = np.linspace(0, 100, 101) - model = StochasticDegradationModel(20) + model = DegradationModel(20) time, mol_count = model.simulate_raw([0.1]) values = model.interpolate_mol_counts(time, mol_count, times) self.assertTrue(len(time), len(mol_count)) @@ -68,27 +68,27 @@ def test_simulate(self): def test_mean_variance(self): # test mean - model = pints.toy.StochasticDegradationModel(10) + model = pints.toy.DegradationModel(10) v_mean = model.mean([1], [5, 10]) self.assertEqual(v_mean[0], 10 * np.exp(-5)) self.assertEqual(v_mean[1], 10 * np.exp(-10)) - model = pints.toy.StochasticDegradationModel(20) + model = pints.toy.DegradationModel(20) v_mean = model.mean([5], [7.2]) self.assertEqual(v_mean[0], 20 * np.exp(-7.2 * 5)) # test variance - model = pints.toy.StochasticDegradationModel(10) + model = pints.toy.DegradationModel(10) v_var = model.variance([1], [5, 10]) self.assertEqual(v_var[0], 10 * (np.exp(5) - 1.0) / np.exp(10)) self.assertAlmostEqual(v_var[1], 10 * (np.exp(10) - 1.0) / np.exp(20)) - model = pints.toy.StochasticDegradationModel(20) + model = pints.toy.DegradationModel(20) v_var = model.variance([2.0], [2.0]) self.assertAlmostEqual(v_var[0], 20 * (np.exp(4) - 1.0) / np.exp(8)) def test_errors(self): - model = pints.toy.StochasticDegradationModel(20) + model = pints.toy.DegradationModel(20) # parameters, times cannot be negative times = np.linspace(0, 100, 101) parameters = [-0.1] @@ -109,7 +109,7 @@ def test_errors(self): self.assertRaises(ValueError, model.variance, parameters_3, times) # Initial value can't be negative - self.assertRaises(ValueError, pints.toy.StochasticDegradationModel, -1) + self.assertRaises(ValueError, pints.toy.DegradationModel, -1) if __name__ == '__main__': diff --git a/pints/toy/__init__.py b/pints/toy/__init__.py index 989c21bca..db53fb8fa 100644 --- a/pints/toy/__init__.py +++ b/pints/toy/__init__.py @@ -35,4 +35,3 @@ from ._simple_egg_box import SimpleEggBoxLogPDF # noqa from ._sir_model import SIRModel # noqa from ._twisted_gaussian_banana import TwistedGaussianLogPDF # noqa -from ._stochastic_degradation_model import StochasticDegradationModel # noqa diff --git a/pints/toy/stochastic/__init__.py b/pints/toy/stochastic/__init__.py new file mode 100644 index 000000000..61e23b3b9 --- /dev/null +++ b/pints/toy/stochastic/__init__.py @@ -0,0 +1,16 @@ +# +# Root of the toy module. +# Provides a number of toy models and logpdfs for tests of Pints' functions. +# +# This file is part of PINTS. +# Copyright (c) 2017-2019, University of Oxford. +# For licensing information, see the LICENSE file distributed with the PINTS +# software package. +# +from __future__ import absolute_import, division +from __future__ import print_function, unicode_literals + +from ._degradation_model import DegradationModel # noqa +from ._markov_jump_model import MarkovJumpModel # noqa +from ._michaelis_menten_model import MichaelisMentenModel # noqa +from ._schlogl_model import SchloglModel # noqa diff --git a/pints/toy/_stochastic_degradation_model.py b/pints/toy/stochastic/_degradation_model.py similarity index 97% rename from pints/toy/_stochastic_degradation_model.py rename to pints/toy/stochastic/_degradation_model.py index 9068458e3..06176ce91 100644 --- a/pints/toy/_stochastic_degradation_model.py +++ b/pints/toy/stochastic/_degradation_model.py @@ -11,10 +11,10 @@ from scipy.interpolate import interp1d import pints -from . import ToyModel +from .. import ToyModel -class StochasticDegradationModel(pints.ForwardModel, ToyModel): +class DegradationModel(pints.ForwardModel, ToyModel): r""" Stochastic degradation model of a single chemical reaction starting from an initial molecule count :math:`A(0)` and degrading to 0 with a fixed rate @@ -64,7 +64,7 @@ class StochasticDegradationModel(pints.ForwardModel, ToyModel): https://doi.org/10.1016/0021-9991(76)90041-3 """ def __init__(self, initial_molecule_count=20): - super(StochasticDegradationModel, self).__init__() + super(DegradationModel, self).__init__() self._n0 = float(initial_molecule_count) if self._n0 < 0: raise ValueError('Initial molecule count cannot be negative.') diff --git a/pints/toy/stochastic/_markov_jump_model.py b/pints/toy/stochastic/_markov_jump_model.py new file mode 100644 index 000000000..6bda241f4 --- /dev/null +++ b/pints/toy/stochastic/_markov_jump_model.py @@ -0,0 +1,188 @@ +# +# Stochastic degradation toy model. +# +# This file is part of PINTS. +# Copyright (c) 2017-2019, University of Oxford. +# For licensing information, see the LICENSE file distributed with the PINTS +# software package. +# +from __future__ import absolute_import, division +from __future__ import print_function, unicode_literals +import numpy as np +from scipy.interpolate import interp1d +import pints +import random + +from .. import ToyModel + + +class MarkovJumpModel(pints.ForwardModel, ToyModel): + r""" + A general purpose Markov Jump model used for any systems of reactions + that proceed through jumps. We simulate a population of N different species + reacting through M different mechanisms. + + A model has three parameters: + - x_0 - an N-vector specifying the initial population of each + of the N species + - V - an NxM matrix consisting of stochiometric vectors v_i specifying + the changes to the state, x, from reaction i taking place + - a - a function from the current state, x, and reaction rates, k, + to a vector of the rates of each reaction taking place + + Simulations are performed using Gillespie's algorithm [1]_, [2]_: + + 1. Sample values :math:`r_0`, :math:`r_1`, from a uniform distribution + + .. math:: + r_0, r_1 \sim U(0,1) + + 2. Calculate the time :math:`\tau` until the next single reaction as + + .. math:: + \tau = \frac{-\ln(r)}{a_0} + + 3. Decide which reaction, i, takes place using r_1 * a_0 and iterating + through propensities + + 4. Update the state :math:`x` at time :math:`t + \tau` as: + + .. math:: + x(t + \tau) = x(t) + V[i] + + 4. Return to step (1) until no reaction can take place + + The model has one parameter, the rate constant :math:`k`. + + Extends :class:`pints.ForwardModel`, :class:`pints.toy.ToyModel`. + + Parameters + ---------- + initial_molecule_count + The initial molecule count :math:`A(0)`. + + References + ---------- + .. [1] A Practical Guide to Stochastic Simulations of Reaction Diffusion + Processes. Erban, Chapman, Maini (2007). + arXiv:0704.1908v2 [q-bio.SC] + https://arxiv.org/abs/0704.1908 + + .. [2] A general method for numerically simulating the stochastic time + evolution of coupled chemical reactions. Gillespie (1976). + Journal of Computational Physics + https://doi.org/10.1016/0021-9991(76)90041-3 + """ + def __init__(self, x0, V, a): + super(MarkovJumpModel, self).__init__() + self._x0 = np.asarray(x0) + self._V = V + self._a = a + if any(self._x0 < 0): + raise ValueError('Initial molecule count cannot be negative.') + + def n_parameters(self): + """ See :meth:`pints.ForwardModel.n_parameters()`. """ + return len(self._V) + + def simulate_raw(self, rates, max_time): + """ + Returns raw times, mol counts when reactions occur + """ + # parameters = np.asarray(parameters) + # if len(parameters) != self.n_parameters(): + # raise ValueError('This model should have only 1 parameter.') + # k = parameters[0] + + current_rates = self._a(self._x0, rates) + a_0 = sum(current_rates) + + # Initial time and count + t = 0 + x = np.array(self._x0) + + # Run gillespie SSA, calculating time until next + # reaction, deciding which reaction, and applying it + mol_count = [np.array(x)] + time = [t] + while a_0 > 0 and t <= max_time: + r_1, r_2 = random.random(), random.random() + t += -np.log(r_1) / (a_0) + s = 0 + r = 0 + while s <= r_2 * a_0: + s += current_rates[r] + r += 1 + r -= 1 + x = np.add(x, self._V[r]) + + current_rates = self._a(x, rates) + a_0 = sum(current_rates) + + time.append(t) + mol_count.append(np.array(x)) + return time, mol_count + + def simulate_approx(self, rates, max_time, tau): + assert tau > 0, "cannot tau-leap with negative tau" + current_rates = np.array(self._a(self._x0, rates)) + # Initial time and count + t = 0 + x = self._x0.copy() + N = len(rates) + # Run gillespie SSA, calculating time until next + # reaction, deciding which reaction, and applying it + mol_count = [x.copy()] + time = [t] + while any(current_rates > 0) and t <= max_time: + # Estimate number of each reaction in [t, t+tau) + k = [np.random.poisson(current_rates[i] * tau) for i in range(N)] + + # Apply the reactions + for r in range(N): + x += np.array(self._V[r]) * k[r] + + # Update rates + current_rates = np.array(self._a(x, rates)) + + # Advance Time + t += tau + time.append(t) + mol_count.append(x.copy()) + return time, mol_count + + def interpolate_mol_counts(self, time, mol_count, output_times): + """ + Takes raw times and inputs and mol counts and outputs interpolated + values at output_times + """ + # Interpolate as step function, decreasing mol_count by 1 at each + # reaction time point + interp_func = interp1d(time, mol_count, kind='previous', axis=0, + fill_value="extrapolate", bounds_error=False) + + # Compute molecule count values at given time points using f1 + # at any time beyond the last reaction, molecule count = 0 + values = interp_func(output_times) + return values + + def simulate(self, parameters, times, approx_tau=None): + """ See :meth:`pints.ForwardModel.simulate()`. """ + times = np.asarray(times) + if np.any(times < 0): + raise ValueError('Negative times are not allowed.') + if approx_tau is None: + # run Gillespie + time, mol_count = self.simulate_raw(parameters, max(times)) + else: + if (not approx_tau) or approx_tau <= 0: + ValueError("You must provide a positive value for approx_tau\ + to use tau-leaping approximation") + # Run Euler tau-leaping + time, mol_count = self.simulate_approx(parameters, max(times), + approx_tau) + # interpolate + values = self.interpolate_mol_counts(np.asarray(time), + np.asarray(mol_count), times) + return values + diff --git a/pints/toy/stochastic/_michaelis_menten_model.py b/pints/toy/stochastic/_michaelis_menten_model.py new file mode 100644 index 000000000..2861717a0 --- /dev/null +++ b/pints/toy/stochastic/_michaelis_menten_model.py @@ -0,0 +1,42 @@ +# +# Stochastic degradation toy model. +# +# This file is part of PINTS. +# Copyright (c) 2017-2019, University of Oxford. +# For licensing information, see the LICENSE file distributed with the PINTS +# software package. +# +from __future__ import absolute_import, division +from __future__ import print_function, unicode_literals + +from . import MarkovJumpModel + + +class MichaelisMentenModel(MarkovJumpModel): + r""" + Simulates the Michaelis Menten Dynamics using Gillespie. + + This system of reaction involves 4 chemical species with + inital counts x_0, and reactions: + - X1+X2 -> X3 with rate k1 + - X3 -> X1+X2 with rate k2 + - X3 -> X2+X4 with rate k3 + """ + def __init__(self, x_0): + mat = [[-1, -1, 1, 0], + [1, 1, -1, 0], + [0, 1, -1, 1]] + super(MichaelisMentenModel, self).__init__(x_0, + mat, self._propensities) + + @staticmethod + def _propensities(xs, ks): + return [ + xs[0] * xs[1] * ks[0], + xs[2] * ks[1], + xs[2] * ks[2] + ] + + def n_parameters(self): + """ See :meth:`pints.ForwardModel.n_parameters()`. """ + return 3 diff --git a/pints/toy/stochastic/_schlogl_model.py b/pints/toy/stochastic/_schlogl_model.py new file mode 100644 index 000000000..7f904e590 --- /dev/null +++ b/pints/toy/stochastic/_schlogl_model.py @@ -0,0 +1,50 @@ +# +# Stochastic Schlogl toy model. +# +# This file is part of PINTS. +# Copyright (c) 2017-2019, University of Oxford. +# For licensing information, see the LICENSE file distributed with the PINTS +# software package. +# +from __future__ import absolute_import, division +from __future__ import print_function, unicode_literals + +from . import MarkovJumpModel + + +class SchloglModel(MarkovJumpModel): + r""" + Simulates the Schlögl System using Gillespie. + + This system of reaction involved 4 chemical species with + inital counts x_0, and reactions: + - A + 2X -> 3X, rate k1 + - 3X -> A+2X, rate k2 + - B -> X, rate k3 + - X -> B, rate k4 + """ + def __init__(self, x_0): + self.a_count = 1e5 + self.b_count = 2e5 + # We are only concered with the change in X concentration + mat = [[1], + [-1], + [1], + [-1]] + super(SchloglModel, self).__init__([x_0], mat, self._propensities) + + def simulate(self, parameters, times, approx=None, approx_tau=None): + return super(SchloglModel, self).simulate( + parameters, times, approx_tau)[:, 0] + + def _propensities(self, xs, ks): + return [ + self.a_count * xs[0] * xs[0] * ks[0], + xs[0] * xs[0] * xs[0] * ks[1], + self.b_count * ks[2], + xs[0] * ks[3] + ] + + def n_parameters(self): + """ See :meth:`pints.ForwardModel.n_parameters()`. """ + return 4