source/analyses/fitting/Massively_Parallel_Fitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parallel Model Fitting\n",
    "\n",
    "In this notebook, we'll fit a simple ageing model to all of the counties in the United States. As before, we'll use `scipy.optimize` to perform the fitting, but we'll use python's `multiprocessing` library to perform these optimizations in parallel.\n",
    "\n",
    "## When to use this technique\n",
    "This technique is appropriate when we are modeling a large number of entirely independent but structurally identical systems. In this example, we're conceptualizing the population of counties to be influenced by aging and exogenous migration patterns. If we were to attempt to link the models together, for instance by specifying that the outmigration from one county needed to be accounted for as the inmigration to another county, we would need to perform a single large-scale optimization, or some form of hybrid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import pandas as pd\n",
    "import pysd\n",
    "import scipy.optimize\n",
    "import multiprocessing\n",
    "import numpy as np\n",
    "import seaborn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ingredients\n",
    "\n",
    "#### Data\n",
    "The first ingredient theat we'll use is census data from the 2000 and 2010 census:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('../../data/Census/Males by decade and county.csv', header=[0,1], skiprows=[2])\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Model\n",
    "The model will be a simple ageing chain that groups individuals into 10 year cohorts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pysd.read_vensim('../../models/Aging_Chain/Aging_Chain.mdl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Recipe\n",
    "As in our other optimization problems, we'll construct an error function that calculates the sum of squared errors between our model prediction and the measured data. We also construct a helper function called `fit` which basically makes the call to the optimizer and formats the result in something that we can aggregate into a Pandas DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "param_names = ['dec_%i_loss_rate'%i for i in range(1,10)]\n",
    "\n",
    "def error(param_vals, measurements):\n",
    "    predictions = model.run(params=dict(zip(param_names, param_vals)),\n",
    "                            initial_condition=(2000, measurements['2000']),\n",
    "                            return_timestamps=2010).loc[2010]\n",
    "\n",
    "    errors = predictions - measurements['2010']\n",
    "    return sum(errors.values[1:]**2) #ignore first decade: no birth info\n",
    "\n",
    "def fit(row):\n",
    "    res = scipy.optimize.minimize(error, args=row,\n",
    "                                  x0=[.05]*9,\n",
    "                                  method='L-BFGS-B');\n",
    "    return pd.Series(index=['dec_%i_loss_rate'%i for i in range(1,10)], data=res['x'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, fitting the model is a simple matter of applying the fit function to the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "county_params = data.apply(fit, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On my 2014 era machine, this optimization takes about half an hour.\n",
    "\n",
    "We can plot the distributions of the fit parameters for each of the counties in a histogram, to get a sense of the result. (Here we're ignoring the first decade, which will not have reasonable parameters, as we have no information about births to the system.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df2 = county_params.drop('dec_1_loss_rate',1)\n",
    "df2.plot(kind='hist', bins=np.arange(-.15,.4,.01), alpha=.4, histtype='stepfilled')\n",
    "plt.xlim(-.15,.4)\n",
    "plt.title('Fit yearly loss rates from each US county\\n by age bracket from 2000 to 2010', fontsize=16)\n",
    "plt.ylabel('Number of Counties', fontsize=16)\n",
    "plt.xlabel('Yearly Loss Rate in 1% Brackets', fontsize=16)\n",
    "plt.legend(frameon=False, fontsize=14)\n",
    "plt.savefig('Loss_Histogram.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Executing the optimization in parallel\n",
    "\n",
    "We can take advantage of the multicore nature of most modern machines by using python's `multiprocessing` module to distribute the various counties between each of the cores we have available for the calculation. The basic structure for this piece of code comes from [this gist](https://gist.github.com/yong27/7869662). We are essentially creating a helper function that will apply the fit function to a subset of the census DataFrame, and calling this function once on each of our worker nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "\n",
    "def _apply_df(args):\n",
    "    df, func, kwargs = args\n",
    "    return df.apply(func, **kwargs)\n",
    "\n",
    "def apply_by_multiprocessing(df, func, workers, **kwargs):\n",
    "    pool = multiprocessing.Pool(processes=workers)\n",
    "    result = pool.map(_apply_df, [(d, func, kwargs) for d in np.array_split(df, workers)])\n",
    "    pool.close()\n",
    "    return pd.concat(list(result))\n",
    "\n",
    "county_params = apply_by_multiprocessing(data[:10], fit, axis=1, workers=4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}