real-data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "centroid_data_fname = \"data/area/statistical-area-2-2018-centroid-true.csv\"\n",
    "hierarchy_data_fname = \"data/area/statistical-area-2-higher-geographies-2018-generalised.csv\"\n",
    "areas = distances.AreaSubset(centroid_data_fname, hierarchy_data_fname)\n",
    "areas.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter by distance from some SA2"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "center = \"Christchurch Central\"\n",
    "D = 5000 # meters but K is in km to avoid numerical issues\n",
    "areas.filter_radius(center, D)\n",
    "areas.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter by territorial authority."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "areas.filter_rc_name(['Auckland Region', 'Waikato Region'])\n",
    "# areas.filter_ta_name(['Nelson City'])\n",
    "areas.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "telco_data_fname = \"data/telco/pop_data_2020-04-01.dat\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "telco = load_csv.load_telco_data(telco_data_fname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter to only data from times of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "date_first = datetime(2020, 2, 18, 7)\n",
    "date_last = datetime(2020, 2, 18, 9)\n",
    "\n",
    "telco_subset_time = telco[\n",
    "        (telco[\"time\"] == date_first) | (telco[\"time\"] == date_last)\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find regions with too few people and remove them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count_threshold = 10\n",
    "\n",
    "current_codes = set(areas.sa2_codes())\n",
    "# less than theshold at some time\n",
    "disallowed_codes = set(telco_subset_time[telco_subset_time['count'] <= count_threshold].sa2_code.to_list())\n",
    "allowed_codes = set(telco_subset_time.sa2_code.to_list()) - disallowed_codes\n",
    "\n",
    "sa2_codes_to_remove = current_codes - allowed_codes\n",
    "\n",
    "areas.remove_sa2(sa2_codes_to_remove)\n",
    "areas.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "areas.sa2_names()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter to only telco data from the regions of interest and put SA2 code in columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "telco_subset = pd.pivot_table(\n",
    "    telco_subset_time[\n",
    "        telco_subset_time[\"sa2_code\"].isin(areas.sa2_codes())\n",
    "    ],\n",
    "    index=\"time\",\n",
    "    columns=\"sa2_code\",\n",
    "    values=\"count\",\n",
    ")\n",
    "\n",
    "telco_subset = telco_subset.reindex(areas.sa2_codes(), axis='columns')\n",
    "\n",
    "N = telco_subset.to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert telco_subset.columns.to_list() == areas.sa2_codes()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = areas.distance_table(units='km')\n",
    "d.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.all(np.isfinite(N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.count_nonzero(N) == N.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 80"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(N.sum(axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimate movement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scale = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = Akagi(N * scale, d, K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.lamda = 1e1 / scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.M_jitter(10/N.max()/scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check that all regions have at least one neighbour\n",
    "assert np.all(a.gamma_exc.sum(axis=1) > 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%time result = a.exact_inference(1e-3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_est = a.M / scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(a.save_options.output_dir + \"/sa2_codes\", \"w\") as f:\n",
    "    for line in areas.sa2_codes():\n",
    "        f.write(line)\n",
    "        f.write(\", \")\n",
    "    f.write(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a.beta_bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check if beta bounds are saturated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.isclose(a.beta, a.beta_bounds[0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.isclose(a.beta, a.beta_bounds[0][1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are there approximately the right number of people in the end?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N.sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_est.sum(axis=(1, 2)).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.rint((M_est.sum(axis=2) - N[:-1])).astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does `M` look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "np.rint(a.M[0]).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(np.where(d < K, d, np.nan))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(np.where(d < K, np.log(M_est[0]), np.nan))\n",
    "plt.colorbar()"
   ]
  }
 ],
 "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
}