# Parallel Model Fitting¶

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.

## When to use this technique¶

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.

```
%pylab inline
import pandas as pd
import pysd
import scipy.optimize
import multiprocessing
import numpy as np
import seaborn
```

```
Populating the interactive namespace from numpy and matplotlib
```

## Ingredients¶

### Data¶

The first ingredient theat we’ll use is census data from the 2000 and 2010 census:

```
data = pd.read_csv('../../data/Census/Males by decade and county.csv', header=[0,1], skiprows=[2])
data.head()
```

Unnamed: 0_level_0 | Unnamed: 1_level_0 | 2000 | 2010 | |||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Unnamed: 0_level_1 | Unnamed: 1_level_1 | dec_1 | dec_2 | dec_3 | dec_4 | dec_5 | dec_6 | dec_7 | dec_8 | dec_9 | dec_1 | dec_2 | dec_3 | dec_4 | dec_5 | dec_6 | dec_7 | dec_8 | dec_9 | |

0 | 1 | 1 | 3375 | 3630 | 2461 | 3407 | 3283 | 2319 | 1637 | 825 | 284 | 3867 | 4384 | 3082 | 3598 | 4148 | 3390 | 2293 | 1353 | 454 |

1 | 1 | 3 | 9323 | 10094 | 7600 | 9725 | 10379 | 8519 | 6675 | 4711 | 1822 | 11446 | 12006 | 9976 | 11042 | 12517 | 12368 | 10623 | 6307 | 2911 |

2 | 1 | 5 | 2002 | 2198 | 2412 | 2465 | 2178 | 1699 | 1026 | 689 | 301 | 1673 | 1739 | 2260 | 2208 | 2233 | 1910 | 1490 | 739 | 324 |

3 | 1 | 7 | 1546 | 1460 | 1680 | 1762 | 1624 | 1237 | 774 | 475 | 187 | 1471 | 1577 | 1798 | 2016 | 1928 | 1581 | 1140 | 579 | 211 |

4 | 1 | 9 | 3741 | 3615 | 3393 | 3901 | 3773 | 3007 | 2227 | 1269 | 550 | 3741 | 4252 | 3312 | 3719 | 4129 | 3782 | 3052 | 1723 | 652 |

### Model¶

The model will be a simple ageing chain that groups individuals into 10 year cohorts.

```
model = pysd.read_vensim('../../models/Aging_Chain/Aging_Chain.mdl')
```

## The Recipe¶

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.

```
param_names = ['dec_%i_loss_rate'%i for i in range(1,10)]
def error(param_vals, measurements):
predictions = model.run(params=dict(zip(param_names, param_vals)),
initial_condition=(2000,measurements['2000']),
return_timestamps=2010,
rtol=1).loc[2010]
errors = predictions - measurements['2010']
return sum(errors.values[1:]**2) #ignore first decade: no birth info
def fit(row):
res = scipy.optimize.minimize(error, args=row,
x0=[.05]*9,
method='L-BFGS-B');
return pd.Series(index=['dec_%i_loss_rate'%i for i in range(1,10)], data=res['x'])
```

At this point, fitting the model is a simple matter of applying the fit function to the data:

```
%%capture
county_params = data.apply(fit, axis=1)
```

On my 2014 era machine, this optimization takes about half an hour.

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.)

```
df2 = county_params.drop('dec_1_loss_rate',1)
df2.plot(kind='hist', bins=np.arange(-.15,.4,.01), alpha=.4, histtype='stepfilled')
plt.xlim(-.15,.4)
plt.title('Fit yearly loss rates from each US county\n by age bracket from 2000 to 2010', fontsize=16)
plt.ylabel('Number of Counties', fontsize=16)
plt.xlabel('Yearly Loss Rate in 1% Brackets', fontsize=16)
plt.legend(frameon=False, fontsize=14)
plt.savefig('Loss_Histogram.svg')
```

### Executing the optimization in parallel¶

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. 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.

```
%%capture
def _apply_df(args):
df, func, kwargs = args
return df.apply(func, **kwargs)
def apply_by_multiprocessing(df, func, workers, **kwargs):
pool = multiprocessing.Pool(processes=workers)
result = pool.map(_apply_df, [(d, func, kwargs) for d in np.array_split(df, workers)])
pool.close()
return pd.concat(list(result))
county_params = apply_by_multiprocessing(data[:10], fit, axis=1, workers=4)
```