No commit message

No commit message

No commit message
662c083c · Mark Penney · 41f8e77b · 662c083c
Commit 662c083c authored 4 years ago by Mark Penney
--- a/network-data/Prem/Canadian-Avgs.ipynb
+++ b/network-data/Prem/Canadian-Avgs.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib notebook\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import pandas as pd\n",
+    "from scipy import optimize\n",
+    "import math\n",
+    "import random\n",
+    "import csv\n",
+    "import itertools\n",
+    "from itertools import product\n",
+    "def Flatten(A):\n",
+    "    return list(itertools.chain.from_iterable(A))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Average mixing numbers for Canadian population\n",
+    "\n",
+    "To adjust the POLYMOD results to fit Canadian demographic structures, we take the distributions extracted from POLYMOD and rescale them so that their averages match those determined by Prem."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Population statistics from Statcan 2019\n",
+    "N = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]\n",
+    "N[0] = 1943175\n",
+    "N[1] = 2039352\n",
+    "N[2] = 2031762\n",
+    "N[3] = 2114635\n",
+    "N[4] = 2476698\n",
+    "N[5] = 2625474\n",
+    "N[6] = 2603938\n",
+    "N[7] = 2580021\n",
+    "N[8] = 2421009\n",
+    "N[9] = 2396406\n",
+    "N[10] = 2502667\n",
+    "N[11] = 2749626\n",
+    "N[12] = 2511888\n",
+    "N[13] = 2096607\n",
+    "N[14] = 1706760\n",
+    "N[15] = 1164277+786704+510828+242554+74086+10795 #75+\n",
+    "Ntot = np.sum(N)\n",
+    "\n",
+    "# Load the Prem files\n",
+    "Malldf = pd.read_csv(\"Prem_Canada.csv\")\n",
+    "Mall = Malldf.values.tolist()\n",
+    "Mhomedf = pd.read_csv(\"Prem_home.csv\")\n",
+    "Mhome = Mhomedf.values.tolist()\n",
+    "Mschooldf = pd.read_csv(\"Prem_school.csv\")\n",
+    "Mschool = Mschooldf.values.tolist()\n",
+    "Mworkdf = pd.read_csv(\"Prem_work.csv\")\n",
+    "Mwork = Mworkdf.values.tolist()\n",
+    "Mrestdf = pd.read_csv(\"Prem_rest.csv\")\n",
+    "Mrest = Mrestdf.values.tolist()\n",
+    "\n",
+    "# Age brackets\n",
+    "Brac = ('Y','M','O')\n",
+    "Bracind = {'Y': (0,1,2,3,4), 'M': (5,6,7,8,9,10,11,12), 'O':(13,14,15)}\n",
+    "\n",
+    "# Determine contact matrices\n",
+    "allCont = {}\n",
+    "homeCont = {}\n",
+    "wsCont = {}\n",
+    "restCont = {}\n",
+    "for x, y in list(product(Brac,Brac)):\n",
+    "    tar = Bracind[y]\n",
+    "    src = Bracind[x]\n",
+    "    allSum = np.sum(Flatten([[N[j]*Mall[i][j] for i in tar] for j in src]))\n",
+    "    homeSum = np.sum(Flatten([[N[j]*Mhome[i][j] for i in tar] for j in src]))\n",
+    "    sSum = np.sum(Flatten([[N[j]*Mschool[i][j] for i in tar] for j in src]))\n",
+    "    wSum = np.sum(Flatten([[N[j]*Mwork[i][j] for i in tar] for j in src]))\n",
+    "    rSum = np.sum(Flatten([[N[j]*Mrest[i][j] for i in tar] for j in src]))\n",
+    "    Num = np.sum([N[i] for i in src])\n",
+    "    allCont[(x,y)] = allSum/Num\n",
+    "    homeCont[(x,y)] = homeSum/Num\n",
+    "    wsCont[(x,y)] = (wSum+sSum)/Num\n",
+    "    restCont[(x,y)] = rSum/Num"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Saving a copy of the mixing matrices\n",
+    "swap = ['Y', 'M', 'O']\n",
+    "\n",
+    "wsMat = np.full([3,3], 0.)\n",
+    "restMat = np.full([3,3], 0.)\n",
+    "for i, j in list(product(range(3), range(3))):\n",
+    "    wsMat[i][j] = wsCont[(swap[i], swap[j])]\n",
+    "    restMat[i][j] = restCont[(swap[i], swap[j])]\n",
+    "\n",
+    "pd.DataFrame(wsMat).to_csv(\"wsMixing.csv\", index=False)\n",
+    "pd.DataFrame(restMat).to_csv(\"restMixing.csv\", index=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[2.22220876 1.27618792 0.15504874]\n",
+      " [0.84971392 3.72295479 0.74083542]\n",
+      " [0.21683524 0.782842   0.66709872]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(restMat)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Rescaling workschool distributions from POLYMOD to Prem\n",
+    "\n",
+    "Since the workschool distributions from POLYMOD are fit to a 2-variable distribution we have to make a choice as to how to adjust those parameters. We first modify only p, keeping alpha from POLYMOD fixed."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Input the POLYMOD parameters\n",
+    "POLY = {('Y', 'Y'): [0.4378787156156721, 0.093112367309826], ('Y', 'M'): [0.4144210645686216, 0.2988493341114929], ('Y', 'O'): [0.9071801708997953, 0.8256463657373431], ('M', 'Y'): [0.76732665579341, 0.41260047034256747], ('M', 'M'): [0.421073517853328, 0.12039762005212719], ('M', 'O'): [0.8424225575069703, 0.5573421127782885], ('O', 'Y'): [0.8976481193486225, 0.40369064730244425], ('O', 'M'): [0.7855029449186991, 0.33072950060358974], ('O', 'O'): [0.9407010943356611, 0.4756835100017134]}\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Premparams = {}\n",
+    "for x, y in list(product(Brac,Brac)):\n",
+    "    paramval = (1-POLY[(x,y)][0])/wsCont[(x,y)]\n",
+    "    if paramval <1:\n",
+    "        Premparams[(x,y)] = ['geo', POLY[(x,y)][0], paramval ]\n",
+    "    else:\n",
+    "        Premparams[(x,y)] = ['pois', POLY[(x,y)][0], 1/paramval ]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(Premparams)"
+   ]
+  },
+  {
+   "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.8.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
+%% Cell type:code id: tags:
+``` python
+%matplotlib notebook
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+from scipy import optimize
+import math
+import random
+import csv
+import itertools
+from itertools import product
+def Flatten(A):
+    return list(itertools.chain.from_iterable(A))
+```
+%% Cell type:markdown id: tags:
+# Average mixing numbers for Canadian population
+To adjust the POLYMOD results to fit Canadian demographic structures, we take the distributions extracted from POLYMOD and rescale them so that their averages match those determined by Prem.
+%% Cell type:code id: tags:
+``` python
+# Population statistics from Statcan 2019
+N = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
+N[0] = 1943175
+N[1] = 2039352
+N[2] = 2031762
+N[3] = 2114635
+N[4] = 2476698
+N[5] = 2625474
+N[6] = 2603938
+N[7] = 2580021
+N[8] = 2421009
+N[9] = 2396406
+N[10] = 2502667
+N[11] = 2749626
+N[12] = 2511888
+N[13] = 2096607
+N[14] = 1706760
+N[15] = 1164277+786704+510828+242554+74086+10795 #75+
+Ntot = np.sum(N)
+# Load the Prem files
+Malldf = pd.read_csv("Prem_Canada.csv")
+Mall = Malldf.values.tolist()
+Mhomedf = pd.read_csv("Prem_home.csv")
+Mhome = Mhomedf.values.tolist()
+Mschooldf = pd.read_csv("Prem_school.csv")
+Mschool = Mschooldf.values.tolist()
+Mworkdf = pd.read_csv("Prem_work.csv")
+Mwork = Mworkdf.values.tolist()
+Mrestdf = pd.read_csv("Prem_rest.csv")
+Mrest = Mrestdf.values.tolist()
+# Age brackets
+Brac = ('Y','M','O')
+Bracind = {'Y': (0,1,2,3,4), 'M': (5,6,7,8,9,10,11,12), 'O':(13,14,15)}
+# Determine contact matrices
+allCont = {}
+homeCont = {}
+wsCont = {}
+restCont = {}
+for x, y in list(product(Brac,Brac)):
+    tar = Bracind[y]
+    src = Bracind[x]
+    allSum = np.sum(Flatten([[N[j]*Mall[i][j] for i in tar] for j in src]))
+    homeSum = np.sum(Flatten([[N[j]*Mhome[i][j] for i in tar] for j in src]))
+    sSum = np.sum(Flatten([[N[j]*Mschool[i][j] for i in tar] for j in src]))
+    wSum = np.sum(Flatten([[N[j]*Mwork[i][j] for i in tar] for j in src]))
+    rSum = np.sum(Flatten([[N[j]*Mrest[i][j] for i in tar] for j in src]))
+    Num = np.sum([N[i] for i in src])
+    allCont[(x,y)] = allSum/Num
+    homeCont[(x,y)] = homeSum/Num
+    wsCont[(x,y)] = (wSum+sSum)/Num
+    restCont[(x,y)] = rSum/Num
+```
+%% Cell type:code id: tags:
+``` python
+# Saving a copy of the mixing matrices
+swap = ['Y', 'M', 'O']
+wsMat = np.full([3,3], 0.)
+restMat = np.full([3,3], 0.)
+for i, j in list(product(range(3), range(3))):
+    wsMat[i][j] = wsCont[(swap[i], swap[j])]
+    restMat[i][j] = restCont[(swap[i], swap[j])]
+pd.DataFrame(wsMat).to_csv("wsMixing.csv", index=False)
+pd.DataFrame(restMat).to_csv("restMixing.csv", index=False)
+```
+%% Cell type:code id: tags:
+``` python
+print(restMat)
+```
+%% Output
+    [[2.22220876 1.27618792 0.15504874]
+     [0.84971392 3.72295479 0.74083542]
+     [0.21683524 0.782842   0.66709872]]
+%% Cell type:markdown id: tags:
+## Rescaling workschool distributions from POLYMOD to Prem
+Since the workschool distributions from POLYMOD are fit to a 2-variable distribution we have to make a choice as to how to adjust those parameters. We first modify only p, keeping alpha from POLYMOD fixed.
+%% Cell type:code id: tags:
+``` python
+# Input the POLYMOD parameters
+POLY = {('Y', 'Y'): [0.4378787156156721, 0.093112367309826], ('Y', 'M'): [0.4144210645686216, 0.2988493341114929], ('Y', 'O'): [0.9071801708997953, 0.8256463657373431], ('M', 'Y'): [0.76732665579341, 0.41260047034256747], ('M', 'M'): [0.421073517853328, 0.12039762005212719], ('M', 'O'): [0.8424225575069703, 0.5573421127782885], ('O', 'Y'): [0.8976481193486225, 0.40369064730244425], ('O', 'M'): [0.7855029449186991, 0.33072950060358974], ('O', 'O'): [0.9407010943356611, 0.4756835100017134]}
+```
+%% Cell type:code id: tags:
+``` python
+Premparams = {}
+for x, y in list(product(Brac,Brac)):
+    paramval = (1-POLY[(x,y)][0])/wsCont[(x,y)]
+    if paramval <1:
+        Premparams[(x,y)] = ['geo', POLY[(x,y)][0], paramval ]
+    else:
+        Premparams[(x,y)] = ['pois', POLY[(x,y)][0], 1/paramval ]
+```
+%% Cell type:code id: tags:
+``` python
+print(Premparams)
+```
+%% Cell type:code id: tags:
+``` python
+```