test_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import optimize as opt\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from flowstock import Akagi\n",
    "from flowstock.data import load_csv\n",
    "from flowstock.data import distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "centroid_data_fname = \"/Users/easther/unixland/thinking/covid_models/data/area/statistical-area-2-2018-centroid-true.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "telco_data_fname = \"/Users/easther/unixland/thinking/covid_models/data/telco/pop_data_2020-04-01.dat\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "center = \"Christchurch Central\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "date_first = datetime(2020, 2, 18, 8)\n",
    "date_last = datetime(2020, 2, 18, 9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = 15000 # meters but K is in km to avoid numerical issues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "centroid_data = load_csv.load_centroid_data_2018(centroid_data_fname)\n",
    "centroid_data_shrunk = distances.remove_furthest_points(centroid_data, D, center)\n",
    "distance_table = distances.calculate_distances(centroid_data_shrunk)\n",
    "d = (distance_table.to_numpy()[:, 2:] / 1000).astype(float)\n",
    "\n",
    "region_codes = distance_table['region_code']\n",
    "\n",
    "telco = load_csv.load_telco_data(telco_data_fname)\n",
    "\n",
    "telco_subset = pd.pivot_table(\n",
    "    telco[\n",
    "        (telco[\"region_code\"].isin(region_codes)) &\n",
    "        ((telco[\"time\"] == date_first) |\n",
    "        (telco[\"time\"] == date_last))\n",
    "    ],\n",
    "    index=\"time\",\n",
    "    columns=\"region_code\",\n",
    "    values=\"count\",\n",
    ")\n",
    "\n",
    "N = telco_subset.to_numpy()*1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(167, 167)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 167)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dtype('int64')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N.dtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 80"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x102ddec2e8>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAAmCAYAAADtNdFAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAKqUlEQVR4nO2de4xcVR3HP7957HYfbXeXZdulS+hDWq0YClYpAkZbnpVAjJJAMEKEkPCHEY0PmiYmxGgCGjUmBiSKMYIoImBpNCgPTUxMaYsUCttCG9qyfW1b6Jbudl+zP//4/e4yM07p6s7u3Mjvk2zu3HPPOfd7vuec3505c+euqCpBEARBesnUWkAQBEHw3kSgDoIgSDkRqIMgCFJOBOogCIKUE4E6CIIg5USgDoIgSDkTCtQicqWIbBeRHSJy51SLCoIgCN5FTnUftYhkgdeAy4AeYCNwg6q+OvXygiAIgom8o34C6ALWqeow8Fvg2ilVFQRBEIyTm0Ce54FR4Gzf7wEueK8CdfkmnVHfgowUABhrsNNoxygA2Z7S64NmBAAZHE4q8ANastWslZPhUa+3DoCM74/MsnL5YyNWbszLjfh+44yS846fL2f6hlutfN2RIc8gXs+YHW+38rkTVq94/TI0UlLPWD7j7bLk7IDp03zW8o8WStpTaMx6Pk/PWfpYzs6fHR7zcmNer1fsm6G2jPvgyS5HCqWfljLl9bieQoOVzx03nZR9yBqrs3yjTVJSf67PfNL6vLfHjx8ftAz1dV6/lc8MWvsS35J6E53jOo4M+AmS4WnHR2fVmx+Dpp/+E1bez5Mw2mR6skPeXh+HOuwGFY2DUfdebUO+P6nb2jDW0mh1jJb1uXs41GYac97krLdRBz2h2cv7GB0fwwU/j4+FzAdsd2yPj30fc8kYHKs3gcOtllz3lm0zQ95nyZgYsf1Cs3s14G32OZTMmcTzkZnu+UBppyftS3Rq1vtq2Do/6YvccR8MSR81W/25fksfqzN/MgNDJTrGyVt7R2ZbvnzfSGl7Cu5n3mOIx4rMcKG0vmSbLXvvmZR3nyUZU2X5xvvF840zHgMS3e6jj/lxHQXvB/ep/PzjY3ks8cnK507Y8ZGmLMPvvMXoYL9QgYkE6m3Ah8vS/mO9RERuA24DmFE/m48vu538/qMADCw+HYCFd20DYN+tXVbGGzHqkyHXvcsqO2OON9IHy6B18uCiDgDqt7wBwNCyBQDM2NELwMHLrN45z+yz8idssoweOAhAZslSS08mwY49lt5mo7/nc2cBMO9B00nGJ7FPuj23fgSAjs02+JMAWt+91/LNagZg+IzZdl4PUE3ddv6xFjuu3vnig+vweZa/rbvfZM+1QDLcZINnZo+1P9973OqZ4QHQJ8XO62zbcMAnkzev4wXTmQT8hr1WPnO4z8rPtXYfWj7LfPuHz35HBqzdJxa1A7D/YpucTT2mu+PJnVb/mdYvA/OsH5uefsXKnzUPgLeXtQHQstXGQzIJB7rMj4aDFnD3XzQTgM57NwOQbW0xvxKfrlpkPm05Zukv2upb5oNLbNtn/vWusnFwmufLHrbt6O43Ld/iD423sedy8yC5qHY9bd6wdYe1/eJzAah7xwLISKNNmaYXra5Xv2PnmvuspbduPgxA4TXzZviijwLQuM3HwMwmq98varLnAAB/+ttzAKxedZ2l+8U/CUzHz7E5NOsOO+/A98zbhm4rX5hjXmV223n6V9jcaN64287nY25wqZXLH7W+ff0L5vn8J5Oru23qjvhF8LhttbnB9vdb+w59xvqi41kb+8kFp+8C82PWy5ZveJ6N7brN5qcWSgOhdNrYObhqLgBz/2JzV/3NmvjFOGlfocHfTO3zseQXjvE3Y+2tJXqk7x0r97blz7afZvmbGkp0FNqsX7LHzJdk7Ce+jQdyvxAOLbT+qN91xI4ftTEms2eWljtm56fVfEguBH3LO90nK997yelse/xHnIyJLH30AJ1F+13AvvJMqnq/qi5X1eX5XNMEqg2CIAgmwkQC9UZgAZAXkTrgemBdeSYRuU1ENonIppHR/irLDIIgeP9yyrs+AETkZuBn2LvrB1T1u6fIfwjoBw5XQeNU0U7omwyhb3KEvsnx/6jvLFU9vdKBiQbq+cB6VT1nomcUkU2qunyi+aeb0Dc5Qt/kCH2T4/2m75RLHyLyMPBPYImI9IjILdU6eRAEQXBqTnnXh6reMB1CgiAIgspM5bM+7p/CuqtB6JscoW9yhL7J8b7SN6E16iAIgqB2xNPzgiAIUk7VA3XanrQnImeKyHMi0i0ir4jIVzy9TUT+KiKv+7a1xjqzIvIvEVnv+wtEZIPr+53fw14rbS0i8qiIbHMfL0yTfyLyVe/brSLysIjMqLV/IvKAiPSKyNaitIqeifETnzMvicj5NdL3fe/jl0TkcRFpKTq2xvVtF5EraqGv6NjXRURFpN33U+Gfp3/ZPXpFRO4pSp+cf6patT8gC+wEFgJ1wBZgaTXP8T9o6gTO99czsScBLgXuAe709DuBu2us82vAb7DbIAEeAa731/cBt9dQ26+AW/11HdCSFv+AecAbQEORbzfX2j/gk8D5wNaitIqeAauBP2M/4F4BbKiRvsuBnL++u0jfUp/L9diP33YC2enW5+lnAk8Bu4H2lPn3aeBpoN73O6rlX7XFXwg8VbS/Blgz1ab9lxr/iD2ydTvQ6WmdwPYaauoCngFWAut9wB0umjQlvk6ztlkeCKUsPRX+eaB+E2jD7mJaD1yRBv+A+WUTuaJn2I/JbqiUbzr1lR37LPCQvy6Zxx4oL6yFPuBR4FxgV1GgToV/2JuDSyvkm7R/1V76SCZNQo+npQL/4c55wAZgjqruB/BtR+2U8WPgm4A/TonTgKOq6o/kqqmPC4FDwC99aebnItJESvxT1b3AD4A9wH6gD9hMevwr5mSepXHefAl7lwop0Sci1wB7VXVL2aFU6AMWA5f4ktvfReRjnj5pfdUO1JUe0ZeK20pEpBn4A3CHqh6rtZ4EEbka6FXVzcXJFbLWyscc9hHvXlU9D3s0QM2/e0jwdd5rsY+UZwBNwFUVsqZiHJ6ENPU3IrIWe7TxQ0lShWzTqk9EGoG1wLcrHa6QVgv/ckArtvzyDeARERGqoK/agboHW0NKqPikvelGRPJYkH5IVR/z5IMi0unHO4HeGsm7CLhGRHZh/5RhJfYOu0VEkh8k1dLHHqBHVTf4/qNY4E6Lf5cCb6jqIVUdAR4DPkF6/CvmZJ6lZt6IyE3A1cCN6p/TSYe+RdjFeIvPlS7gBRGZmxJ9uI7H1Hge+4TcXg191Q7UG4Gz/Rv3kz5pbzrxK9ovgG5V/WHRoXXATf76JmztetpR1TWq2qWq8zG/nlXVG4HngM+nQN8B4E0RWeJJq4BXSYl/2JLHChFp9L5O9KXCvzJO5tk64It+98IKoC9ZIplORORK4FvANao6UHRoHXC9iNSLyALsn4g8P53aVPVlVe1Q1fk+V3qwmwQOkBL/sP+GtRJARBZjX7wfphr+TcEC+2rszoqdwNqpXtCfgJ6LsY8ZLwEv+t9qbB34GeB137alQOunePeuj4XemTuA3+PfJNdI1zJgk3v4BPbxLjX+AXdh/+BiK/Br7Nv1mvoHPIytmY9gQeWWk3mGfTT+qc+Zl4HlNdK3A1tLTebJfUX517q+7cBVtdBXdnwX736ZmBb/6oAHfRy+AKysln/xy8QgCIKUE79MDIIgSDkRqIMgCFJOBOogCIKUE4E6CIIg5USgDoIgSDkRqIMgCFJOBOogCIKUE4E6CIIg5fwb7hSYoq5FiVsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1032134d30>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEDCAYAAAA849PJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3hUddr/8fdNCT30HkLvRYUAIisi6y4IiiK4q2tD3UW3+qyrFMW119VVt6m4a9fVlSKIBRuKBUtiSaH3Tug9kGTu3x+ZfX48McokTHJmMp/XdeVi5nzPzLm/JJzPKZkbc3dERCTxVAm6ABERCYYCQEQkQSkAREQSlAJARCRBKQBERBKUAkBEJEHFXQCY2RNmlmtm2RGsm2pm883sKzPLNLORFVGjiEg8iLsAAJ4CRkS47lTgP+5+EnAB8I/yKkpEJN7EXQC4+wJg59HLzKyjmb1pZhlm9qGZdfvv6kBy+HF9YFMFlioiEtOqBV1AlEwDrnb35WY2kKIj/WHALcBbZvZboA5wRnAliojElrgPADOrC5wCvGxm/11cI/znhcBT7v6AmQ0CnjWzXu4eCqBUEZGYEvcBQNFlrN3ufmIJY1cSvl/g7gvNrCbQBMitwPpERGJS3N0DKM7d9wKrzex8ACtyQnh4HfDD8PLuQE1gWyCFiojEGIu3bqBm9m9gKEVH8luBm4H3gEeAlkB14EV3v83MegCPA3UpuiE80d3fCqJuEZFYE3cBICIi0RH3l4BERKRs4uomcJMmTbxdu3ZBlyEiElcyMjK2u3vT4svjKgDatWtHenp60GWIiMQVM1tb0nJdAhIRSVAKABGRBKUAEBFJUAoAEZEEpQAQEUlQCgARkQSlABARSVAKABGRGLZx9yFufTWHgsLod7GPqw+CiYgkilDIef6ztdzzxhJCDmNOak2flAZR3YYCQEQkxqzatp/JM7L4fM1OTu3chLvG9KZNo9pR344CQEQkRhQUhnj8w9U8+M4yalarwp/G9WFcvxSO+t8Oo0oBICISA3I27WHSjEyyN+5lRM8W3HZuT5rVq1mu21QAiIgEKC+/kL++t5xHP1hFw9pJPHJRX87s3bJCtq0AEBEJSMbanUycnsnKbQcY2zeFm87qToPaSRW2fQWAiEgFO3C4gD/NW8rTC9fQqn4tnr5iAKd1+Va7/nKnABARqUALlm1jyswsNu05xGWD2nHd8K7UrRHMrlgBICJSAXYfPMIdry1mesYGOjStw8tXDSKtXaNAa1IAiIiUszeyNnPT7Bx2HTzCr0/vyG+HdaZm9apBl6UAEBEpL7n78rh5dg5vZG+hZ6tknr6iPz1b1Q+6rP+lABARiTJ3Z3rGBu54bTGH8guZOKIrvzi1A9Wrxlb7tYirMbOqZvaVmc0tYWyImX1pZgVmNu6o5W3NLMPMvjazHDO7+qixfmaWZWYrzOwvVl4fdRMRqUDrdx7k0ic+5/rpmXRpXpc3rjmVXw3tFHM7fyjdGcA1wGIguYSxdcB44LpiyzcDp7j7YTOrC2Sb2Rx33wQ8AkwAPgVeB0YAb5SufBGR2BAKOc8sXMN985ZiwO3n9OSigW2pUiV2j20jCgAzSwFGAXcC1xYfd/c14fVCxZYfOeppDcJnHGbWEkh294Xh588A56IAEJE4tCJ3H5NmZJGxdhendWnKnWN6kdIw+s3boi3SM4CHgIlAvdJuwMzaAK8BnYDr3X2TmaUBG45abQPQ+jteP4GiMwVSU1NLu3kRkXKTXxhi2oJVPPzOcmrXqMqff3ICY05qXW7N26LtmAFgZmcBue6eYWZDS7sBd18P9DGzVsArZjYdKOlvx7/j9dOAaQBpaWklriMiUtGyN+5h4vRMFm3ey6jeLblldE+a1qsRdFmlEskZwGBgtJmNBGoCyWb2nLtfXJoNhY/8c4BTgY+BlKOGU4BNpXk/EZEg5OUX8vC7y5m2YBWN6iTx6MX9GNGrRdBllckxb0u7+xR3T3H3dsAFwHuR7vzNLMXMaoUfN6QoTJa6+2Zgn5mdHP7tn0uB2WWdhIhIRfhizU5GPvwhj7y/knF9U3jn96fF7c4fjuNzAGZ2G5Du7nPMrD8wC2gInG1mt7p7T6A78ICZOUWXfe5396zwW/wSeAqoRdHNX90AFpGYtP9wAfe9uYRnFq4lpWEtnrtyID/o3CToso6bucfPZfW0tDRPT08PugwRSSDzl+Zy48wsNu/N4/JT2nPd8C7UToqvz9CaWYa7pxVfHl+zEBGpILsOHOH2uYuY+dVGOjWry/SrT6Ff24ZBlxVVCgARkaO4O69nbeHmOdnsPpjP74Z14tfDOlGjWvDN26JNASAiEpa7N4+pr2Tz1qKt9G5dn2euGEiPViU1P6gcFAAikvDcnZfTN3D7a4s4UhBiypnduPIH7akWg/17okkBICIJbd2Og0yZlcnHK3YwoH0j7h3bh/ZN6gRdVoVQAIhIQioMOU99sob75y2lahXjjnN78bMBqTHdvC3aFAAiknCWb93HxBmZfLVuN6d3bcqdY3rTqkGtoMuqcAoAEUkYRwpCPPrBSv723grq1KjKQz89kXNObBU3zduiTQEgIgkhc8NuJk7PZMmWfZx9QituPrsHTerGV/O2aFMAiEilduhIIQ+9s4zHP1xF03o1ePzSNH7Uo3nQZcUEBYCIVFqfrtrB5BmZrNlxkAsHtGHKyO4k16wedFkxQwEgIpXOvrx87nljCc9/to7URrV54ecDOaVT/DdvizYFgIhUKu8t2cqNs7LZujePn/+gPX/4cVdqJVW+Ng7RoAAQkUph54Ej3PZqDq98vYkuzevyj4tO4aTUytW8LdoUACIS19ydVzM3c8ucHPbl5fM/Z3TmV0M7kVStcrdxiAYFgIjErS178pj6ShbvLM7lhDYNuG9sH7q2qBd0WXFDASAiccfdefGL9dz12mLyQyGmjurO5YPbUzWB2jhEgwJAROLK2h0HmDwji4WrdjCoQ2PuGdubto0To3lbtCkARCQuFIacJz9ezf1vLaV6lSrcfV5vLujfJmHbOESDAkBEYt7SLUXN275Zv5szujfjjnN706J+zaDLinsKABGJWUcKQvx9/gr+8f4K6tWszl8uPImz+7TUUX+UKABEJCZ9vX43E6d/w7Kt+zn3xFb88eyeNKqTFHRZlYoCQERiyqEjhTzw1lKe+Hg1zZNr8sT4NIZ1U/O28hDxJyXMrKqZfWVmc0sYG2JmX5pZgZmNO2r5iWa20MxyzCzTzH561NhTZrbazL4Of514/NMRkXj2ycrtDH9oAf/8aDUXDkjlrd8P0c6/HJXmDOAaYDGQXMLYOmA8cF2x5QeBS919uZm1AjLMbJ677w6PX+/u00tZs4hUMnvz8rn79cX8+/P1tGtcmxcnnMzJHRoHXValF1EAmFkKMAq4E7i2+Li7rwmvFyq2fNlRjzeZWS7QFNiNiAjw9qKtTH0li237DnPVkA78zxld1LytgkR6BvAQMBEo82eszWwAkASsPGrxnWb2R+BdYLK7Hy7hdROACQCpqall3byIxJjt+w9zy5wc5mZupluLejx+aRp9UhoEXVZCOeY9ADM7C8h194yybsTMWgLPApe7+3/PEqYA3YD+QCNgUkmvdfdp7p7m7mlNmzYtawkiEiPcnVe+2siP/vwBb+Vs5Q8/6sKc3/xAO/8ARHIGMBgYbWYjgZpAspk95+4XR7IBM0sGXgOmuvun/13u7pvDDw+b2ZN8+/6BiFQym3YfYuor2by3JJeTUouat3VuruZtQTlmALj7FIqO1jGzocB1pdj5JwGzgGfc/eViYy3dfbMVfaLjXCC7lLWLSJwIhZwXPl/HPW8soTDk/PGsHlx2Sjs1bwtYmT8HYGa3AenuPsfM+lO0o28InG1mt7p7T+AnwBCgsZmND790vLt/DTxvZk0BA74Grj6OeYhIjFq9/QCTZ2Ty2eqdDO7UmLvH9CG1ce2gyxLA3D3oGiKWlpbm6enpQZchIhEoKAzxr49W8+e3l5FUrQo3jerB+WkpauMQADPLcPe04sv1SWARibpFm/YyaUYmWRv38OMezbn93F40T1bztlijABCRqDlcUMjf3lvBI++vpEHt6vz9Z30Z2buFjvpjlAJARKIiY+0uJs3IZEXufs7r25qbRvWgoZq3xTQFgIgcl4NHCvjTvKU89ckaWibX5MnL+3N612ZBlyURUACISJl9tHw7k2dmsmHXIS4d1JaJI7pRt4Z2K/FC3ykRKbU9B/O58/VF/Cd9Ax2a1OE/Vw1iQPtGQZclpaQAEJFSeTN7CzfNzmbngSP8cmhHrvlhZ2pWV/O2eKQAEJGIbNtX1LzttazN9GiZzJPj+9Ordf2gy5LjoAAQke/l7sz8ciO3zV3EoSOFXD+8KxOGdKB61Yj/PymJUQoAEflOG3cf4oaZWXywbBv92jbk3rF96NSsbtBlSZQoAETkW0Ih57nP1nLvG0tw4Jaze3DpoHZUUfO2SkUBICL/x8pt+5k8I5Mv1uzi1M5NuGtMb9o0UvO2ykgBICIA5BeGePzDVTz0znJqVa/K/eefwNi+rdXGoRJTAIgI2Rv3MGlGJjmb9nJmrxbcek5PmtVT87bKTgEgksDy8gv563vLefSDVTSsncQjF/XlzN4tgy5LKogCQCRBpa/ZycQZmazadoBx/VKYOqo7DWqreVsiUQCIJJgDh4uatz29cA2t6tfimSsGMKRL06DLkgAoAEQSyAfLtnHDzCw27TnEZYPacf3wrtRR87aEpe+8SALYffAIt89dzIwvN9CxaR1evmoQae3UvC3RKQBEKrk3sjZz0+wcdh08wm9O78RvhnVS8zYBFAAilVbu3jz+ODuHN3O20LNVMk9f0Z+erdS8Tf4/BYBIJePuTM/YwO1zF5FXEGLSiG784tT2VFPzNilGASBSiazfeZAbZmXx4fLtDGjXiLvH9qZjUzVvk5JFfEhgZlXN7Cszm1vC2BAz+9LMCsxs3FHLTzSzhWaWY2aZZvbTo8bam9lnZrbczF4yM/0CskgZFYacJz9ezfCHFvDl2l3cfk5PXpxwsnb+8r1KcwZwDbAYSC5hbB0wHriu2PKDwKXuvtzMWgEZZjbP3XcD9wIPuvuLZvYocCXwSGknIJLoVuTuY9KMLDLW7uK0Lk2567zetG5QK+iyJA5EFABmlgKMAu4Eri0+7u5rwuuFii1fdtTjTWaWCzQ1sz3AMOBn4eGngVtQAIhELL8wxGMfrOQv766gdo2q/PknJzDmJDVvk8hFegbwEDARqFfWDZnZACAJWAk0Bna7e0F4eAPQ+jteNwGYAJCamlrWzYtUKtkb93D99EwWb97LqD4tueXsnjStVyPosiTOHDMAzOwsINfdM8xsaFk2YmYtgWeBy9w9ZCUfonhJr3X3acA0gLS0tBLXEUkUefmFPPTOch7/cBWN6iTx2CX9GN6zRdBlSZyK5AxgMDDazEYCNYFkM3vO3S+OZANmlgy8Bkx190/Di7cDDcysWvgsIAXYVPryRRLHZ6t2MHlmFqu3H+CnaW24YWR36teuHnRZEseOGQDuPgWYAhA+A7iuFDv/JGAW8Iy7v3zUe7qZzQfGAS8ClwGzS129SALYl5fPfW8u5dlP19KmUS2e//lABndqEnRZUgmU+ZMhZnabmY0OP+5vZhuA84HHzCwnvNpPgCHAeDP7Ovx1YnhsEnCtma2g6J7Av8o8C5FKav7SXIY/uIDnPlvLFYPbM+9/hmjnL1Fj7vFzWT0tLc3T09ODLkOk3O06cITb5y5i5lcb6dysLveO60Pf1IZBlyVxyswy3D2t+HJ9Elgkhrg7r2Vt5ubZOew5lM/vhnXi18M6UaOamrdJ9CkARGLE1r15TH0lm7cXbaV36/o89/OBdG9Z0ucuRaJDASASMHfnP+nrueO1xRwpCHHDyG5cMVjN26T8KQBEArRux0Emz8zkk5U7GNi+EfeO7UO7JnWCLksShAJAJACFIeepT9Zw/7ylVK1i3DmmFxf2T6VKFbVxkIqjABCpYMu27mPi9Ey+Xr+bYd2aceeYXrSsr+ZtUvEUACIV5EhBiEfeX8nf5i+nbo1qPHzBiYw+oZWat0lgFAAiFeCb9buZNCOTJVv2cfYJrbjl7B40rqvmbRIsBYBIOTp0pJAH31nGPz9cRdN6NXj80jR+1KN50GWJAAoAkXKzcOUOpszMZM2Og1w4IJUpI7uRXFPN2yR2KABEomxvXj73vLGEFz5bR9vGtXnhFwM5paP690jsUQCIRNF7S7Zyw8xscvfl8YtT23Ptj7pSK0ltHCQ2KQBEomDH/sPcNncRs7/eRNfm9Xj0kn6c2KZB0GWJfC8FgMhxcHfmfLOJW19dxL68fP7njM78amgnkqqpjYPEPgWASBlt3nOIqbOyeXdJLie0acB9Y/vQtUWZ/9tskQqnABAppVDIefGL9dz9+mLyQyGmjurO5YPbU1VtHCTOKABESmHN9gNMnpnJp6t2MqhDY+4Z25u2jdW8TeKTAkAkAgWFIZ78eA0PvL2U6lWqcM95vflp/zZq4yBxTQEgcgxLtuxl0vRMvtmwhzO6N+OOc3vTon7NoMsSOW4KAJHvcLigkL/PX8k/5q+gfq3q/PXCkzirT0sd9UuloQAQKcFX63YxaUYmy7buZ8xJrbnprB40qpMUdFkiUaUAEDnKwSMFPPDWMp74eDUtkmvyxPg0hnVT8zapnBQAImGfrNjO5JlZrNt5kItPTmXSiG7UU/M2qcQi/riimVU1s6/MbG4JY0PM7EszKzCzccXG3jSz3cVfZ2ZPmdlqM/s6/HVi2achUnZ7DuUzeUYmP/vnZ1QxeHHCydxxbm/t/KXSK80ZwDXAYiC5hLF1wHjguhLG/gTUBq4qYex6d59eihpEouqtnC1MfSWb7fsPc9VpHfj9GV2oWV3N2yQxRBQAZpYCjALuBK4tPu7ua8LrhUoYe9fMhh5XlSJRtn3/YW6Zk8PczM10a1GPf16WRp8UNW+TxBLpGcBDwEQg2o1O7jSzPwLvApPd/XCU31/k/3B3Xvl6I7e+uoiDhwv5w4+6cPXQjlSvquZtkniOGQBmdhaQ6+4ZUT6SnwJsAZKAacAk4LYStj8BmACQmpoaxc1Lotm0+xA3zspi/tJtnJRa1Lytc3M1b5PEFckZwGBgtJmNBGoCyWb2nLtffDwbdvfN4YeHzexJSr5/gLtPoyggSEtL8+PZpiSmUMh5/vN13PvGEgpDzh/P6sFlp7RT8zZJeMcMAHefQtHROuEzgOuOd+cffq+W7r7Zij5WeS6QfbzvKVLcqm37mTwji8/X7OQHnZpw93m9adOodtBlicSEMn8OwMxuA9LdfY6Z9QdmAQ2Bs83sVnfvGV7vQ6AbUNfMNgBXuvs84HkzawoY8DVw9XHOReR/FRSG+OdHq3nw7WXUqFaF+8b14fx+KWrjIHIUc4+fqyppaWmenp4edBkS4xZt2svEGd+QvXEvw3s25/ZzetEsWc3bJHGZWYa7pxVfrk8CS6VxuKCQv723gkfeX0mD2tX5x0V9ObNXCx31i3wHBYBUChlri5q3rcjdz3l9W3PTqB40VPM2ke+lAJC4duBwAfe/tZSnPllDq/q1eOry/gzt2izoskTiggJA4taHy7cxZWYWG3Yd4tJBbZk4oht1a+hHWiRS+tcicWfPwXzueG0RL2dsoEOTOvznqkEMaN8o6LJE4o4CQOLKm9lbuGl2NjsPHOFXQzvyux92VvM2kTJSAEhcyN2Xxy1zcng9aws9Wibz5Pj+9GpdP+iyROKaAkBimrsz48uN3D53EYfyC7l+eFcmDOmg5m0iUaAAkJi1YddBbpiVzYJl2+jXtiH3ju1Dp2Z1gy5LpNJQAEjMCYWcZz9dy71vLgHg1tE9ueTktlRR8zaRqFIASExZuW0/k6Znkr52F0O6NOWuMb1IaajmbSLlQQEgMSG/MMS0Bat4+N3l1KpelfvPP4GxfVurjYNIOVIASOCyN+5h0oxMcjbtZWTvFtwyuifN6ql5m0h5UwBIYPLyC/nLu8t5bMEqGtZO4tGL+zKiV8ugyxJJGAoACcQXa3YyaXomq7Yf4Px+KUwd1YP6tasHXZZIQlEASIXaf7iA+95cwjML19K6QS2euWIAQ7o0DboskYSkAJAK88GybdwwM4tNew4x/pR2XD+8K3XUvE0kMPrXJ+Vu98Ej3DZ3ETO/3EjHpnWYfvUg+rVV8zaRoCkApFy9nrWZP87OZvfBfH5zeid+M6yTmreJxAgFgJSL3L153DQ7m3k5W+nVOpmnrxhAz1Zq3iYSSxQAElXuzssZG7hj7iLyCkJMGtGNX5zanmpq3iYScxQAEjXrdx5kyswsPlqxnQHtGnHP2N50aKrmbSKxSgEgx60w5DyzcA33vbmUKga3n9uLiwakqnmbSIyL+LzczKqa2VdmNreEsSFm9qWZFZjZuGJjb5rZ7uKvM7P2ZvaZmS03s5fMLKns05CgrMjdx/mPfsKtry5iYIdGvHXtaercKRInSnNh9hpg8XeMrQPGAy+UMPYn4JISlt8LPOjunYFdwJWlqEUCll8Y4q/vLmfkwx+xavsBHvzpCTw5vj+tG9QKujQRiVBEAWBmKcAo4J8ljbv7GnfPBEIljL0L7Cv2fgYMA6aHFz0NnBt52RKkrA17OPuvH/HA28v4Uc/mvHPtaYw5KUWdO0XiTKT3AB4CJgL1orTdxsBudy8IP98AtI7Se0s5ycsv5MF3lvH4glU0qVuDxy7px/CeLYIuS0TK6JgBYGZnAbnunmFmQ6O03ZIOFf07tj8BmACQmpoapc1LaX22ageTZ2axevsBLujfhikju1O/lpq3icSzSM4ABgOjzWwkUBNINrPn3P3i49judqCBmVULnwWkAJtKWtHdpwHTANLS0koMCSk/+/LyuffNJTz36TraNKrF8z8fyOBOTYIuS0Si4JgB4O5TgCkA4TOA645z54+7u5nNB8YBLwKXAbOP5z0l+uYvyeXGWVls3pvHlT9ozx9+3IXaSfrNYZHKoswfzzSz28xsdPhxfzPbAJwPPGZmOUet9yHwMvBDM9tgZsPDQ5OAa81sBUX3BP5V1lokunYeOMLvX/qay5/6gjo1qjHjl6dw01k9tPMXqWTMPX6uqqSlpXl6enrQZVRa7s7czM3cMieHPYfy+dXpnfj16R2pUU3N20TimZlluHta8eU6pBMAtu7N48ZZ2byzeCt9Uurz3M8H0r1lctBliUg5UgAkOHfnpS/Wc+frizlSEOKGkd24YrCat4kkAgVAAlu34yCTZ2byycodDGzfiHvH9qFdkzpBlyUiFUQBkIAKQ86TH6/m/reWUq1KFe4a05sL+rdR/x6RBKMASDBLt+xj4oxMvlm/m2HdmnHnmF60rK/+PSKJSAGQII4UhPjH+yv4+/wV1KtZnYcvOJHRJ7RS/x6RBKYASADfrN/NxOmZLN26j9EntOLms3vQuG6NoMsSkYApACqxQ0cK+fPbS/nXR6tpVq8m/7w0jTN6NA+6LBGJEQqASmrhyh1MnpnJ2h0H+dnAVCaf2Y3kmmreJiL/nwKgktmbl8/dry/h35+vo23j2rzwi4Gc0lHN20Tk2xQAlcg7i7Zy4ytZbNt3mAlDOvD7M7pQK0ltHESkZAqASmDH/sPc+uoi5nyzia7N6/HYJWmc2KZB0GWJSIxTAMQxd2fON5u4ZU4O+w8X8PszuvDLoR1JqqY2DiJybAqAOLV5zyGmzsrm3SW5nNimAfeN60OX5tH6HztFJBEoAOJMKOT8+4t13P36EgpCIaaO6s7lg9tTVW0cRKSUFABxZPX2A0yekclnq3dySsfG3HNeH1Ib1w66LBGJUwqAOFBQGOKJj1fzwFvLSKpahXvO681P+7dRGwcROS4KgBi3ePNeJs3IJHPDHs7o3pw7zu1Fi/o1gy5LRCoBBUCMOlxQyN/nr+Qf81dQv1Z1/vazkxjVu6WO+kUkahQAMejLdbuYND2T5bn7GXNSa/54Vg8a1kkKuiwRqWQUADHk4JECHnhrGU98vJoWyTV5cnx/Tu/WLOiyRKSSUgDEiI9XbGfyzEzW7zzExSenMmlEN+qpeZuIlCMFQMD2HMrnrtcW81L6eto3qcNLE05mYIfGQZclIglAARCgt3K2MPWVbLbvP8xVpxU1b6tZXc3bRKRiRNw0xsyqmtlXZja3hLEhZvalmRWY2bhiY5eZ2fLw12VHLX/fzJaa2dfhr4S52L1t32F+/cKXTHg2g0Z1knjl14OZcmZ37fxFpEKV5gzgGmAxkFzC2DpgPHDd0QvNrBFwM5AGOJBhZnPcfVd4lYvcPb20Rccrd+eVrzdy66uLOHi4kOt+3IWrTutI9apq3iYiFS+iADCzFGAUcCdwbfFxd18TXi9UbGg48La77wyPvw2MAP5d9pLj08bdh7hxVhbvL91G39Si5m2dmql5m4gEJ9IzgIeAiUBp91itgfVHPd8QXvZfT5pZITADuMPdvfgbmNkEYAJAampqKTcfvFDIef6ztdzzxhJCDjef3YNLB7VT8zYRCdwxA8DMzgJy3T3DzIaW8v1L2sv9dyd/kbtvNLN6FAXAJcAz31rZfRowDSAtLe1bARHLVm3bz+QZWXy+Zic/6NSEu8/rTZtGat4mIrEhkjOAwcBoMxsJ1ASSzew5d784gtduAIYe9TwFeB/A3TeG/9xnZi8AAyghAOJRQWGIxz9czYPvLKNmtSrcN64P5/dLURsHEYkpxwwAd58CTAEInwFcF+HOH2AecJeZNQw//zEwxcyqAQ3cfbuZVQfOAt4pbfGxaNGmvUyc8Q3ZG/cyvGdzbj+nF82S1bxNRGJPmT8HYGa3AenuPsfM+gOzgIbA2WZ2q7v3dPedZnY78EX4ZbeFl9UB5oV3/lUp2vk/fnxTCVZefiF/e28Fj36wkga1k3jkor6c2btl0GWJiHwnK+G+a8xKS0vz9PTY+63RjLU7mTg9k5XbDjC2bwo3ndWdBrXVvE1EYoOZZbh7WvHl+iTwcThwuIA/zVvK0wvX0Kp+LZ6+YgCndWkadFkiIhFRAJTRgmXbmDIzi427D3HZoLZcP6IbdWvor1NE4of2WKW051VOqp4AAAXuSURBVGA+t7+2iOkZG+jQtA4vXz2I/u0aBV2WiEipKQBK4c3szdw0O4edB47wq6Ed+d0PO6t/j4jELQVABHL35XHz7BzeyN5Cj5bJPDm+P71a1w+6LBGR46IA+B7uzvSMDdzx2mIO5Rdy/fCuTBjSQc3bRKRSUAB8h/U7D3LDrCw+XL6dtLYNuWdsHzo1qxt0WSIiUaMAKCYUcp5ZuIb75i0F4NbRPbnk5LZUUfM2EalkFABHWZG7n8kzMklfu4shXZpy15hepDRU8zYRqZwUAEB+YYhpC1bx8DvLqZVUlQfOP4Hz+rZW8zYRqdQSPgCyN+5h4vRMFm3ey8jeLbh1dC+a1qsRdFkiIuUuYQMgL7+Qh99dzrQFq2hUJ4lHL+7LiF5q3iYiiSMhA+CLNTuZND2TVdsPcH6/FKaO6kH92tWDLktEpEIlVADsP1zAfW8u4ZmFa0lpWItnrxzAqZ3VvE1EElPCBMD7S3O5cVY2m/Yc4vLB7bjux12po+ZtIpLAEmIPOGVmFv/+fB2dmtVl+tWn0K9tw2O/SESkkkuIAGjXuDa/HdaJ3wzrRI1qat4mIgIJEgBXndYx6BJERGKOupqJiCQoBYCISIJSAIiIJCgFgIhIglIAiIgkKAWAiEiCUgCIiCQoBYCISIIydw+6hoiZ2TZgbRlf3gTYHsVy4oHmnBg058rveOfb1t2/1fkyrgLgeJhZurunBV1HRdKcE4PmXPmV13x1CUhEJEEpAEREElQiBcC0oAsIgOacGDTnyq9c5psw9wBEROT/SqQzABEROYoCQEQkQVW6ADCzEWa21MxWmNnkEsZrmNlL4fHPzKxdxVcZXRHM+VozW2RmmWb2rpm1DaLOaDrWnI9ab5yZuZnF9a8MRjJfM/tJ+PucY2YvVHSN0RbBz3Wqmc03s6/CP9sjg6gzmszsCTPLNbPs7xg3M/tL+O8k08z6HtcG3b3SfAFVgZVAByAJ+AboUWydXwGPhh9fALwUdN0VMOfTgdrhx79MhDmH16sHLAA+BdKCrrucv8edga+AhuHnzYKuuwLmPA34ZfhxD2BN0HVHYd5DgL5A9neMjwTeAAw4GfjseLZX2c4ABgAr3H2Vux8BXgTOKbbOOcDT4cfTgR+amVVgjdF2zDm7+3x3Pxh++imQUsE1Rlsk32eA24H7gLyKLK4cRDLfXwB/d/ddAO6eW8E1Rlskc3YgOfy4PrCpAusrF+6+ANj5PaucAzzjRT4FGphZy7Jur7IFQGtg/VHPN4SXlbiOuxcAe4DGFVJd+Yhkzke7kqIjiHh2zDmb2UlAG3efW5GFlZNIvsddgC5m9rGZfWpmIyqsuvIRyZxvAS42sw3A68BvK6a0QJX23/v3qmz/KXxJR/LFf881knXiScTzMbOLgTTgtHKtqPx975zNrArwIDC+ogoqZ5F8j6tRdBloKEVneB+aWS93313OtZWXSOZ8IfCUuz9gZoOAZ8NzDpV/eYGJ6v6rsp0BbADaHPU8hW+fFv7vOmZWjaJTx+875Yp1kcwZMzsDuBEY7e6HK6i28nKsOdcDegHvm9kaiq6VzonjG8GR/lzPdvd8d18NLKUoEOJVJHO+EvgPgLsvBGpS1DStMovo33ukKlsAfAF0NrP2ZpZE0U3eOcXWmQNcFn48DnjPw3dX4tQx5xy+HPIYRTv/eL82DMeYs7vvcfcm7t7O3dtRdN9jtLunB1PucYvk5/oVim72Y2ZNKLoktKpCq4yuSOa8DvghgJl1pygAtlVolRVvDnBp+LeBTgb2uPvmsr5ZpboE5O4FZvYbYB5Fv0XwhLvnmNltQLq7zwH+RdGp4gqKjvwvCK7i4xfhnP8E1AVeDt/vXufuowMr+jhFOOdKI8L5zgN+bGaLgELgenffEVzVxyfCOf8BeNzMfk/RZZDxcX4wh5n9m6LLeE3C9zZuBqoDuPujFN3rGAmsAA4Clx/X9uL870tERMqosl0CEhGRCCkAREQSlAJARCRBKQBERBKUAkBEJEEpAEREEpQCQEQkQf0/xX8YSukQRQ8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(N.sum(axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimate movement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = Akagi(N, d, K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.lamda = 1e-3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/covid/lib/python3.6/site-packages/numba/np/ufunc/parallel.py:355: NumbaWarning: The TBB threading layer requires TBB version 2019.5 or later i.e., TBB_INTERFACE_VERSION >= 11005. Found TBB_INTERFACE_VERSION = 10005. The TBB threading layer is disabled.\n",
      "  warnings.warn(problem)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L  -11516978500.0 -17200931235.105717\n",
      "step #  0 , L =  -17200931235.105717\n",
      "beta  0.1\n",
      "pi  [0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02\n",
      " 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02]\n",
      "L  -11516978500.0 -17200931235.105717\n",
      "L  -8473499778.938085 -14199828460.118093\n",
      "L  -2140474910.8437152 -8044522380.397306\n",
      "L  -2089432034.342115 -8000100785.80077\n",
      "L  -1864672499.2419987 -7789447069.662265\n",
      "L  -1507591583.165478 -7430554060.344786\n",
      "L  -891673638.1543653 -6773466825.376482\n",
      "L  -377855970.50362873 -6180410134.941915\n",
      "L  -175060506.08331066 -5921281003.530142\n",
      "L  -107228611.64134148 -5834782793.284738\n",
      "L  -67610704.63863763 -5795044343.42051\n",
      "L  -48584166.672355354 -5777124600.41204\n",
      "L  -40175369.04467631 -5768978920.716447\n",
      "L  -34626027.94237438 -5762946414.555107\n",
      "L  -31476811.30787311 -5759030859.551863\n",
      "L  -28507796.161398772 -5754257095.78191\n",
      "L  -28231852.8124086 -5752942949.169406\n",
      "L  -29541542.68464518 -5752593278.361451\n",
      "L  -29984949.21142276 -5752810785.73913\n",
      "L  -29595383.762074787 -5752583620.392287\n",
      "M done\n",
      "pi done\n",
      "beta done\n",
      "L  -29595383.762074787 -5722227818.650305\n",
      "step #  1 , L =  -5722227818.650305\n",
      "beta  0.06712901997863827\n",
      "pi  [5.67741126e-03 3.83212510e-02 5.21597260e-02 1.34169815e-03\n",
      " 8.23747941e-02 7.38773994e-02 1.82975173e-02 9.72987098e-02\n",
      " 9.39488169e-02 7.93480246e-02 1.11557709e-01 4.93016036e-03\n",
      " 1.22497463e-01 3.71459944e-02 1.25014553e-01 5.40271242e-02\n",
      " 5.59575341e-02 9.14220482e-02 1.31245082e-01 6.64566114e-02\n",
      " 6.85059620e-02 1.16874998e-01 1.20957052e-01 4.38104531e-02\n",
      " 1.31078607e-01 1.31452157e-01 1.37386064e-01 2.56086848e-03\n",
      " 2.74630916e-02 6.37587060e-03 7.25137190e-02 1.02557324e-01\n",
      " 7.21281825e-02 7.92640764e-02 4.28184446e-02 9.70897423e-02\n",
      " 6.41981525e-03 6.09523556e-02 2.45222556e-02 6.47707001e-03\n",
      " 4.03850514e-03 2.36781787e-02 3.06656298e-02 6.34820585e-02\n",
      " 6.93446132e-02 3.59116984e-02 2.26374932e-02 2.58837932e-03\n",
      " 2.62505322e-03 9.60427166e-02 6.64869156e-02 7.66911361e-02\n",
      " 1.14679287e-01 7.46194136e-02 3.11124659e-03 3.29157506e-02\n",
      " 7.27419091e-03 1.04176593e-01 2.64802058e-02 2.34977969e-03\n",
      " 3.20417519e-03 4.50582758e-02 6.70023724e-02 4.20329906e-03\n",
      " 6.91896740e-02 7.60847771e-02 3.90357562e-03 8.62030595e-02\n",
      " 6.49921229e-03 4.29442744e-03 7.19077455e-02 2.71257618e-03\n",
      " 1.24788731e-01 1.00201766e-01 9.11376513e-02 2.76947057e-02\n",
      " 8.67301866e-02 1.23803814e-01 6.90723029e-03 1.04732756e-01\n",
      " 3.23806328e-03 8.02235598e-02 9.39698652e-02 1.71270020e-01\n",
      " 6.16030532e-02 2.58802935e-03 2.71277343e-03 1.13521842e-01\n",
      " 6.53100582e-02 8.80398337e-02 7.53719940e-02 4.50408390e-03\n",
      " 2.13506391e-01 1.90087767e-03 3.75118478e-03 9.02028130e-02\n",
      " 1.78969826e-01 4.62032406e-03 1.20828920e-01 1.08014459e-01\n",
      " 6.75053291e-03 8.08800258e-02 1.35718707e-04 8.38960437e-02\n",
      " 6.97964035e-02 1.02166576e-01 4.16885113e-03 1.34387273e-03\n",
      " 2.57492373e-02 1.21423426e-01 6.45705870e-03 7.04262043e-02\n",
      " 4.81540648e-02 1.50282170e-01 4.87098549e-02 5.15315801e-02\n",
      " 9.22014496e-02 2.06680812e-03 8.58715885e-02 6.55879168e-02\n",
      " 4.43415152e-02 5.07960325e-02 8.95918636e-02 3.27345535e-02\n",
      " 3.27306796e-03 5.17827763e-02 7.88927654e-02 7.45284137e-02\n",
      " 7.93947703e-03 3.98263837e-02 9.07159905e-02 1.21203864e-01\n",
      " 9.10060425e-02 7.88313824e-02 5.40601803e-02 5.34408853e-02\n",
      " 1.34060328e-01 8.51695672e-02 7.28916580e-02 7.45683708e-02\n",
      " 8.76591376e-02 8.70245476e-02 1.06723728e-01 3.79020976e-02\n",
      " 8.19203681e-02 8.49708692e-02 1.06542372e-01 5.47695813e-02\n",
      " 8.78394689e-02 1.11518666e-01 4.10977190e-02 7.39733619e-03\n",
      " 9.39845067e-02 1.44060542e-01 8.28133721e-03 3.26964176e-02\n",
      " 5.59937769e-02 1.02767987e-01 9.05136745e-02 1.03255255e-01\n",
      " 5.83132245e-02 8.32930833e-02 9.53000509e-02 5.84547728e-02\n",
      " 8.21367456e-02 1.15569099e-01 1.41886351e-02]\n",
      "L  -29595383.762074787 -5722227818.650305\n",
      "L  -30132153.900643434 -5722135306.621993\n",
      "L  -31987555.261505242 -5721793379.1140375\n",
      "L  -31987155.994708143 -5721551048.910359\n",
      "L  -32068441.450851556 -5721024798.9147415\n",
      "L  -31827488.350007504 -5720124691.200269\n",
      "L  -31334673.970704194 -5719123879.915715\n",
      "L  -30976034.367043156 -5718246183.391383\n",
      "L  -30917123.345097795 -5717735972.984603\n",
      "L  -30848615.270055264 -5716699087.514412\n",
      "L  -30182462.498772744 -5715819225.394375\n",
      "L  -31758313.824637793 -5715508227.263897\n",
      "L  -30124298.2211542 -5714676191.3873415\n",
      "L  -30237156.934457663 -5714405392.539141\n",
      "L  -30363402.60839596 -5714298113.001392\n",
      "L  -28884602.82948749 -5713725686.791171\n",
      "L  -29141241.672577582 -5713686978.224254\n",
      "M done\n",
      "pi done\n",
      "beta done\n",
      "L  -29141241.672577582 -5711316285.553089\n",
      "step #  2 , L =  -5711316285.553089\n",
      "beta  0.055742346112640756\n",
      "pi  [4.70794561e-03 3.05479635e-02 4.81041783e-02 9.98654210e-04\n",
      " 7.67422322e-02 6.78808966e-02 1.41319101e-02 8.71359941e-02\n",
      " 8.89192932e-02 7.31276369e-02 1.04719229e-01 3.75099342e-03\n",
      " 1.13347015e-01 3.30093061e-02 1.21247832e-01 5.40163981e-02\n",
      " 4.89298090e-02 8.32133605e-02 1.10298020e-01 6.11485176e-02\n",
      " 6.29759807e-02 1.04912534e-01 1.14555698e-01 4.18228155e-02\n",
      " 1.21779147e-01 1.23264636e-01 1.26868509e-01 2.03601505e-03\n",
      " 2.33810956e-02 4.43638605e-03 6.45576555e-02 9.78320547e-02\n",
      " 6.89743889e-02 8.26294752e-02 3.49691594e-02 9.65039869e-02\n",
      " 4.96039047e-03 5.56743355e-02 2.16108232e-02 4.94810774e-03\n",
      " 3.19590258e-03 2.07406747e-02 2.61214686e-02 5.71112063e-02\n",
      " 6.40550257e-02 3.14770474e-02 2.06744110e-02 2.06802592e-03\n",
      " 2.14542687e-03 1.05545757e-01 6.46248510e-02 7.12450684e-02\n",
      " 1.07115399e-01 6.95169492e-02 2.49672416e-03 2.65470166e-02\n",
      " 6.03381980e-03 9.76059659e-02 2.20941281e-02 1.79823629e-03\n",
      " 2.59674431e-03 3.60176874e-02 6.19795288e-02 3.32335372e-03\n",
      " 5.63654394e-02 7.01677446e-02 3.26032584e-03 8.09656866e-02\n",
      " 5.18995900e-03 3.40070240e-03 6.69678471e-02 2.12332840e-03\n",
      " 1.11731636e-01 9.20309231e-02 8.23074474e-02 2.46923400e-02\n",
      " 8.25613052e-02 1.06840267e-01 5.53051891e-03 9.77197689e-02\n",
      " 2.59265924e-03 7.55538063e-02 8.67349076e-02 1.54393713e-01\n",
      " 5.62330117e-02 2.04581032e-03 2.18322016e-03 1.21932738e-01\n",
      " 5.23074999e-02 9.60685266e-02 7.26271004e-02 3.53192580e-03\n",
      " 2.46361750e-01 1.77741894e-03 3.19794423e-03 8.42539223e-02\n",
      " 1.63355688e-01 3.72154774e-03 1.15792540e-01 1.09668707e-01\n",
      " 5.48021773e-03 7.07026041e-02 2.18562554e-04 7.79823797e-02\n",
      " 6.28176118e-02 9.62855042e-02 3.34979170e-03 1.30924800e-03\n",
      " 2.10150671e-02 1.15637770e-01 5.33915913e-03 6.15813167e-02\n",
      " 4.44929976e-02 1.51062608e-01 3.96854733e-02 4.74414992e-02\n",
      " 9.07189854e-02 1.70066899e-03 8.94002258e-02 6.19437197e-02\n",
      " 3.89627056e-02 4.48218499e-02 8.57449246e-02 2.54259753e-02\n",
      " 2.63606248e-03 4.50997800e-02 7.09399466e-02 7.90032315e-02\n",
      " 6.49592832e-03 3.55570216e-02 8.05399961e-02 1.09503445e-01\n",
      " 8.05558564e-02 7.24351464e-02 5.47184129e-02 4.71950541e-02\n",
      " 1.43895305e-01 7.76526922e-02 6.76694205e-02 7.11507050e-02\n",
      " 8.09779020e-02 7.97003752e-02 9.84087418e-02 3.10266449e-02\n",
      " 7.57929769e-02 7.64430432e-02 9.35983255e-02 4.65192984e-02\n",
      " 8.17005134e-02 1.09255613e-01 3.65024963e-02 6.14958123e-03\n",
      " 8.69350290e-02 1.53070653e-01 6.86122753e-03 2.57141891e-02\n",
      " 5.17181479e-02 9.60811374e-02 8.10382310e-02 9.63196176e-02\n",
      " 5.16209002e-02 7.42096622e-02 8.86702726e-02 5.12803551e-02\n",
      " 8.99497105e-02 1.17450417e-01 1.09256202e-02]\n",
      "L  -29141241.672577582 -5711316285.553089\n",
      "L  -29448396.497190867 -5711343290.172014\n",
      "L  -29182764.086130787 -5711319242.311115\n",
      "L  -29147114.038275596 -5711316660.568181\n",
      "L  -29142091.023353066 -5711316334.676899\n",
      "L  -29141366.93303004 -5711316292.076255\n",
      "L  -29141260.503474932 -5711316286.426983\n",
      "L  -29141244.55879497 -5711316285.670848\n",
      "L  -29141242.12369929 -5711316285.568992\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L  -29141241.744503103 -5711316285.55523\n",
      "L  -29141241.684277583 -5711316285.553374\n",
      "L  -29141241.674519982 -5711316285.553126\n",
      "L  -29141241.672906354 -5711316285.553093\n",
      "L  -29141241.672635652 -5711316285.55309\n",
      "L  -29141241.67258743 -5711316285.55309\n",
      "L  -29141241.672578674 -5711316285.553089\n",
      "L  -29141241.672579914 -5711316285.553089\n",
      "L  -29141241.672583677 -5711316285.553089\n",
      "L  -29141241.67258613 -5711316285.55309\n",
      "L  -29141241.67258384 -5711316285.553089\n",
      "L  -29141241.672585413 -5711316285.553089\n",
      "Error minimizing M b'ABNORMAL_TERMINATION_IN_LNSRCH'\n",
      "M done\n",
      "pi done\n",
      "beta done\n",
      "L  -29141241.672577582 -5711316285.514339\n",
      "CPU times: user 2min 58s, sys: 2.99 s, total: 3min 1s\n",
      "Wall time: 1min 48s\n"
     ]
    }
   ],
   "source": [
    "%time result = a.exact_inference(1e-4,True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2382076,       0,       0, ...,       0,       0,       0],\n",
       "       [    529, 1561044,       3, ...,       0,       0,      97],\n",
       "       [    924,     195, 2977775, ...,       0,       0,     215],\n",
       "       ...,\n",
       "       [   3216,    1419,     222, ..., 3485691,       0,    1466],\n",
       "       [   1830,     921,      30, ...,       0, 1874823,     962],\n",
       "       [    371,      94,      14, ...,       0,       0, 4161394]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.rint(a.M[0]).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are there approximately the right number of people in the end?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([410154000, 413233000])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N.sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2393343, 1610233, 3128256, 8216389, 2193455, 2044757, 2428033,\n",
       "        1138554, 2137396, 1977739, 1711580, 2232090, 1309123, 2904764,\n",
       "        1749410, 4439540, 1681872, 1442243,  586199, 2307676, 2216707,\n",
       "         971431, 1652628, 4888401, 1291020, 1434690, 1118360, 4942985,\n",
       "        3011256, 2706343, 1401653, 2023456, 2942911, 3354712, 1611506,\n",
       "        2403116, 1793516, 2287775, 4150479, 2320390, 2979346, 3863091,\n",
       "        2694301, 1966365, 2328520, 2707561, 6421125, 5013989, 3949955,\n",
       "        3307334, 3352768, 2259435, 1599655, 2392531, 4144109, 1969397,\n",
       "        1846465, 1804518, 2931985, 5021361, 4006170, 1298516, 2443473,\n",
       "        3075263,  936072, 2072714, 3695731, 2234357, 1958727, 3013270,\n",
       "        2493388, 4121580,  851923, 1439285, 1366787, 4087319, 2385273,\n",
       "         656210, 1710797, 1723721, 3769030, 2448318, 1660945,  678782,\n",
       "        2203824, 4286448, 4849017, 2561756,  790852, 3519329, 2885914,\n",
       "        2378953, 1813268, 5362721, 3625915, 2035447,  734891, 2763425,\n",
       "        1741545, 2273141, 1740852, 1135128, 8773767, 2073503, 1766428,\n",
       "        1913464, 3143297, 6902756, 2956438, 1689619, 2184078, 1277453,\n",
       "        3405213, 1581486, 1400417, 2965578, 2473125, 5088762, 3077838,\n",
       "        3054184, 2204776, 1764026, 2359207, 1273478, 3883063, 1646462,\n",
       "        1524825, 3768484, 1705325, 2834631, 1079667, 1014113,  953778,\n",
       "        1822061, 4579483, 1676045, 2156889, 1557119, 2408275, 2815982,\n",
       "        1814766, 1591185, 1420086, 1985966, 1978676, 1390014,  836018,\n",
       "        1456477, 1988551, 2016270, 2501120, 1805497, 1715880, 1936071,\n",
       "        1610616, 1529243, 2926383, 1807500, 1200593, 1755638, 1709355,\n",
       "        1227782, 1863527, 1633072, 3830218, 2124326, 4207362]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.M.sum(axis=(2)).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2343,   234,  -743,  -611,  -544,  -242,    34,   554,  -604,\n",
       "         -261,  -420,  1090,   124,  -236,  -590, -1459,  -127,   244,\n",
       "         1199,  -323,  -293,   431,  -372, -1598,    21,  -309,   360,\n",
       "          986,  -744,   343,   654,  -544, -1088, -1287,   507,  -883,\n",
       "         1516,  -224,  -520,  1390,  1347,    91,  -698,  -634,  -479,\n",
       "         -438,  -874,   989,   -44, -1666, -1232,  -565,  -345,  -468,\n",
       "         1110,   398,  2465,  -482,   -14,   361,  1171,   517,  -527,\n",
       "         1264,  1072,  -285,  1732,  -642,  1727,  1270,  -611,   580,\n",
       "          924,   285,  -212,  -680,  -726,  1210,  1798,  -279,  1031,\n",
       "         -682,   -55,   782,  -176,   448,  1017, -1244,   853, -1671,\n",
       "        -1086,   954, -1731,  -278,  1915,  -553,   892,  1426,  -454,\n",
       "         -859,  1853,   129, -3233,  -496,  -572,  -535,  1297, -1244,\n",
       "         -562,  -380,  2078,   453,  -787,  -514,   418,  -422,  -875,\n",
       "         -237, -1161,  -815,  -223,  1027,  -792,   478,  1064,   462,\n",
       "         -174, -1516,  2326,  -369,   667,   114,  1778,    61, -1516,\n",
       "         1045, -1111,   119,  -725, -1018,  -234,   185,    87,   -34,\n",
       "         -323,    14,  1018,   478,  -448,  -729,   120,  2498,  -120,\n",
       "         -928,  2617,   244,  -617,  -500,   593,  -361,   355,   782,\n",
       "         -473,    72, -1782,  -674, -1638]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.rint((a.M.sum(axis=2) - N[:-1])).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([410154000, 413233000])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N.sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "x and y must have same first dimension, but have shapes (3,) and (2,)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-24-8f7e0d8c1201>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m plt.plot(\n\u001b[1;32m      2\u001b[0m     \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m3\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----> 3\u001b[0;31m     \u001b[0mN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\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      4\u001b[0m )\n\u001b[1;32m      5\u001b[0m plt.plot(\n",
      "\u001b[0;32m/anaconda3/envs/covid/lib/python3.6/site-packages/matplotlib/pyplot.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(scalex, scaley, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m   2794\u001b[0m     return gca().plot(\n\u001b[1;32m   2795\u001b[0m         *args, scalex=scalex, scaley=scaley, **({\"data\": data} if data\n\u001b[0;32m-> 2796\u001b[0;31m         is not None else {}), **kwargs)\n\u001b[0m\u001b[1;32m   2797\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2798\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/envs/covid/lib/python3.6/site-packages/matplotlib/axes/_axes.py\u001b[0m in \u001b[0;36mplot\u001b[0;34m(self, scalex, scaley, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1663\u001b[0m         \"\"\"\n\u001b[1;32m   1664\u001b[0m         \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcbook\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormalize_kwargs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmlines\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mLine2D\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_alias_map\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1665\u001b[0;31m         \u001b[0mlines\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_lines\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdata\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[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1666\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mline\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlines\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1667\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_line\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/envs/covid/lib/python3.6/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    223\u001b[0m                 \u001b[0mthis\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0margs\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    224\u001b[0m                 \u001b[0margs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\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--> 225\u001b[0;31m             \u001b[0;32myield\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_plot_args\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mthis\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    226\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    227\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget_next_color\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\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/anaconda3/envs/covid/lib/python3.6/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36m_plot_args\u001b[0;34m(self, tup, kwargs)\u001b[0m\n\u001b[1;32m    389\u001b[0m             \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mindex_of\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtup\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\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    390\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 391\u001b[0;31m         \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_xy_from_xy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\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    392\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    393\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcommand\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'plot'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/envs/covid/lib/python3.6/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36m_xy_from_xy\u001b[0;34m(self, x, y)\u001b[0m\n\u001b[1;32m    268\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    269\u001b[0m             raise ValueError(\"x and y must have same first dimension, but \"\n\u001b[0;32m--> 270\u001b[0;31m                              \"have shapes {} and {}\".format(x.shape, y.shape))\n\u001b[0m\u001b[1;32m    271\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\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[1;32m    272\u001b[0m             raise ValueError(\"x and y can be no greater than 2-D, but have \"\n",
      "\u001b[0;31mValueError\u001b[0m: x and y must have same first dimension, but have shapes (3,) and (2,)"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAANgElEQVR4nO3ccYjfd33H8efLxE6mtY7lBEmi7Vi6Gsqg7ug6hFnRjbR/JP8USaC4SmnArQ5mETocKvWvKUMQsmm2iVPQWv1DD4nkD1fpECO50lmalMAtOnNE6Fm7/lO0Znvvj99P77hcct/e/e4u3vv5gMDv+/t9fr9758PdM798f/f7paqQJG1/r9rqASRJm8PgS1ITBl+SmjD4ktSEwZekJgy+JDWxavCTfC7Jc0meucLtSfLpJHNJnk7ytsmPKUlaryHP8D8PHLjK7XcB+8Z/jgL/tP6xJEmTtmrwq+oJ4GdXWXII+EKNnALekORNkxpQkjQZOyfwGLuBC0uO58fX/WT5wiRHGf0vgNe+9rV/dMstt0zgy0tSH08++eRPq2pqLfedRPCzwnUrfl5DVR0HjgNMT0/X7OzsBL68JPWR5L/Xet9J/JbOPLB3yfEe4OIEHleSNEGTCP4M8N7xb+vcAbxYVZedzpEkba1VT+kk+TJwJ7AryTzwUeDVAFX1GeAEcDcwB7wEvG+jhpUkrd2qwa+qI6vcXsBfTWwiSdKG8J22ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNTEo+EkOJDmXZC7Jwyvc/uYkjyd5KsnTSe6e/KiSpPVYNfhJdgDHgLuA/cCRJPuXLfs74LGqug04DPzjpAeVJK3PkGf4twNzVXW+ql4GHgUOLVtTwOvHl28ALk5uREnSJAwJ/m7gwpLj+fF1S30MuDfJPHAC+MBKD5TkaJLZJLMLCwtrGFeStFZDgp8Vrqtlx0eAz1fVHuBu4ItJLnvsqjpeVdNVNT01NfXKp5UkrdmQ4M8De5cc7+HyUzb3A48BVNX3gNcAuyYxoCRpMoYE/zSwL8lNSa5j9KLszLI1PwbeBZDkrYyC7zkbSbqGrBr8qroEPAicBJ5l9Ns4Z5I8kuTgeNlDwANJfgB8Gbivqpaf9pEkbaGdQxZV1QlGL8Yuve4jSy6fBd4+2dEkSZPkO20lqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0MCn6SA0nOJZlL8vAV1rwnydkkZ5J8abJjSpLWa+dqC5LsAI4BfwbMA6eTzFTV2SVr9gF/C7y9ql5I8saNGliStDZDnuHfDsxV1fmqehl4FDi0bM0DwLGqegGgqp6b7JiSpPUaEvzdwIUlx/Pj65a6Gbg5yXeTnEpyYKUHSnI0yWyS2YWFhbVNLElakyHBzwrX1bLjncA+4E7gCPAvSd5w2Z2qjlfVdFVNT01NvdJZJUnrMCT488DeJcd7gIsrrPlGVf2yqn4InGP0D4Ak6RoxJPingX1JbkpyHXAYmFm25uvAOwGS7GJ0iuf8JAeVJK3PqsGvqkvAg8BJ4Fngsao6k+SRJAfHy04Czyc5CzwOfKiqnt+ooSVJr1yqlp+O3xzT09M1Ozu7JV9bkn5TJXmyqqbXcl/faStJTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITg4Kf5ECSc0nmkjx8lXX3JKkk05MbUZI0CasGP8kO4BhwF7AfOJJk/wrrrgf+Gvj+pIeUJK3fkGf4twNzVXW+ql4GHgUOrbDu48AngJ9PcD5J0oQMCf5u4MKS4/nxdb+W5DZgb1V982oPlORoktkkswsLC694WEnS2g0Jfla4rn59Y/Iq4FPAQ6s9UFUdr6rpqpqempoaPqUkad2GBH8e2LvkeA9wccnx9cCtwHeS/Ai4A5jxhVtJurYMCf5pYF+Sm5JcBxwGZn51Y1W9WFW7qurGqroROAUcrKrZDZlYkrQmqwa/qi4BDwIngWeBx6rqTJJHkhzc6AElSZOxc8iiqjoBnFh23UeusPbO9Y8lSZo032krSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWpiUPCTHEhyLslckodXuP2DSc4meTrJt5O8ZfKjSpLWY9XgJ9kBHAPuAvYDR5LsX7bsKWC6qv4Q+BrwiUkPKklanyHP8G8H5qrqfFW9DDwKHFq6oKoer6qXxoengD2THVOStF5Dgr8buLDkeH583ZXcD3xrpRuSHE0ym2R2YWFh+JSSpHUbEvyscF2tuDC5F5gGPrnS7VV1vKqmq2p6ampq+JSSpHXbOWDNPLB3yfEe4OLyRUneDXwYeEdV/WIy40mSJmXIM/zTwL4kNyW5DjgMzCxdkOQ24LPAwap6bvJjSpLWa9XgV9Ul4EHgJPAs8FhVnUnySJKD42WfBF4HfDXJfyaZucLDSZK2yJBTOlTVCeDEsus+suTyuyc8lyRpwnynrSQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0MCn6SA0nOJZlL8vAKt/9Wkq+Mb/9+khsnPagkaX1WDX6SHcAx4C5gP3Akyf5ly+4HXqiq3wc+Bfz9pAeVJK3PkGf4twNzVXW+ql4GHgUOLVtzCPi38eWvAe9KksmNKUlar50D1uwGLiw5ngf++EprqupSkheB3wV+unRRkqPA0fHhL5I8s5aht6FdLNurxtyLRe7FIvdi0R+s9Y5Dgr/SM/Vawxqq6jhwHCDJbFVND/j62557sci9WOReLHIvFiWZXet9h5zSmQf2LjneA1y80pokO4EbgJ+tdShJ0uQNCf5pYF+Sm5JcBxwGZpatmQH+Ynz5HuDfq+qyZ/iSpK2z6imd8Tn5B4GTwA7gc1V1JskjwGxVzQD/CnwxyRyjZ/aHB3zt4+uYe7txLxa5F4vci0XuxaI170V8Ii5JPfhOW0lqwuBLUhMbHnw/lmHRgL34YJKzSZ5O8u0kb9mKOTfDanuxZN09SSrJtv2VvCF7keQ94++NM0m+tNkzbpYBPyNvTvJ4kqfGPyd3b8WcGy3J55I8d6X3KmXk0+N9ejrJ2wY9cFVt2B9GL/L+F/B7wHXAD4D9y9b8JfCZ8eXDwFc2cqat+jNwL94J/Pb48vs778V43fXAE8ApYHqr597C74t9wFPA74yP37jVc2/hXhwH3j++vB/40VbPvUF78afA24BnrnD73cC3GL0H6g7g+0Med6Of4fuxDItW3YuqeryqXhofnmL0noftaMj3BcDHgU8AP9/M4TbZkL14ADhWVS8AVNVzmzzjZhmyFwW8fnz5Bi5/T9C2UFVPcPX3Mh0CvlAjp4A3JHnTao+70cFf6WMZdl9pTVVdAn71sQzbzZC9WOp+Rv+Cb0er7kWS24C9VfXNzRxsCwz5vrgZuDnJd5OcSnJg06bbXEP24mPAvUnmgRPABzZntGvOK+0JMOyjFdZjYh/LsA0M/nsmuReYBt6xoRNtnavuRZJXMfrU1fs2a6AtNOT7Yiej0zp3Mvpf338kubWq/meDZ9tsQ/biCPD5qvqHJH/C6P0/t1bV/238eNeUNXVzo5/h+7EMi4bsBUneDXwYOFhVv9ik2TbbantxPXAr8J0kP2J0jnJmm75wO/Rn5BtV9cuq+iFwjtE/ANvNkL24H3gMoKq+B7yG0QerdTOoJ8ttdPD9WIZFq+7F+DTGZxnFfruep4VV9qKqXqyqXVV1Y1XdyOj1jINVteYPjbqGDfkZ+TqjF/RJsovRKZ7zmzrl5hiyFz8G3gWQ5K2Mgr+wqVNeG2aA945/W+cO4MWq+slqd9rQUzq1cR/L8Btn4F58Engd8NXx69Y/rqqDWzb0Bhm4Fy0M3IuTwJ8nOQv8L/Chqnp+66beGAP34iHgn5P8DaNTGPdtxyeISb7M6BTervHrFR8FXg1QVZ9h9PrF3cAc8BLwvkGPuw33SpK0At9pK0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDXx/4aZaro1YsjCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(\n",
    "    np.arange(3),\n",
    "    N.sum(axis=1)\n",
    ")\n",
    "plt.plot(\n",
    "    np.arange(0.5, 2.5),\n",
    "    a.M.sum(axis=(1,2))\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cycler import cycler\n",
    "import matplotlib.pyplot as plt\n",
    "import itertools\n",
    "\n",
    "color_list = ['b', 'orange', 'r', 'green', 'k', 'gray']\n",
    "color_cycle = []\n",
    "for i in range(len(color_list)):\n",
    "    color_cycle.append(color_list[i])\n",
    "    color_cycle.append(color_list[i])\n",
    "\n",
    "num_regions = 6\n",
    "plt.rc('axes', prop_cycle=(cycler('color', color_cycle)))\n",
    "\n",
    "for i in range(num_regions):\n",
    "    plt.plot(\n",
    "        np.arange(3),\n",
    "        N[:, i],\n",
    "        linestyle='-',\n",
    "    )\n",
    "    plt.plot(\n",
    "        np.arange(0.5, 2.5),\n",
    "        a.M[:, i].sum(axis=1),\n",
    "        linestyle=':',\n",
    "    )"
   ]
  },
  {
   "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}