gistfile1.txt

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A primer for working with the <u>*Sp*</u>atial <u>*Int*</u>eraction modeling (SpInt) module in the python spatial analysis library (PySAL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spatial interaction modeling involves the analysis of flows from an origin to a destination either over physical space (i.e., migration) or through abstract space (i.e., telecommunication). Despite being a fundamental technique to regional science (citation), there is relatively little software available to carry out spatial interaction modeling and the analysis of flow data. Therefore, the purpose of this primer is to provide an overview of the recently develped spatial interaction modeling (SpInt) module of the python spatial analysis library (PySAL). First, the current framework of the module will be highlighted. Next, the main functionality of the module will be illustrated with a common regional science example, migration flows, with a dataset previously used for spatial interaction modeling tutorials in the R programming environment. Finally, some auxilliary functionality and future additions are discussed. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 2 The SpInt Framework"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Modeling framework"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The core purpose of the SpInt module is to provide the functionality to calibrate spatial interaction models. Since the \"family\" of spatial interaction models put forth by Wilson (Wilson, 1971) are perhaps the most popular, they were chosen as the starting point of the module. Consider the basic gravity model (Fotheringham & O'Kelly, 1989), \n",
    "\n",
    "$$T_{ij} = k\\frac{V_{i}^\\mu W_{j}^\\alpha}{d_{ij}^\\beta} \\quad(1)$$\n",
    "where \n",
    "\n",
    "$T_{ij}$ = an $n \\times m$ matrix of flows between $n$ origins (subscripted by $i$) to $m$ destinations (subscripted by $j$)\n",
    "\n",
    "$V$ = an $n \\times p$ and vector of $p$ origin attributes\n",
    "\n",
    "$W$ = an $m \\times p$ vector of $p$ destination attributes\n",
    "\n",
    "$d$ = an $n \\times m$ matrix of the costs to overcome the physical separation between  $i$ and $j$ (usually distance or time)\n",
    "\n",
    "$k$ = a scaling factor to be estimated\n",
    "\n",
    "$\\mu$ = a $p \\times 1$ vector of parameters representing the effect of $p$ origin attributes on flows\n",
    "\n",
    "$\\alpha$ = a $p \\times 1$ vector of parameters representing the effect of $p$ destination attributes on flows\n",
    "\n",
    "$\\beta$ = an exponential parameter representing the effect of movement costs on flows. \n",
    "\n",
    "When data for $T$, $V$,  $W$, and $d$ are available we can estimate the model parameters (also called calibration), which summarize the effect that each model component contributes towards explaining the system of known flows ($T$). In contrast, known parameters can be used to predict unknown flows when there are deviations in model components ($V$, $W$, and $d$) or the set of locations in the system are altered.\n",
    "\n",
    "Using an entropy-maximizing framework, Wilson derives a more informative and flexible \"family\" of four spatial interaction models (Wilson, 1971). This framework seeks to assign flows between a set of origins and destinations by finding the most probable configuration of flows out of all possible configurations, without making any additional assumptions. By using a common optimization problem and including information about the total inflows and outflows at each location (also called constraints), the following \"family\" of models can be obtained:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "&Unconstrained \\ (Gravity) \\\\\n",
    "&Tij = V_{i}^\\mu W_{j}^\\alpha  f(d_{ij}) \\quad & (2) \\\\\n",
    "\\\\\n",
    "&Production-constrained \\\\\n",
    "&T_{ij} = A_{i}O_{i}W_{j}^\\alpha f(d_{ij}) \\quad & (3) \\\\\n",
    "&A_{i} = \\sum_{j} W_{j}^\\alpha f(d_{ij}) \\quad & (3a) \\\\\n",
    "\\\\\n",
    "&Attraction-constrained \\\\\n",
    "&T_{ij} = B_{j}D_{j}V_{i}^\\mu f(d_{ij}) \\quad & (4) \\\\\n",
    "&B_{j} = \\sum_{i} V_{i}^\\mu f(d_{ij}) \\quad & (4a) \\\\\n",
    "\\\\\n",
    "&Doubly-constrained \\\\\n",
    "&T_{ij} = A_{i}B_{j}O_{i}D_{j}f(d_{ij}) \\quad & (5) \\\\\n",
    "&A_{i} = \\sum_{j} W_{j}^\\alpha B_{j} D_{j} f(d_{ij}) \\quad & (5a) \\\\\n",
    "&B_{j} = \\sum_{i} V_{i}^\\mu A_{i} O_{i} f(d_{ij}) \\quad & (5b)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "where \n",
    "\n",
    "$O_{i}$ = an $n \\times 1$ vector of the total number of flows emanating from origin $i$\n",
    "\n",
    "$D_{j}$ = an $m \\times 1$ vector of the total number of flows terminating at destination $j$\n",
    "\n",
    "$A_{i}$ = an $n \\times 1$ vector of thevorigin balancing factors that ensures the total out-flows are preserved in the predicted flows\n",
    "\n",
    "$B_{j}$ = an $m \\times 1$ vector of the destination balancing factors that ensures the total in-flows are preserved in the predicted flows\n",
    "\n",
    "$f(d_{ij})$ = a function of cost or distance, referred to as the distance-decay function. Most commonly this an exponential or power function given by,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "&Power\\\\\n",
    "&f(d_{ij}) = d_{ij}^\\beta \\quad & (6) \\\\\n",
    "\\\\\n",
    "&Exponential \\\\\n",
    "&f(d_{ij}) = exp(\\beta*d_{ij}) \\quad & (7) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $\\beta$ is expected to take a negative value. Different distance-decay functions assume different responses to higher costs associated with moving to more distant locations. Of note is that the unconstrained model with of power function distance-decay is equivalent to the basic gravity model in equation (2), except that the scaling factor, $k$, is not included. In fact, there is no scaling factor in any of the members of the family of maximum entropy models because there is a total trip constrained implied in their derivation and subsequently incorporated into their calibration (Fotheringham & O'Kelly , 1989). Another aside is this in the doubly-constrained maximum entoropy model the values for $A_{i}$ and $B_{j}$ are dependent upon each other and may need to be computed iteratively depending on calibration technqiue. It is also usually assumed that $n=m$ for doubly-constrained models.\n",
    "\n",
    "The family provides different model strucutre depedning on the available data or the research question at hand. The so-called unconstrained model does not conserve the total inflows or outflows during parameter estimation. The production-constrained and attraction-constrained models conserve either the number of total inflows or outflows, respectively, and are therefore useful for building models that allocate flows either to a set of origins or to a set of destinations. Finally, the doubly-constrained model conserves both the inflows and the outflows at each location during model calibration. The quantity of explanatory information provided by each model is given by the number of parameters it provides. As such, the unconstrained model provides the most information, followed by the two singly-constrained models, with the doubly-constrained model providing the least information. Conversely, the model's predictive power increases with higher quantities of built-in information (i.e. total in or out-flows) so that the doubly-constrained model usually provides the most accurate predictions, followed by the two singly-constrained models, and the unconstrained model supplying the weakest predictions (Fotheringham & O'Kelly, 1989). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Calibration framework"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spatial interaction models are often calibrated via linear programming, nonlinear optimization, or increasingly more often through linear regression. Given the flexibility and extendability of a regression framework it was chosen as the primary model calibration technqiue within the SpInt module. By taking the natural logarithm of both sides of a spatial interaction model, say the basic gravity model, is is possible to obtain the so-called log-linear or log-normal spatial interaction model:\n",
    "\n",
    "\n",
    "$$\\ln{T_{ij}} = k + \\mu \\ln{V_{i}} + \\alpha \\ln{W_{j}} - \\beta \\ln{d_{ij}} + \\epsilon \\quad (8)$$\n",
    "\n",
    "\n",
    "where $\\epsilon$ is a a normally distributed error term with a mean of 0. Constrained spatial interaction models can be achieved by including a fixed effects for the origins (production-constrained), a fixed effect for the destinations (attraction-constrained) or both (doubly-constrained). However, there are several limitations of the log-normal gravity model, which include:\n",
    "\n",
    "1. flows are often counts of people or objects and should be modeled as discrete entities; \n",
    "2. flows are often not normally distributed;\n",
    "3. Bias\n",
    "4. zero flows \n",
    "\n",
    "Therefore, Flowerdew (1984) proposes the Poisson log-linear regression specification for the family of spatial interaction models. This specification assumes that the number of flows between $i$ and $j$ is drawn from a Poisson distribution with mean, $\\lambda_{ij} = T_{ij}$, where $\\lambda_{ij}$ is assumed to be logarithmically linked to the linear combination of variables,\n",
    "\n",
    "$$\\ln{\\lambda_ij} = k + \\mu V_{i} + \\alpha W_{j} - \\beta d_{ij}) \\quad (9)$$\n",
    "\n",
    "and exponentiating both sides of the equation yields the following family of Poisson log-linear spatial interaction models:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "&Unconstrained \\\\\n",
    "&T_{ij} = \\exp(k + \\mu V_{i} + \\alpha W_{j} - \\beta d_{ij}) \\quad & (10) \\\\\n",
    "\\\\\n",
    "&Production-constrained \\\\\n",
    "&T_{ij} = \\exp(k + \\mu_{i} + \\alpha W_{j} - \\beta d_{ij}) \\quad & (11) \\\\\n",
    "\\\\\n",
    "&Attraction-constrained \\\\\n",
    "&T_{ij} = \\exp(k + \\mu V_{i} + \\alpha_{j} - \\beta d_{ij}) \\quad & (12) \\\\\n",
    "\\\\\n",
    "&Doubly-constrained \\\\\n",
    "&T_{ij} = \\exp(k + \\mu_{i} + \\alpha_{j} - \\beta d_{ij}) \\quad & (13) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $\\mu_{i}$ are origin fixed effects and $\\alpha_{i}$ are destination fixed effects that achieve the same constraints as the balancing factors in equations (2-5). Using Poisson regression automatically alleviates limiations (1) and (2) and since we no longer need to take the logarithm of $T_{ij}$, issue (4) is also absolved. Using fixed effects within Poisson regression to calibrate the doubly-constrained model also avoids the need for iterative computation of the balancing factors that exists in other calibration methods. \n",
    "\n",
    "Calibration of Poisson regression can be carried out within a generalized linear modeling framework (GLM) using iteratively weighted least sqaures (IWSL), which converges to the maximum likelihood estimates for the parameter estimates (Nelder & Wedderburn, 1972). To maintain computational efficiency with increasingly larger spatial interaction datasets, SpInt is built upon a custom GLM/IWLS routine that leverages sparse data structures for the production-constrained, attraction-constrained, and doubly-constrained models. As the number of locations in these models increases, so the does the number of binary indicator variables needed to construct the fixed effects that enforce the constraints. Therefore, larger spatial interaction datasets become increasingly sparse and the utilization of sparse data structures take advantage of this feature. As a metric, constrained models with $n = 3,000$ locations, which implies $n^2 = 9,000,000$ observations can still be calibrated within minutes on a standard macbook pro notebook.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 An Illustrative example: migration in Austria"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 The data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following example purposefully utilizes data that has was previously used to demonstrate spatial interaction modeling the R programming language for validation and consistency. The data is migration flows between Austrian NUTS level 2 municipalities in 2006. In order to use a regression-based calibration, the data has to be transformed from the matrices and vectors described in equations (1-5) to a table where each row represents a single origin-destination dyad, $(i,j)$ and any variables associated with either location. Details on how to do this are outlined further in (LeSage & Pace, 2008), though this has already been done in the example data. Let's have a look! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Origin</th>\n",
       "      <th>Destination</th>\n",
       "      <th>Data</th>\n",
       "      <th>Oi</th>\n",
       "      <th>Dj</th>\n",
       "      <th>Dij</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AT11</td>\n",
       "      <td>AT11</td>\n",
       "      <td>0</td>\n",
       "      <td>4016</td>\n",
       "      <td>5146</td>\n",
       "      <td>1.000000e-300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AT11</td>\n",
       "      <td>AT12</td>\n",
       "      <td>1131</td>\n",
       "      <td>4016</td>\n",
       "      <td>25741</td>\n",
       "      <td>1.030018e+02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AT11</td>\n",
       "      <td>AT13</td>\n",
       "      <td>1887</td>\n",
       "      <td>4016</td>\n",
       "      <td>26980</td>\n",
       "      <td>8.420467e+01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>AT11</td>\n",
       "      <td>AT21</td>\n",
       "      <td>69</td>\n",
       "      <td>4016</td>\n",
       "      <td>4117</td>\n",
       "      <td>2.208119e+02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>AT11</td>\n",
       "      <td>AT22</td>\n",
       "      <td>738</td>\n",
       "      <td>4016</td>\n",
       "      <td>8634</td>\n",
       "      <td>1.320075e+02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Origin Destination  Data    Oi     Dj            Dij\n",
       "0   AT11        AT11     0  4016   5146  1.000000e-300\n",
       "1   AT11        AT12  1131  4016  25741   1.030018e+02\n",
       "2   AT11        AT13  1887  4016  26980   8.420467e+01\n",
       "3   AT11        AT21    69  4016   4117   2.208119e+02\n",
       "4   AT11        AT22   738  4016   8634   1.320075e+02"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAACjCAYAAABmOs66AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnWV0VFcXhp8bdyNKEiBo0ODubsVdW+CjOLRQoIKWosWt\nWHGKFXeXUCyQAEECASIQd89kcs/3IymFEkoSYsB91mIxM/fIvqxhz7n77LNfSQiBgoKCgsLngUZ+\nG6CgoKCgkHcoTl9BQUHhM0Jx+goKCgqfEYrTV1BQUPiMUJy+goKCwmeE4vQVFBQUPiO0MttQkiQN\n4BbgL4ToIEmSC7AaMAR8gL5CiLgM+vkA0YAMpAghauaA3QoKCgoK2SArK/2xwP3X3q8DJgohXID9\nwMR39JOBxkKIKorDV1BQUMhfMuX0JUlyANoC61/7uLQQwjX99Rmg67u6Z3YeBQUFBYXcJbPOeDHw\nHfD68V1PSZI6pL/uATi8o68ATkuSdFOSpP9lz0wFBQUFhZzgvU5fkqR2QLAQwoO0VfvfDAZGSpJ0\nk7S4vuodQ9QTQlQl7UlhpCRJ9T/QZgUFBQWFbJKZjdx6QAdJktoC+oCxJElbhBADgFYAkiSVAtpl\n1FkIEZj+d6gkSfuBmoDrv9tJkqQUAVJQUFDIIkII6f2t/uG9K30hxA9CiCJCiOJAL+CcEGKAJElW\n8Cqr5yfgt3/3lSTJQJIko/TXhkBLwPM/5ipQf6ZNm5bvNig2fTo2FVS7FJs+Xpuyw4dssPaWJMkL\neAC8FEJsApAkyU6SpCPpbWwAV0mS3IFrwGEhxKkPmFNBQUFB4QPIdJ4+gBDiInAx/fUyYFkGbQKB\n9umvnwOVP9xMBQUFBYWcQEml/A8aN26c3ya8hWJT5iiINkHBtEuxKXMURJuyg5TduFBOI0mSKCi2\nKCgoKHwMSJKEyOmNXAUFBQWFTwfF6SsoKCh8RihOX0FBQeEzQnH6CgoKCp8RitNXUFBQ+IxQnL6C\ngoLCZ4Ti9BUUFBQ+IxSnr6CgoPAZkWmnL0mShiRJ7pIkHUp/7yJJ0l+SJN2RJOng34XVMujXWpKk\nR5IkPZYkaVJOGa6goKCgkHVyVS4xvQLnCtJKMJcnrUibc/bNVVBQUFD4EHJbLrEm8EQI4SuESAF2\nAh0/wF4FBYUcRpbl/DZBIQ/JbblEe8D/tfcv0j9TUFDIZ1JTU3G0t0dTU5PiRRzp1rkzV69e5fTp\n0wQFBb3RVghBUlJSPlmqkJO8t7Ty63KJkiQ1fu3SYGCZJElTgEO8Wy4x00yfPv3V68aNG38yVe0U\nFAoi169f50VAAGeG98c/KoaV190Z1KMbFgYG/PXYGz8/PxwdHQHYunUrAwcOZED//kQGB1GsREnG\njh9PiRIl8vkuPi8uXLjAhQsXPmiM91bZlCRpNtAPUJMulwjsE2lyiX+3KQVsFULU/lff2sB0IUTr\n9PeTASGEmJfBPEqVTQWFPGTMyBFYPnvI1FaN3vhclgWa42dSuXJlyhR34rGXF+73H+BsY0VXl7KU\ns7bkUVgEq67epl27dpQo44yjgwM9e/XCwMAgn+7m8yQ7VTazVFpZkqRGwHghRAdJkqxEmu6tBrAR\nOP+3etZr7TUBL6AZEAjcAHoLIR5mMLbi9BUU8pAfvv+ePb+v58n3I9+61mvHQR4Gh9GzkjMtSzlx\nPyiEvtUqoaX5T0Q4Ij6RLW538AyNwC8qhqexidz38kJPTy8vb+OzJq+d/hhgJGlx/n1CiB/S29gB\n64QQ7dPftwaWkrZ/sEEIMfcdYytOX0EhD7lx4wa1atVCLJ6WI+N1336AY55e1K9di+Zt2qKhqUnn\nzp0pXrx4joyv8Da57vRzE8XpKyjkLV8NHMi5o4fx/WlMjo0Zm5TMiUfeXPELICIhkZ237tGoXl1+\n37oNY2NjjI2NmTd7Nn4+PkiamlhbWfHjtGno6OjkmA2fE4rTV1BQyDSmJiZ0KVuCjb1zL4s6VZb5\n+Ywrv5y6hK6ONro6urgUtqFH+VIcvv+YYw+fANCwQQO2bttGkSJFcs2WTxHF6SsoKGSKqKgozM3N\n2dG/C72rVsz1+VJlGQkJv6ho7E2N0dbUBCAsLoFaS9bzLDwSAKdixXj2/PkHz5eUlERQUBDFihXD\nx8eHqKgoKleu/MHjFjSy4/Tfm7KpoKDwcRMZGcmDBw8oXLgwpqamWFhYsHfvXio62NG1Urk8sUFT\nI20DuJiF2RufWxoZ8PSnMajUqbRbt4Mzj5+xa+dOuvfogYZG9kqDybJM+bJleebjg76uLonJyVia\nmtJ/4EBatmmDjo4OlpaWVKpU6YPv62NEWekrKOQgkZGR9P+yP6oUFV07d2XIoCFopq9q8wtNTU1k\nWUZXW5vklBRaNGyASSFL/ty/n8WdWjGuUe33D5JHXPT2YcLxiwTGJ9C5cxe69OhBgwYN0NLK/Pp0\n3bp1zJv6Ex7jBqOWZRJUKcQlq1h88RoX/AIw1tPF0/8lU6dMZfKPP+bi3eQ+SnhHQSGf6dOvD96h\n3rT9ui0HFh+ABDi47yBFixbNF3uCg4OxtbUlcMZ4bE2MeBIazuOQcNxeBqGhITG2fk1M9HTzxbb/\n4mFwKPvvebHv0VN8wiMp5uBAXEICcfHxOBYujIODI1a2tgS+eIFzhQocPXaMGjVqkJKSwrHDhzgw\nsBv1i797f2Dd1VsM3X2EyRMnUrJUKfr07Yu+vv6r60IIJClLvjRfUJy+gkI+cuToEb5o/wVzT86l\nesvqCCHYt2QfGyZvoHPXzlSpXIWh/xuKubl5rtuSlJTEpo0bWbNyBSYpyZwc0gs97Y8zmvsyKoaX\n0bGoUlNxMDN59T4wJhY7E2P+8g/g1CNv4lQpVHWwY3brxpSztfrPMR+HhDPuyFlqFbZm1z0vEjQ0\nqVq1GqbGxqSmqvnz4CHq1qzJ9t27sba2zqM7zTqK01dQyEeePXv2qizBGXHm1edBPkFcP3qdNRPW\noEpS4ebmRtWqVXN1JdmtYwe877jzbb3q9KtWCQ2Ngr9qzS9SUlM59uAJbi8CMdLRRl9bG1tjI6ae\nvEjtlq3ZtGVLfpv4ThSnr6CQzzRs0pAqParQYXiHt66lpqYypvYYvNy8uHTpEg0aNMjS2EIIwsPD\nMTIy+s9Tr/fv36dK5cpE/zIRfR3tLN+DQhpf/nGQzTc8KMh+KTtOX1HOUlDIQQoVKkREYESG185s\nPYOXmxdOTk7UqFEjy2MfP34cKysripcszuHDh0lJSQHAzc2NAwcOvGrn5eWFvo624vA/kD5VKgCf\nXunpD1HOqixJ0tX0z25IklT9Hf180tW13CVJupFThiso5BWJiYmZahcfH8+zp8/Q0894Fd6wW0PK\n1ynP8+fPadGiBUeOHOHkyZPcvn2bp0+fkpCQkGE/lUqFj48PPj4+WDtY029WP/43/H8YGhpiZGRE\ni9Yt6DegH81bNcfd3Z24uDhK2RbcOPTHQmFTYwAWL1qUz5bkLFnZ2flbOcsk/f08YJoQ4pQkSW2A\nBUCTDPrJQGMhROQHWaqgkIcIIZBlmXXr1jF8+HCatWzGiqUriIyMpEqVKm+FV+Lj49m+fTsP7j/g\n10m/ZjimvpE+iy4vYs34Nexbug9XV1dKu5QmPiYeBESERFC7bm109XTp16sf7dq14/79+0ycOJG/\n/voLs0JmfP/H91RrUY1WX7YiPiaeG8dvUK9TPQKfBXJg+QGqVq0KQCFDAyLiEzDT18t2vvvnTgU7\nawbUqU5CJn/0PxYyFdNPV87aCPwCfJtecO048LsQYo8kSb2BdkKIfhn0fQ5UF0KEv2cOJaavUGDo\n1LUTB/cdBGDS5knERsSy6adNaOtqY25uTscOHZkxbQampqacOnWKVq1aATB191Qadm+YrTmTEpK4\neugqiXGJXNlzheunrgNQt11d2o9qT83WNf+zvzpFzbN7zziw7ACnNp969bm+jjYdypXm144tcTAz\n+Y8RFP7NTycucF2lwelz5/LblAzJtY1cSZL2kObwTfmnyqYzcBKQ0v/UFUL4Z9D3GRAFpAJrhRDr\n3jGH4vQVCgzde3Tngf8DZh2dhYlFmqNUJavQ0tbi4bWHHFl9BK8rXvTp3Yf9+/bj9ciLHX47sHbM\nubDKs3vPkCQJpwpOWe772O0xIf4h6Bno4X7uNnsX7iU1VcbRwowB1Sowq22zHLPzUyU8PoGWv23D\nNz6JsIiM92nym1xx+unKWW2EEKPSlbP+XukvJa2G/gFJkroBXwshWmTQ304IEShJkhVwGhj1mrbu\n6+0Up69QYNi2bRtzls5h2c1lGV4XQuDl5sWVP6/wxO0JY9eNxc7JLo+tzDxxUXFc3H2RJcOXIGTB\n77068FWtKvltVoElLC4BqykLgDSFsZo1//spK7/Irdo79YAOkiS1JV05S5KkrUB7IcRYACHEXkmS\nNmTUWQgRmP53qCRJ+0kTS3/L6YMil6hQMEhNTWXBogV0ntT5nW0kScK5hjPONZzz0LLsY2RmxJHf\nDiNkwYDqlehcsWx+m1RgCYqJo+GqLWhpaXHjxg2qVCk4P455Ipf4RuM3RVTuAyOEEBclSWoGzBVC\n1PhXewNAQwgRJ0mSIXAKmCGEOJXB2MpKX6FA4OPjQ/kK5fkz4k+0P6G0x60zt7J52mbali/N0SG9\n89ucAkuxWcswsbPn2vXrBV7+Ma/z9IcCCyVJcgdmpb9HkiQ7SZKOpLexAVzT21wDDmfk8BUUChJG\nRkZoaGgQEVQw47jZpf/U/iy8uJATD73R+HYmRzy98tukAom2pibt2rYt8A4/u2SpGIcQ4iJwMf31\nFeCt3Pz0cE779NfPgU+viLXCJ83ylctp0rsJNkVs8tuUHMeloQuF7CwIfRnGiD+P0b5Cmfw2qUAR\nFpeAd0gYRZ2yvnn+saAk8CoopBMZGcnEyRNZunQpbb9um9/m5BqtBrWmcMnCvIyOzW9TChyuz/0A\naNeuXT5bknsoTl9BAbh16xYuVV14EPaA1XdWU6pqqfw2Kdf4cuaX6BvqYW1siFr9aZUYiEtSERaX\n8cnmzNC0lBPFbawoUqQI69ete1Xq4lNCKbim8FmRmppKbGwssiwjyzI3btxgxeoVXL16lTG/jaFh\nt+wdrPrY6FesL0G+wVibGBE8Y3x+m5MjLLpwlfEH07YMS9ta8fC7YQBZPpEshGDDdXf+t+sw9+7d\no0KFCjlua06hVNlUUPgXISEh7N+/n4SEBO49uMfBAwdRqVRoamoiaUjYF7en9bDWNOnVBD2Dd1eu\n/BQ5s+0MSwcvJGHeD/ltSrbxDAxml/t9Zp2+DICeiTH9929nXbN/qpx2qFCGg4N7ZWlclTqV0vN/\nY8KUqfTq3RtLS8sctTunUJy+gkI6x44fY8r0KXg99KJeh3qYWJlg6WhJ/S71sS1mm9/mFQiWjVzG\noVWHEIun5bcpmSYqIYlivyzDysgQfU0N7gWGAGBb3hnTIo60WzgLm7Jl8Nx/mKPf/EBMUAjq5GRs\nTY0JnP5tlua6FxDM8EOnueL1lMTExP8sZ51fKMLoCp89siwzfeZ01q5fy6jVo6javCq6+gVPDvBD\nOLT6EBsmriM5SYW5pSnNvmzJkDlD/rOPLMv80HIyPp7PSVGpUatTiY9NwKXwx5Wh1Hz1FqITEolX\np6JWqShatxZDzx9GS0fnjXYVOn9Bhc5fALCtx5fc23MAWZazFOqpWNiGFV80p2lgKCqVqkA6/eyg\nOH2Fj5qUlBSmzZiG9zNvQkND8XrkReFShVnuthwLW4v8Ni/HkWWZFaOW82UNFzpXdObAPS82L9iN\nzz0fZh6amaFTUyWp2Lt4L25nbzOsbjUsDQ1ehUM29e6Y17fwQRQyTNOx7bl9LeU6tc+UYLqBmRnG\nerrZqjbqUtiG1s7FWfzrr0ybOTPL/QsiSnhH4aNFCMG3E77lgtsFWv+vNaZWptgWs8WxjGN+m5Zr\nRIZE0t2mOykLpqCllebEXJ/50fy3rega6tH5m66kqlM5/ftJtPW0CXsZRlKSCh0tTSY1rcvMNk0B\niElKosHyTezo14Xydh9P7X2TyXOITVYxJdQbo0zE2dVqNT9qp7V7NHkkZWyyHpt38wug2+6juN/z\nzBN946ygxPQVPhvOXzjP+O/GExgSyMLLCz/Jg1QZsWrcKg4uP4B64ZQ3Po9JSmLc/pPs9LhPoiqF\nRiWK4mBmQuMSxWhfvjS2Jkb5ZHHOovHNDASgqaPNZF9PtPTSQncvb7rj1KQBWlpaxIWFcXXlei7O\nWURKsuqN/tnZv1Cnygw/cJL9nl6M/eZbyjg706NHj5y4nQ9GcfoKnzxubm7UqFED80Lm9PyxJ+2/\nbv9ZZd1813QC7uc9ePzDKEpZFXrreloqKq+eAj41ZFnGaPIcElPUb3ye7vwwMDYiITYOgL7VK7Ks\nU2tW/+XG8ituJKWkEPXLpGzNe+zBE7Q0JA4+fMpvl6/z5759dOrU6YPv50PJVacvSZIGcAvwTy+4\nVhlYDegBKaQVX3PLoF9rYAlpB8E2CCHmvWN8xekrvJcBAwawdetWVtxY8dFUuMxJOpp2JD4mnmWd\nWzO6Ya38NidfeB4eyaWnvvStWomA2Fhmn77Eok6tCIyOY6eHJ/cCQ2hUvCjD62ddhzgjph4/z8+n\nLr31eVBQEDY2+fuEmdtO/xugGmCS7vRPAgtfk0ucKIRo8q8+GsBjoBkQANwEegkhHmUwvuL0Fd7L\nvXv3aNS0EePWjaNep3r5bU6eIssyLTVb0r5cKQ7/r09+m/PJsvD8X0w/eYlUIRjbsBb9qlag88bd\nPAkNB0nC1MqR6BA/EhIS0NfXz1dbc63KZrpcYltg/Wsfy6QpaQGYAS8z6FoTeCKE8BVCpAA7gY8r\nXUChQFGxYkUO/HmAzRM3M7vnbCJD/pFevnrkKnMHzMXjggdbZmxhUouJzO47Ox+tzVnU6rSQRhFz\n0/e0VMgO9wNDqPrrGiYcOo2tS1OM7Usz98xlqixci3dYOMZmlvx4IgFN7bT0UB8fn/w1OJtkNmVz\nMfAd/zh5gG+Ak5IkLSRdLjGDfvbA6xKKL0j7IVBQyDYNGzbE844nP037iWGVhjHj0Ay2z9zG1aPX\ncLQ15czWM+jpapGSkkqqLPB74MuKmyszTO+7euQqF3ZeYOxvYzEwKtildHV0dHBp5MLqi250qVSW\nZqWL57dJHzUnH3oz68xlIhOTCY6LJyx9L6BI+br0nXcMgFDfh+yY3JaoIB86fL8NLR09KjXrxYUt\ns6hRuy5aGhJ/7t1Ls2ZN8/NWskSuyiVKktQVaCWE+LvWfj+gphBiTAbziGnT/tlZV5SzFDLDzJ9n\nMm1q2vfm7KoBNK3u9MYhnI7j/+DQ5cev2tfpUIcWA1rw5697eH7vOQnxSQBoamnStHcT2g//gvJ1\nyuf9jWSBMbVH8+D6Q4pbWTCohgvfN6ufrRz0z50ai9bhERROUZdGWBevRNX2Q7B0zLjU9L8PdsWE\nvuTmvqW47lzA6t/WMuzr/+WJzf9WzpoxY0auaOTOBvoBatLlEoH9pMklmr/WLloIYfqvvrWB6UKI\n1unvJwMio81cJaavkB2io6OxtraiTkVHvh9Qh1Z1Sr7VRqVSc9btOe2/+QNZCLS1NSlia8bo7tUJ\nDo9n/WF3EpNSiEtQYWFtxu7gvflwJ5lHlmVun73N8fXHcd3nipYksaV3J3pUKdg/VgWFS94+DP/z\nOA+DQrAu4cKw9R5ZHuPwgiHcPraB58+fU6xYsZw3MpPkespmNuQSNQEv0jZyA4EbQG8hxMMMxlac\nvkK2uHfvHqdPnWL8hAmM7FGHFRNaZrqvWdN5oK9LfFQ86hQ1bQa3Yfz6j6fqpCzLLB6yiBObTiKE\nwERfj+IWZkQlJXPwq55Usv88zi+8zsPgUPps249nYAjq1FS0NDVxH/8/KtjZEBQTh920hegamFCr\n2zfU6f4tekYmWZ5jzeCKJIT5ER0dnQt3kHny2unXA5YCmkASaT8A7pIk2QHrhBDt0/u0Tm/3d8rm\n3HeMrTh9hQ/Cw8ODVi2asemn1rSpm7l6+MaN51KukQuzjswi8Hkgdk52H2WoRJZl4qLi2LdkH1cP\n/kXA0wBSklTc/W4YZW2s8tu8PEOWZTTH/4y2rgGtRy/FsogzZ9dM4sXDawyrU401V29RyKEkwzd/\nmFTkvp97oxX9lNtuN3LI8uyhHM5S+OxxdXWlW+cvuLN1MDaF/vsUqkqlRrf+L7Qf1p5xq8flkYV5\nR0fTDlhq6+Dz09j8NiXXCYiOYeml66x0dSNepWL0Nm8s7Eu8ur7126Y8cz9Pvd6TaD40w3Xne0mI\nDiMqyIeXj25yfOlIli9fwciRI3LqFrKF4vQVFIDvJ33Hw5snODCv63vbmjWdS3RcMmfEmTywLG84\ntfkk879cgCRJTG3ZkOmtG+e3SbnO4J2H+P26O9p6hrT7ZhUuLQe81SYmLAATy8LZnmP1QGdC/Lxo\n36EzTRs3ZNy4sUhSlvxtjpNrefoKCh8TcXFxBITGkplFRHH7tFyEg6sO8uzus9w2LU+4tDetgqZq\n/k856vC3u92jxeotOMxYTJk5K5l6/Hy+yy2mpqbNP7ZB2unklKT4DB0+8EEOH6DbzP2YWjtw985t\nRo4cke8OP7sopZUVPhlSU1NZsngxJ48d4NbGLzP1n3LLtE50mbybVWNWIssy5lZmaOtokxSfSHxc\nIjZFbNDR16HlwJZ0n1Awimy9Dwu7tJLSJ7ye0L58ximImWXy4TMce/iEao6F2XLzDqYGBtQuVYqI\nuDh+PnWJK8/9mNe+BdWLfJhDzQ6zTl1iyvHzr95r6+rz1TLXXJvPqmhZev58kLVfV2Pr1m0MHjwo\n1+bKTRSnr/DJMKB/P3b8sZNHe0ZibJg54ZQKJW14vHc0AEFhcYxbfAIjfR3cHgZwJzKOl08DAFg3\nef1H4fRVKhVeN9M2KR1Ms56VAhASG8eLqFgiExKZd+4K1ZycOPXEl8IWFqwcNAgTg7RDbDP27OGy\nlxe1l6ynQmEb7IwNWdalTYaF4HIDHU3NV6+LVmpAnznH0DF49z6OKjGeqCAfdPSNMLMtmq05NTS1\nMDA0olOnj7ewgBLTV8h1kpOTWbVqFb1798bWNnekCsPDw7G0tKR9w7Ic/vXDnbNF8wVU61iXTqM7\ncf6P89TpUIfKjSvngKU5jyzLqFVqrh29xsrRK4kOiUKdmgrA+REDWXXlJld9X9CslBPzv2iOtbER\nIbFxjNl3Ap+IKB6GhKOjpYGNkSGBMXHEJqtISe8PcHjSJIz+QzVq4/nzPA0O5opX2o+Nka4uqbKM\npobExZFfUtXRLtfuvc6SDTxO1mL0H8/f2SbipTdrBldCQ1MDIyNjwkKCqNb2KyydKmFqXYSStdqg\nrZv5Gjp/fNuQ1b9Op2nT/D+Fq2zkKhQohBBcvHiRqVOnc/nyRVauXMmIEbmT7TBm9EjWrVvLwrEt\nqVDcmvVH7qFWy4zrWY2a5e2zPF7twb9zw9MfbR1tDE0McCjjyMjloyhZ+e3DXwDu59zZMHk9eoZ6\nBD0PItgvBH1DPdoNa08hu0LUaFMDOyc7dPR0MuyfVcICwpjd6xeCngURHR5NclJa3fiStrYs7N+f\nK15erDhxggSVCktjY6ISEtDU0AAh+KqmC7/f8EBXW4fiVlZUdnLiZUQEd319AWjh4sKAhg1fzaWT\nCXUqAA8fH+76+pIqy9hbWLDvxg38w0JpWaYEVext0dbQ4PsWDXLk/v/GYcYSdIpXZ8DCjDfiY8MC\nOLZ0JI9cDxAZGYmZmRm+vr7s3r2H575+bNq4gWbDFlPti6FA2nf23pkdRAd602DA27X31apkdkxs\ngYOpBhfOn8v39F7F6SsUGO7cuUPlypUpXLg4+voWREQ859tvv+H69ZucPn2CTp064+DggKfnQ1JT\n1ezbtwdjY+O3xklJSSEyMhJPT0+srKwwNjamaNGib8XrAwIC+HnGVMLDQnnxwp+q1WpSxrkcM2dM\nZd7Ixgz6Iuur9IQkFct23eD3Q+488Y/ghx0/0LT326u7rTO3snnaZgx1dXGytsbM0JDi1tb4h4fz\n15PHaGpokJCUDIC5lRmLLi/+IHWvq4eusmLMCkJ8g+lQowZGenr0q1+f6MREbEzfLMamUqtfOe3n\nISEsPHyYsNhYyjk4MLVbt2zbkBkSVCp+PXSI697eJCSn3X/qwik56ignHT7N/HN/0eXHbdiXrf0q\nTfPArJ68fHKH+PCX1Kpdl769ejJkyNsxeFdXVzp06UHrb9YC4H//Kn/9MRdZlqnfZzKpyfEkhPkT\nE+pPRLA/8TGR6BsYokpKJCQkBBOT7IXQcgrF6SsUGGJjYylevAQNGvyPMmXqcujQL+jpGePkVIuE\nhGhiYgLR0THC0rIo+/bNYtGiRYwZM/pV/5iYGBrWr8Odew8wNjKgallHXgRHER2XSFhEDM2aNubM\n2fPvNiCd48ePM2b4IB7vHpqtbIshvxxiw0F3arepyaxjb1fs/L7FJG6euUU5BwdWDh78znFCY2JQ\np6Yyc9+fPA4IZHfwbswszV5dVyWpuHvxLiWrlXzjc0gL35zbcY4Xj19wbM1RIkIi0dXWZmizZnSp\nVfBr6qvUamRZpsOC+Sxo35yxjWrn2NiyLFN+/m88Cg5FS1uHH08lc3bNRFx3LsDUzIyTJ05Q6z3/\nRr8uWsKRYyfQ0tJCR1ub2T9Pw83NjeDgYAwMDHB0dHz1x8bGBs3X9hLyG8XpK+Q5QghiYmIwNX27\n3K+HhweNGzdjxIht6Oi8O2YaGurL9u3jad++La1bt8Td3YO/XC9RvjAs/aY5iclqLM3SNg8jYxKx\naD4fJycn3N3dM5z3dZKTk6lUoSxz/1eDzk3KZunewqISsGq5AFsnW7Z4b3lrhRrkG0S/Yv2Y0rUr\nTStUyNSYoTEx9Fi8mEOxh3A/68657WfxvuVNsF8wanUqeoZ6HIk78qq922k3lo9Yzkvvl+hoa1HM\n2pqfOnfBsVDebJbmJJ1//ZWo+HisjY0w1NHG0cyEk0P7oafz4fkkfpHRFJ25hC4/7cDQzIprO+cQ\n8eIJ+jrPNZoTAAAgAElEQVSa6OrpM37cGL78ciBCCDw9PSlTpsx7vzsfA9lx+p9t9o5KpUIIga5u\n5rI8FODhw4dcunSJW7fcuXXLnejoaGJjY4iICKN48ZKoVCpq166Nrq4ulpYWnDlzlujoCJ4+daNs\n2XfHcq2sijJs2CZOnVqKm9tSLCyKc+3GTb4Y3pSA0FhKFfnHwZmb6OP5x3DafbsTMzMznjx5QsmS\nGcfZAS5fukRoWDjaWllfnVmY6FGqaCGePA/i8G+H6TjizYyNuxfvoqOtlWmHDxAdHw9AB+MOaGhI\nWJmaUsnBkS7NK2FhbMzMvXvpZNoRAyMDhCaEvQijiJUlY9u0oVPNj7sq+crBg5m8fTs2pqbo6+jg\n7uOD/qRf2NS7IwNrZn+TXJZlGq3YlPZGyBSv1pzi1ZqTmqLC754ryQkxrPpjI+MnTEAIgaV9caKC\n/bh54zrOzp+f+tqHyCXuBEqnXzYHIoUQVTPo5wNEkya6kiKEyPCbm5cr/VWrfmPUqBFoamrx3Xff\nM3v2jByf48SJE/z22zouXjyPr6/Pq9ifEIKkpKR8V9x5HSEE8+cvYPLkSdSoURstLW1CQ0NRqZJR\nqZIxMjKhbFlnrlxxpUSJWpiY2OPoWB49PSM0NDTQ0dEnOjoESFu1p6QkoVIlEhj4CH19I5o2/Roj\nI4tM2yPLajZtHE1omC9JSYloaWpwekV/GlcrBqSt9ot3WUnvPv1ZuGjRO/8tL128SKPGafH8DVOy\nnmInyzI9f/iTo9ef8WfEPnR03tyE3TJjC1umb6FRuXKMbdMGc6PMiY8/DwlBQ5KwMzd/a5PUNzQU\n10ePCIiM5EVEOBO+6PBRruozw9Hbt/n18GHW9GjPgGou2V7xL7l4jW8OnKTtuJXU6PjuRIHk+Bg0\ntXXQ0tFjRhOJCRO+Y8GC+dk1v0CQp3KJ/7r2KxAlhJiVQb9nQDUhROS/r/2rXZ44/Xv37lGpUiWG\nDLmOn58rly79RFhYKO7u7lSqVAkTExNCQkKwtrb+oHm++24iv/66AEjbjNTS0iIpKYkmTZpx7dpf\nVK9eC3NzC6ytC7Fp08YMBT7yijVr1jJt2mwqV26PmZkt2tp6GBlZoKWljaamNklJcYSG+qKpqUW5\nco3y1LY9e6bz4MFFNk/rxIB2LgA88QundLcVbN+2jZKlSlG9evUMNwenTpnCz7NmEX1uMiZGWX+i\nm7vZle9XnmXYomF0++btTU+f+z5832oyMRGx6GlocWDChKzf4GfMhQcPmLFnDwAmenpEz8meaLn+\nxNk4VG5K3/knMt1n6/imPLt9Hm9vb0qUKPH+DgWUXHP66XKJG4FfSBdR+dd1P6CJEOJpBn2fA9WF\nEOHvmSNXnL4QAj8/P27fvs3Nm7dZvnwplSuPolmz2Tx4sJc9e7q/amtpaUtYWBAAwcHBH+T4Hz58\nSKVKLtjaFiU09AVmZhbo6+sRE5PA11//jp/fPVSqRM6fX0d8fCS2tvb873+DsLQsRFhYGCdOnMbX\n15evvhrI5MmTciwMFR0d/UYsUwhB+/Yd0dIqS5UqbXJkjpxkzpy2qFSJxF74HiODf1baF2758POm\n65y/4cXkSRMZP+E7Cv1rRWxf2I4apc04sKBntuYeMusQO84+5PBrMfaM+HvFf/LHHzOd3qiQhkqt\nJjg6mgErVtCwRFF6Vi6PkY4OVR1tqWCXVhZ6zunLlLIqRLfK5d7qL8sy2hNmYW5TlBHbn2Y6Mygh\nOpwDs/vx5MaJV4uyj5HcdPp7SHP4pqSXVn7tWgPSBNLfFbZ5BkQBqcBaIcS6d7TLMacvhGDcuG+5\nccOdBw/uoaWlR+HCVSlUqAq2tlUpXrw5Ojr/PIoHBaWJKCQlRaGtbcCuXV8gywlERISjqamZrd36\nb78dz9mzHnTq9BMpKUkkJESTmBiLmZkN+vpvpnklJsbw4sVDjh37FXNzW/T0zHFwqMDdu8cJDvbh\nwIEDdOz4YScAo6OjGTZsBDt37qBUqbKYmpoQEBBAbGwMhoZm9OgxB3Pz3DtE8yHMmtWCUo5mXF0/\nGDOTNw8JeT4NYczCk/iEpbB5y3bq1av36j9+v379SAq6y945XbI1r9uDAGp+tQ5rByt6/dCHdkPb\nZehUPC54MKHJBLaMHImjpWW25vrcOXr7Nvtv3OBlRAQCSE5JoYSlBaFx8cSkp7uKxW/nzQP8cOQs\nc866Mn5fMEbm71+oJcVFs3NyK4KfezJp4iSmT5uSk7eSp+SK089ALnG8EOKL166vIk38fPE7+tsJ\nIQIlSbICTgOjhBBvFcjIabnEUqWc8fb2olu3XZQvn7UTmklJUcybZ46engGpqWo8PNwpV+7tVcZ/\nsWHD74wcOZLWrUdRuXK7LPUFePTIlV27/vkyxsfHY2CQPQ3XyMhIrKysKFmyOk2bfo0spxIfH4WJ\niSWGhuYYGBTcLAYfnzu4um7n6dObOBez4uHut2O2siz4/bA7i3fdJiQygbDwSFq1bMnly5eYMbQx\nE/pmP60xLCqBFqO38eB5KLIsMDAxoHLTyjTo1hArBytunb7F1hlbqVGyJPP79v2QW1V4jbCYGAat\nXk1sUhL1nZ1xffSIexOHvVr9/83fefrFXBoxcMmFTI2tSoxjTltjli9fzqhRo3LB+twjP+US9wkh\nBqQrY70EqgohAt47mSRNA2KFEIsyuJaj4Z3k5GT69BnIvn270NMzon37tZQv3zvT/VNSEpBlNbdu\nLeP06SkMHDiETZsyfEh5J2vXrmXFiu106ZL1jWJZTiUgwAtdXUM2bBhOixYt2Lp1M2ZmZu/so1ar\nAV49qsbGxjJ16jQSEuJZuzbt8MmUKWfQ0Cg4ecbv45dfWqJWp1Dc3oI5I5rSo8W7JQGFENx5Eszv\nh9y5fDeA5d+2oH7lIjlihyzLtBy9jbM3n2Nrbk5EXCzqVBlZltHV1ubEDz/kyDwK/yDLMmpZJkGl\novOCBVgZGRDy83dAmuRhty37CI2NxdSyMOP2vMzUmKnqFDyOb+TIoq9xdXWlXr16uXkLuU6eKmel\nv28NTBJCNHlHewNAQwgRJ0mSIXAKmCGEOJVB21yJ6aekpKCjo0O/ficpUSLzMnqvc+PGSo4fH0VA\nQAB2dpkLgcTExNC9e09Onz7J1KnnsjXv36jVKrZuHY+BgczDh/ffuObt7c3x48dxdb3K0aOHiY+P\nw8WlCp07d2b69Kk4O9fG2roUenpGRES8pFWrkWhp5UwpgNzmypWdnDmzhjkjmzF5YP18tSUmLolC\nLRbQvmo1xrZtm6+2fE4kqVScuHMHt6dPueLlhb62Nno6OkQlJFCqZmsaD5qFXem3kgaJDQsgNjwQ\nq6JliQl9ga6hCadWj+fxlYOUcXZm7i8/06ZNwdvDyir5kaffE/jjX0a8LpdoA+yXJEmkz7U9I4ef\nm2hra9O8eVu8vPbg5NQsy6tctTqZ48dH0axZc2xtbYmNjWXdunXUrVsXXV1dqlSp8lYfIQRly5Yj\nIOAl48f/+cH3oKWlg7V1EdzcjtGzZx+EEAgB2tpa/PHHNqpUaYmDQyW6d/+FTZvGEhISzPz58+nd\nezalS9f54Pnzitu3j3L37iliY4KIT4hFnZJM3zaV8t3hA0TFJaNOlXHJRxHsz43nISEMWr361XsL\nOyfMCpekqEsjilVuRJGK7/5ebP22EaH+3q/eGxqbEB8bw9q1axk8eHC+18zJTz7pE7lCCDZu3Mzx\n4yc4duwEY8b4o6v7dn2X93HlylxevNiGqWkhwsICePrUG0nSQAgZlUqFtrb2q7bPnz+nQYNGvHzp\nzxdfTKBq1azH8zMiISEaHx8PkpPTDvdoaGghhIyxcSFKlEjTo09NVePpeQ5NTS3s7EpTqJBDjsyd\nF5w6tZqrV3dT0tGC8k5WpKSmMm9kcyqULDjC3pV6r0adpM2qIUPy25RPnpazZpGSmoqBsTlGlva0\nH78Wx/KZW8CkqlNYN7QqDWuUw8LcgmnTpmb6Cf1jQzmR+y8WLVrCkiWbqFBhMAMH/pQthw9Qrdow\nvL0PU7bs/7CxSaF9+ybs2tWWXr3avOHwHz16RMeOnSlatDaDBm3M0di5gYHpe3PkNTW1cHHJXggr\nL1CrVdy9expHxwpYWf1Tzzw83P/V601TO1LPJWfi8DnN+L51GPTzISLj4jJ9EEshe1iamBAYGcng\n1Tff0LrNDMeWjCT4uScVv+rBlCkfb2ZObvHJrvSvX79O69bt6dv3MpaWOXvU+uXLG6xfX4vAwEBs\nbW1JSUnhypUrdO/eg7JlW9KkyeCPVkrtQ5FlmUePXHF2rk94uD/x8VFcvryNwoVLc8fjOLFxaWf0\njAxNMLdwJDT0GUlJiQBYWxjxZO9ITIzeXbs9v5FqzmBAw4Z81STDbSyFHEKlVtN2zhwKFSnL8I2e\nWeorp6Zy/8JuTi8fie/zp5ibm+eSlfmPstJPZ/v2HYwePZ7mzVfmuMOPiwvi4MFeFC5chLt377J0\n6TLWrVuPvr4pZcs2o1GjgZ+tw09KimPt2sFERoagoaGBLKfplxob6BLw8g7F7c05v3c8sgy/bLzE\nvachtKrqzPLxbZGRMcihWvO5iSTB5UePqF6iBCb6+tiZm3PCw4OStrb8ccUV10de1Hd25uee2TsQ\nppCGjpYWFR0d8fC5jyzLWYrBa2hqUrFZbzxPbWLZ8hVMm6qs9l/nk1vpb9u2jf79+9Oq1RJq1x6b\nA5a9yeXLszl37sdX72vW7ICzcyOcnN7OIPjcOHBgLnfunGTdD18QFp3AiK7VC/SqPTvM2+zK3C1/\nEROXhJz+fdXU0CBVflMgvF+DBgwuAMpKHzNxSUl0WrAA50bd6TZ1Z6b7qRLjuHdmB67bZ9G9UzvW\n/rb6/Z0+UpTSyqRt3i5fvoKxY8cwcWI4+vqZL/SVGdTqJNavr0NwsAcVKzanS5cf39/pMyEpKY41\na4aiLcUQcvLTr0PzxC+cH1ad48dB9dHV1sLJzpxktZrWY7ZxzfMlG4YNo7hNwdmI/hjptWQJEUkq\nfjiRkOk+j1wPsGtKZ0aP/YafZ0z7JEoov4vsOP1PLm9JkiQaNkwr43v79socH19LS48+fY4C4Oz8\ncR/syGlkWU1UVCBLv22V36bkCaWKFGLP3O5ULm1HWScr9PS0MDXSY2B6YbiVJ08Sk5B5Z6XwNoXN\nzUlJTsxSn1sHV+JY1IllSxZ90g4/u3xyTh+gcuXK1K1bH2/vk7kyvptb2uNisWIFUyg7v3j27DYA\nvVtVzGdL8pdhXWtwfvVAbj9/TscFC3jw4kV+m/TRMjG95tTNA6sy3efFwxv4+z5nx44/3t/4M+ST\ndPqQVobBx+cKBw58SWJiRI6OXb36cABu3TqWo+N+7JQunSaD5+ru9962N++/pOOEnfT4fg+9f9xL\nrx/3Yt5sAQ2GbiQhXeT7Y6ZxtWIcmN8TYwNdRv2+gf4rVuDh45PfZn102JqZUd7RkWu7F2S6z7jd\nLyjfqCujx4xm7dq1HD58OBct/Pj4JLN3ANzcbuDt7U2pUqUwMrKiefPMf2n+ixcvrrNtW1r44r8k\nAD9HdHQMcLAvRfPRW0m49P07My5kWab+0E0YG1uiVquQNDSRU1Owsa3EzQcemDebT4papqSDBWN6\n1sTe2pjVe93wDY4BBF57Po4iWTo6mliZGxKfqOJFeDjfbN7M+WkZV4pUeDflHRx46nEv0+11DYwp\nXqUpCWH+fP311wDvVVj7nMgL5azWwBLSnio2CCHmvWP8XKm9U6VKdRISbOjd++gHj/XkyTF27GiH\ngYEZ9er1om5dJS3v36jVKubOacvQzlUY1b0mL0JiqFvJESMDHTYcvM3CHdd57BuWplo28fBbdYBS\nUpK5fn0fpqbW3Lp1CF/fu6+u1apTlOtXfSlTzJLihc1QpaQyf3QLbAsZYlvIKMtH62/ef8nqfW5U\nKW3LnSfBeD4LoVZ5e5aO//CaLHM2XuaH1eeo6uSES9GiePr707xiRVq6uHzw2J8TkXFxdFm4kNrd\nxtFqZIaFfDMkOSGWue3SSpgvXryYESNGvKV89imQ23n6Y4H7gAmAEKLXaxP/SlrN/H8bpAGsAJoB\nAcBNSZIOCiEeZcXI7PLixQs8PG5RqVIPgoLuUKhQKbS1s1eeGODQoa+wtS3J119nrdrm50BiYiyR\nkQHcuLGLVDmVNftvccDVhxLFnXAdsw1dbU1UahlLS3s6dvqB8uUboaHx9tdPW1uX+vXTqqFWrNgM\nf39PDh1aQFiYH1Y2xvwwowV33F/idT8UIQTVBqx91VdDQ2Jo52qsnvT+0hdP/MKp+dV67B3M2HX2\nPkZGejiXt2bZrhv8ef4RwRFxSEhYmOpT0tGCjg1KM7xrDXS0NAiKiKeI7X9vENpbmyBJUNbengGN\nMj5JLcsy/uHhFDI2Zuafe1GnprJowMD32v45cfjWLQCaDP4lS/1SU9JChPUbNKBly5afpMPPLply\n+unKWW1JV87KoEkPIKMjijVJq7Xvmz7OTqAjkCdOPzg4GAAvr8OEhV0nIMCXChV6ULx4K/T1Ld4S\nU3kflpbl8PG5gEqVgI5O9n88PnbUahX+/p68fPmIwEBPNDW18fW9h46ONj169GD69DE4OTlRvnx5\n3NzcqF+/AZrahnzV/xccHd9dGjkjHB0rMHLkZvz9Pdm4cTQR4fE8exrO4XNDKVXGitCQODS1NIiN\nTeLo/vv88M0RgsJiiU9KYdAXVejVMk20XJZlNh+9w3XPFzz0DeeJf9o+z6mrw7F3+Kdc9aljD1m3\n4io/9qyEubkBVy4/x+0vX2b8fomJy8+8amdnZUzz6sXo3KQsnRuXfcvuCiWs0ZCkV7n8/yZBpeKb\nTZt4HBgIgKWlEZGRCWw8f57udepgpPfh5xvu+Phw+dEjQmNjuOn9lL7169O3wbsF6gsiLV1c2HTh\nAsv7ODF+X3Cm+gghePHgGjYlKuF6+TK9+/TjjsftXLb04yFXlbMkSeoKtBJCDE1/3w+oKYQYk0Hb\nXAnvhIamrQitra25desWR48e4+pVN8LCwvD2fkqFCgPw9j6Ks3MXGjX6+Z3jxMWFsHChDbq6Rkyc\nePCjrtIXEvKclJQkbG1LoqmpTWpqCgEBXqjVKTg5vV019G+SkxPw9DyLu/shTEz0KVzYlh49ulGo\nUCFcXFxexUzv37/P6NHjSEhI4MGD+9St24fatbMmZJMRf5daBth99Etatn3b2U6ffIyt691IVqmJ\ni03GxsKI0Kh4ZFmgqSFR2MEMe0dTSjpbUatOEfp+lbG+bkY8fxqOnr422joadGq2Hs+7gZgZ6xF5\n9m1tV8OGcyhn78CCfv3f+NwrIIDgqChmHziArr4WI76pR+eelShT1oY+HTdz7NADgFc/GNpamhjp\n6bH8q0HYW1iw4sQJDt26RWpqKt1q10ZDkihsYYGpvj7HPTyISUwkLjkJv9AwAKysjTHQ1yYuPpnw\nsHgqFHFkXt9+GHxEK9/1586x/fJl2o5bRY2Ow9+45n//KiaW9ugamhAT+gKrYuWJCfFnSa+02k6N\nGjfm0MGDmJiYZDT0R0+BU84qCE7/v3j48CGbNm3l3LnTuLm5MXbsc8zMimXYVpZlVq8uS1jYYwYM\nWJjhCVx392P4+3sQHR1I06YjsLd/2ynlBbGx4fj63iEqKoC4uAgMDc2pVasbOjr63L59iMOHF1Oq\nlDNxccmUKdOIW7cO4OjoiJ+fD4ULl0ZTUwtjYyuaNx+Fv78njx9fIjVVzcOHrjRp0piBA/vTtWvX\nd5abaNfuC27cuEvz5kOxsyuNiYlVjtzXvXtn2LfvF7bs7UeHru9PC/3L9Tkrf72Eta0xXw2rRcnS\nVhgY5JyzW7HoMlO/O0pJx0KYGOpwdkV/TIz0mLbmPDM3XOLo99+/4Vz3XrvGypMn0dCQaNikJPtO\nDXrjB0eWZeLiVERHJRL4Mhora2NuXPXl6/67sDAyIiIuDoAiRc1JTRXExSajq6tFSEgsAOUq2GFl\nY4iVtRHx8cksWdMFG9t/nN3p416M/HIPYWFxWJuaUrdUaUZ/BDXlZVmm2c9pC7JhG+5hUzzt6S0l\nKYHZbQzfaCtJEgZGJsTHRr/6bPHiJYwbl/On8wsCuRXTrwd0kCSpLenKWZIkbXlNOasL8K4aBC+B\n10smOqR/liHTp09/9fpD5RIzQ9myZZk3bzb799ekT58B6Oq+O06roaGBtnbaF8zY+G0n9vjxVW7d\n2sWPP05GT0+PQYMG0afPXEqWrEHa1kbuERj4hCtXthATE4osq4mICKRJk2YUK2aLiUkxnjx5ypo1\nX2JhYUdcXChXrlyhbt267Nixg7/+usbMmYdo0KABSUlJnDp1Ci0tLX77bR3r1g0mOTmeqVN/QldX\nlwYNZuOSiY3I+vXrcuzYEUqUqJGjgi1FilREkiRcLz7LlNOvW9+JuvWdcmz+fzPq2wbY2Bmxf9dd\nzhz3ot7QjYzqVoOZGy7Ru169Vw7/6O3b7LtxnWfBIVSobIer+7gMx9PQ0MDERA8TEz0ci6QVCXMq\nUYht69144BnE8CH1KFnaisHD3ywxHBWVQMCLaMpV+O/ywS3alOFR4A+cPPKIIwfu88fmm1ibmtKz\nbt0c+NfIPX7cuRMNTU1K1GiNsWVhAGLDAzk4dwB29o543HYjMjISR0dHVCoVL1++5MWLF5ibm2Nn\nZ4e9vX0+30HO8W+5xOyQ28pZmoAXaRu5gcANoLcQ4mEGbfN8pf83s2bNYsqUKYwc+fCdBdrOnv2R\nK1fm0LXrFMqXT7vdwMAnpKamoFar2Lz5G/bt20fnzp2RZZnhw0dy9uw5EhJUtG37HQ4OWdPYzQwq\nVSLPnt3iypU/iIl5wYYNG3BycsLZ2RlDw39WQLIsc+vWLaKjo6lXrx76+u9PNU1NTeX+/fvY29tT\nqFChLNn18uVLmjdvgUqlTf/+S7N8X//FhQubuHhxM5bWRrh5jcfMrGDsrSxfeJEpE9LObXStVYtR\nrVtz8OZN9t64zouwcGxsjRn9XSOGjqqDjk7ByJSe8t1Rflt8hcUDB1Le0TG/zcmQxUeOcOjWLfrM\nOUqp2m2Jiwji+K+DeHD1OJUquTBv3lxat26d32bmG/khl7gRuCqEWPtam9eVs/7+YVjKPymbc98x\ndr45fW9vbwYN+h+PHj3G2NgRXV0jChdugbNzRywtnVGp4pgzx5iqVdvxxRdpNWXSavz0Ji4ukpQU\nFZaWVgQFBaKp+U8NfSEEVapUo2jRFlSpkvOP0du2fUtEhC9t27bj559n4OSUe6varOLi4sLdu3dz\nRb3r6VM3du78HrVaTcMmpTh0rmCImvRov5HTxx6hr6NDcoqaVFmmanUHOnavyP9G1c3R0FJOIMsy\nXVr+zsVz3qz/+usCVSdILcuM37KFu75+lKjenH4LThEd4s+p5aMoY6NDvz696dChwys96M8VpeDa\nB6JSqbh48SLPnj1jx459eHjcpFy5njRvvpRly0qgoaGmb9/5WFo68uLFAzZsGElKSgpCiDfEVP7m\n2rVrNGzYiO++O4S2tu5b14OCvLG0LJKtEIhKlcicOW25ePEiDRs2zNb95iYnTpzg3LlzLFiwgBIl\nqiNJKdjbV6Jx40E5Nsfx48u5cWMfURkf/cgX7tx+ybpVVylT1opO3Su9CtMUVO56vKRhlWWsGDSo\nwKz21bJM32XLCYuLo8uP2ynXuDu+dy5yZG5/vmjXmiWLflVq6qSjOP0cJjo6mhIlSpGUpKZSpf5c\nvboMc/PC2NkVxc/vPnp6eoSGvjuN7MSJE7Rp04Y+feagp2eMEDKRkS958uQywcHPSUqKIy4ulv79\nf6V48WpZsu3UqeUUKqTmwIF9bzxdFDTOnj3LunUb6Ny5I7169aJVq5HUrt0tR8ZWqRKYM6cdg0fU\nZtqcNpiYfFplnPOC//X9g6N/PuDo5O/z25RXnL9/n5///JNRW5/w4sFVLv3+I9qagp+nT2PIkMH5\nbV6BQhFRyWFMTU25f/8ely5dYsyYb3FyKo2BgR4DBnRm0KBDGBj8dzy5Tp06DBkyhNu39+Hv70to\naAiVK1dj2rRJVKtWjWLFijF27LesWjWRYsXKUadOf0qWfCvzNUOuXt3H/PnzC7TDB2jWrBnNmjUD\nIDExka+++oqqVdvlSAkLHR0DGjYcwI6Nf7J7mwe/7+xDizZlPnjczwVZlnnhF0VicsGqdTRz714A\nlvcriZ29Iwtmz6Rv3z7KAascQlnp5yHJycno6r4d5jl27Bjt2rWjcuWWNGgwAAuL/8422L37Bx4+\nvMquXbvp0aN7bpmb4zRv3gI/vxj69MnZcIxarWLjxpEEBj6lkKURqzd3V5x/Jvh9zTW+G3GQb9q1\no13VgiMC1GTGDEaMGMHQoUMpV65chqFThTSUlX4BJyOHD3D9+k0AHjy4yIMHFylbtj46OvoUK1ad\n2NhwYmOD0NIywMqqGE+eXMPf35Po6OiP7sCJr68f3t6P2bHjGxwcKtOgwYAckZaMigpGT8+E+vX7\ncPnydlYuuqw4/Uywa0vaKdWLDx5QuVgx7C1yVnAou5QpWpRevXplKj1YIesoK/0CQmJiIjo6Onh5\neXH8+HF8fX05dOgoVapUpl69OoSHR+DhcZe4uFimTZtC8+bN89vkLBMXF0d0dDReXl4MGjSEpk2/\noUiRzNfeDwjw4tatI/j7uxMTE4YsBFqa2iQmxmPvaEZYaDwlS1tx/ubIApMWWZCJiEhgzTJX5s04\nS62SJZnbt29+m4QQgk4LF3L7zh1KlCiR3+YUeJSNXIWPhgULFrBx40E6d/4Jbe3MbcD+8ktLQKbN\nF2VxrmCDn08Exsa6TJ/XBqNPTIs3L2nTYDVXXX1yreyzEOKNJzpZCNyePuWevz+psoypvj4VixSh\nrL09kiQxctMm6rdsyZy5c7EoIE8fBRUlvKPw0fD1119z7dpNVq8eQMWKLXB0dKFoURe0tHQyDPls\n2gxEnEYAAAyTSURBVDQatTqFwSNqs/D/7d19dFT1ncfx95dMnhMyCQ8JKxCSQIClHrBSgRaUqkUq\nZ0Vxj1sfuoKnW6mVZekq6rZnqa5nK0hdu7u6ttuWtbtSRJCIK7Ho4oHoNjyYIAFJICEoCIEkkOeE\nPMx3/5gLO4SQBEzmDsz3dU5OZu7MnXwyM/nOze/e3/2+dJcLia9e8YnRxET17bh5+cmTnGpo4I38\nfLYfPMjo4cPJTE2lpqmJkzU1JKWkcMdddxEbF8cXR47wQm4u3shI7p48mQenT2f5mjWUlZby/pYt\nfZrL2Ja+cdnevXt58cV/Ji8vjwMHivF6BzFt2n1MnHgbp04do66ukrS0LF588TuMTE8m/9Mfhdwk\npytdWuyPaWlpZ+V3v8s1KSlEeTzUNDYyNCmpy7N9dvh85753+HzEBhxVc/D4cZ7NyeGIc6LD9BEj\nePmVV/B6vZSWlpKamkpSUhJTpkw578O9o6ODl196iQ3r1rFn717GpqWRX1JCfX19j0fJhTMb3jFX\nLFWlrq6Ow4cPs2TJY+TlbSUzczR1dbXU1JzG4/Hg0zbuuX8SZQdqGZoWz4KF1/OVicPweq2D2aVq\nbm6juqqRxIHRXJf1PKeqG/FERBAfF0djUxPD0tI4VlHB1HHjGOH1UtnURFlFBZU1NTQ2+xuVR0VG\n0trWxnemT+eeqVNJjo/nX3JzeXPHDtra2i57tmxDQwOLFy2itLSUrXl5fflrX3Ws6Jurhs/nO3cG\nyrNjwuXl5fzsuWe59ZbbOF5xnCeWLsWbHMszz88mJjaS2ppm0jNSyBw9iIFJMXg8EcTEeGhoOENS\nkn0wVFU2sPa1QrZ/9Dn/8+4Bzpxpp62tnQfnP8CK5T+nsrKSYcOGER8fT3R0NLW1teTm5lJWVsao\nUaPIysoiLS2NoUOHEhERQXR0NNu3b+fuefM4XV3N+PR0Cg4eZOljj/HcihVu/7phoV+LvtMFaxdw\nNODcO4uAR4B24B1VfbKL9Q4DtYAPaOvqvPvO/azom0ty9OhRHl74kHMct5KQkEhBQQHVVadoamqh\nra0Nn0+JjPQwZuwwvjIxleiYASxcPI3scUMveLz29g72FVWQNWYwCQldH157paqvP8O0Cb9g/PhJ\n7Nq5i+997/uM/9NxvPHG67y9cdOXmuRXVlbG22+/TVpaGnPmzCExMbEPk5vu9HfRXwJcDwx0euR+\nE3gKuF1V20VksKpWdbHeIeB6VT3dw+Nb0Td9qra2ltjYWESEnTt3sm/fPj76KI+cnA3ExUdx/4JJ\n1Na0kpAYhc+n/GLFB+fWvX/BZLLGDOZQaQ3XDB/IN2aOYviIJDKyBvXJ3IJge+gv1hIfM5ZX/+O/\nrsj8pmv9VvSddomrcNolOkX/deCXqtrt7nURKQcmq2p1D/ezom+CoqOjg61bt7Ju/VoyMrKor6/j\n+PGjLJj/V0yYMIEtW7ZQXFJMbe1pfD7Y+NabJHkT2bF9N7Gx0XiTE5g1Zyz3PHAtSd5YEhKjGZVx\n/qGFra3tVFU2MjApxvX/Gv43r5zbb3yFqqqqSz5Ntglt/Vn0L2iXKCKFwFvAbKAZeFxVd3Wx7iH8\nTdM7gF+papddxa3om1C3f/9+Bg0aRH19Pb9fs5oNG9bS2trKyRNVNDY2kZrmZfLUEbyTU4THE4HH\nE0lzcwvXTU5n4d/cwG1zxhEZGeFv3RjRf411mppa2fxOMXsKj5HkjSF3437GjZnBqlW/67efadwR\nrHaJZ7f0i4AtqrpYRL4GvK6qmV2sP0xVj4vIEOA94FFV/bCL++mygMkhweicZUxf6OjooLm5mfz8\nfIqLi5k7dy5RUVGkpqbS0tLChpwNLF/+LPX1p0hOieOTws8YMXIIN96SwX3zJzHl66P6LEt7eweP\nzF/H2tcKGDMmi+kzvsHgwUNY+PAjZGZe8OdprjCdO2c9/fTT/VL0/xF4AP/O2lggEXgTGAwsV9Wt\nzv1KgSndDeOIyDKgXlVf6OI229I3V7Vt27ZRU1PDrFmzKC8vZ/2b63nh5yvYWvhD4uIiqa87w4h0\nL7U1LTQ2tjIyPRlV5URFPaUHqqiva+Fb3x6Lx3PhTte2tg6KPz3BjEn+TmVnW2Kaq1tQO2eJyMPA\nn6jqMhHJBt5T1fRO948DBqhqg4jEA5uBp1V1cxePbUXfhJ1n/mEZK59fyYABA1AFj0dob/cRGRlJ\nckoslSfriIqKInvsaD76cCfXTR7Okz+9lZ35R7lh2ghm3jqa22/8DYW7PiMj4xpKSz+npKSE7Oxs\nt381EwTBLvqRwG+BScAZZ/nWwHaJIpIBbAAU/ykfXgvFdonGuKm6upqoqCgSEhI4dOgQ6enptLS0\nUFRURHZ29rmdrx9//DG//NVLbHpnE9+8+WZyN+Xy1DMzeeyHObS3t4d8bwXT92xyljFh5GxntszM\ndMrKDrsdx7jAir4xYURVWb9+PTfddBNDhgxxO45xgRV9Y4wJI5dT9PvvYGFjjDEhx4q+McaEESv6\nxhgTRqzoG2NMGLGib4wxYcSKvjHGhBEr+sYYE0Z6XfRFZICIFIjIxoBli0Rkv4gUicjFTq8wW0SK\nReSAiDzRF6GNMcZcnkvZ0l8MfHr2itM568+Aa1X1WmBl5xWcFov/CtwGTADuFZFxXypxEAWewjRU\nWKbeCcVMEJq5LFPvhGKmy9Grou90zrod+HXA4oXAc6raDtBVq0TgBuCgqn6mqm3AGmDul4scPKH4\nIlum3gnFTBCauSxT74RipsvR2y39fwIex3+2zLOygRtFJF9EPhCRyV2sdw1wJOD6UWeZMcYYF/RY\n9J3OWSdUdTcQeI4HD5CsqlOBpcDa/olojDGmr/Rr5ywRmQr8VFVnO9efBFRVl3fxc+xsa8YYc4lC\nrXNWBFAC3AIcB3YA96rq/ksJaYwxpm98meP0fwtkOg3SVwN/Cf5G6CLy3wCq2gE8ir9N4j5gjRV8\nY4xxT8icT98YY0z/c31GrogsdiZ3FYnIX7uY4zcickJE9gQsSxaRzSJSIiJ/EJGkEMj05yKyV0Q6\nROSrwczTTaYVziS93SKyXkQGhkCmZ0TkExEpFJF3RSTN7UwBt/2tiPhEJCWYmS6WS0SWichRZ/Jl\ngYjMdjuTs7zHyZ/BzCQiawKeo3IRKQiBTBNF5I/O+3zHRY6iPJ+quvaFf8LWHiAaiMA/DJTpUpbp\n+Ju87wlYthxY6lx+Av+8BLczjQXGAFuAr4bI83QrMMC5/BzwsxDIlBBweRHwb25ncpYPB94FyoGU\nEHn9lgE/CnaWHjLNdOqBx7k+2O1MnW5fCfzE7UzAH4BZzuVvAx/09Dhub+mPB7ar6hn1j/9vA+a5\nEURVPwROd1o8F3jVufwqcKfbmVS1RFUPcv7hs25nel9Vfc7VfPyFze1MDQFX4wEfQXSR9xP8/5wX\nV3STy5X3E1w00w/oefJnsDMFugf4fZDiABfN5APOjkB4gS96ehy3i/5eYIYzjBKHf9bvCJczBRqq\nqicAVLUCGOpynivBQ0Cu2yEARORZEfkcuA/4+xDIcwdwRFWL3M7ShUed4blfB3sY8yJ6M/nTFSIy\nA6hQ1TK3swBLgJXO+3wF8FRPK7ha9FW1GP8QynvAJqAQ6HAzUw9sr3c3ROTHQJuqrnY7C4Cq/kRV\nRwKv4R/icY2IxAJ/h38o5dxil+J09jL+YdVJQAXwgst5ILQnf95LkLfyu/EDYLHzPl+C/6jKbrm9\npY+qrlLVyao6E6gBDrgcKdAJEUkFcHYEnnQ5T8gSkfn4/1O7z+UoXVkN3O1yhixgFPCJiJTjHwL7\nWERc/+9RVSvVGRQG/h34mpt5HEfwTwJFVXcCPhEZ5G6kc3OP5gGvu53F8aCq5gCo6jr85zvrlutF\nX0SGON9HAnfh/wN1LQ7nb31tBOY7lx8E3gp2IC7M1Pk2N5yXyTna43HgDlU9EyKZRgfcdifgxvyQ\nc5lUda+qpqlqpqpm4D8P1XWq6saGROfnKvDIpnn4h12DrfP7PAe4GcCZ/BmpAbP9XcoE8C1gv6oe\nC3KWszpn+sKZNIuI3EJvNpqDuff5Inukt+F/kxUCM13MsRo4BpwBPgcWAMnA+/hnFW8GvCGQ6U78\nW0HN+Gc554ZApoPAZ0CB8/VyCGRaBxQBu/F/WA9zO1On2w/hztE7XT1Xv8N/FN1u/MU2NQQyeYD/\ndF7DXcBNbmdylq8Cvh/s162b5+nrzvNTCPwR/4ZEt49jk7OMMSaMuD68Y4wxJnis6BtjTBixom+M\nMWHEir4xxoQRK/rGGBNGrOgbY0wYsaJvjDFhxIq+McaEkf8DIqDwGKJMiLsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1196e4e90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gp\n",
    "%pylab inline\n",
    "austria_shp = gp.read_file('austria.shp')\n",
    "austria_shp.plot()\n",
    "austria = pd.read_csv('http://dl.dropbox.com/u/8649795/AT_Austria.csv')\n",
    "austria.drop(['Oi2007', 'Dj2007', 'OrigAT11', 'OrigAT12', 'OrigAT13', 'OrigAT21', 'OrigAT22', 'OrigAT31', 'OrigAT32', 'OrigAT33', 'OrigAT34', 'DestAT11', 'DestAT12', 'DestAT13', 'DestAT21', 'DestAT22', 'DestAT31', 'DestAT32', 'DestAT33', 'DestAT34', 'beta', 'Offset'], axis=1, inplace=True)\n",
    "austria.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **Origin** and **Destination** columns refer to the origin, $i$, and destination, $j$, location labels, the **Data** column is the number of flows, the **Oi** and **Dj** columns are the number of total out-flows and total in-flows, respectively, and the **Dij** column is the distance between $i$ and $j$. In this case we use the total out-flow and total in-flow as variables to describe how emissive an origin is and how attractive a destination is. If we want a more informative model we can replace these with application specific variables that pertain to different hypotheses. Next, lets format the data into arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "austria = austria[austria['Origin'] != austria['Destination']]\n",
    "flows = austria['Data'].values\n",
    "Oi = austria['Oi'].values\n",
    "Dj = austria['Dj'].values\n",
    "Dij = austria['Dij'].values\n",
    "Origin = austria['Origin'].values\n",
    "Destination = austria['Destination'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Oi and Dj vectors need not be $n^2 \\times 1$ vectors. In fact, they can be $n^2 \\times k$ where $k$ is the number of variables that are being used to describe origin or desitnation factors associated with flows. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Calibrating the models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, lets import the main SpInt functions and calibrate some models. The main SpInt functions are found within the gravity namespace of the SpInt module ans the estimated parameters can be accessed via the **params** attribute of a successfully instantiated spatial interaction model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pysal.contrib.spint.gravity import Gravity, Production, Attraction, Doubly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unconstrained (basic gravity) model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.44314667  0.51739961 -0.00979932]\n"
     ]
    }
   ],
   "source": [
    "gravity = Gravity(flows, Oi, Dj, Dij, 'exp')\n",
    "print gravity.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Production-constrained model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.90285448 -0.0072617 ]\n"
     ]
    }
   ],
   "source": [
    "production = Production(flows, Origin, Dj, Dij, 'exp')\n",
    "print production.params[-2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Attraction-constrained model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.90037216 -0.00695034]\n"
     ]
    }
   ],
   "source": [
    "attraction = Attraction(flows, Destination, Oi, Dij, 'exp')\n",
    "print attraction.params[-2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doubly-constrained model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.00791533]\n"
     ]
    }
   ],
   "source": [
    "doubly = Doubly(flows, Origin, Destination, Dij, 'exp')\n",
    "print doubly.params[-1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that for the constrained models we have limited the params attribute to print only those associated variables (i.e., not fixed effects), though it is still possible to access the fixed effect parameters too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.16851884  0.52128801  0.98284063 -0.56934181 -0.28515686  0.0381801\n",
      " -0.47906115 -0.0141766  -0.1583821   0.90285448 -0.0072617 ]\n"
     ]
    }
   ],
   "source": [
    "print production.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first parameter is always the overall intercept with the subsequent 8 parameters representing the fixed effects in this case. You might ask, \"why not 9 fixed effects for the 9 different municipalities?\". Due to the coding scheme used in SpInt, and many popular statistical programming languages, you would use $n - 1$ binary indicator variables in the design matrix to include the fixed effects for all 9 municipalities in the model. While the non-zero entries in these columns of the design matrix indicate which rows are associated with which miunicipality, where a row has all zero entries then refers to the $nth$ municipality that has been left out. In Spint this is always the first origin or destination for the production-constrained and attraction-constrained models. For the doubly-cosntrained model, both the first origin and the first destination are left out. In terms of interpetting the parameters, these dropped locations are assumed to be 0. Since the fixed effects parameters are interpretted as deviations from the overall intercept, this essentially means the intercept acts as the fixed effect for the first location. Therefore, we could also say first 9 parameters are the origin fixed effect parameters and that the last two parameters are for destination attractiveness and distance, respectively. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Interpretting the parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Assessing model fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.812858632041\n",
      "0.910155702595\n",
      "0.909354566583\n",
      "0.943539912198\n",
      "\n",
      "0.834893355219\n",
      "0.913167247331\n",
      "0.912363460094\n",
      "0.946661920896\n",
      "\n",
      "0.68757357636\n",
      "0.740913740348\n",
      "0.752155260541\n",
      "0.811852110904\n",
      "\n",
      "1.01757072219\n",
      "0.464519773826\n",
      "0.58404798182\n",
      "0.37928618533\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print gravity.pseudoR2\n",
    "print production.pseudoR2\n",
    "print attraction.pseudoR2\n",
    "print doubly.pseudoR2\n",
    "print ''\n",
    "print gravity.D2\n",
    "print production.D2\n",
    "print attraction.D2\n",
    "print doubly.D2\n",
    "print \"\"\n",
    "print gravity.SSI\n",
    "print production.SSI\n",
    "print attraction.SSI\n",
    "print doubly.SSI\n",
    "print \"\"\n",
    "print gravity.SRMSE\n",
    "print production.SRMSE\n",
    "print attraction.SRMSE\n",
    "print doubly.SRMSE\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 Local models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.6 Testing for overdispersion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Additional functionality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Existing features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Future additions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, there are several paradigms for incorporating spatial effects into spatial interaction models (competing destinations, spatial autoregressive, eigenvector spatial filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python (spint)",
   "language": "python",
   "name": "spint"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

SpIntPrimer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A primer for working with the <u>*Sp*</u>atial <u>*Int*</u>eraction modeling (SpInt) module in the python spatial analysis library (PySAL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spatial interaction modeling involves the analysis of flows from an origin to a destination either over physical space (i.e., migration) or through abstract space (i.e., telecommunication). Despite being a fundamental technique to regional science (citation), there is relatively little software available to carry out spatial interaction modeling and the analysis of flow data. Therefore, the purpose of this primer is to provide an overview of the recently develped spatial interaction modeling (SpInt) module of the python spatial analysis library (PySAL). First, the current framework of the module will be highlighted. Next, the main functionality of the module will be illustrated with a common regional science example, migration flows, with a dataset previously used for spatial interaction modeling tutorials in the R programming environment. Finally, some auxilliary functionality and future additions are discussed. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 2 The SpInt Framework"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Modeling framework"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The core purpose of the SpInt module is to provide the functionality to calibrate spatial interaction models. Since the \"family\" of spatial interaction models put forth by Wilson (Wilson, 1971) are perhaps the most popular, they were chosen as the starting point of the module. Consider the basic gravity model (Fotheringham & O'Kelly, 1989), \n",
    "\n",
    "$$T_{ij} = k\\frac{V_{i}^\\mu W_{j}^\\alpha}{d_{ij}^\\beta} \\quad(1)$$\n",
    "where \n",
    "\n",
    "$T_{ij}$ = an $n \\times m$ matrix of flows between $n$ origins (subscripted by $i$) to $m$ destinations (subscripted by $j$)\n",
    "\n",
    "$V$ = an $n \\times p$ and vector of $p$ origin attributes\n",
    "\n",
    "$W$ = an $m \\times p$ vector of $p$ destination attributes\n",
    "\n",
    "$d$ = an $n \\times m$ matrix of the costs to overcome the physical separation between  $i$ and $j$ (usually distance or time)\n",
    "\n",
    "$k$ = a scaling factor to be estimated\n",
    "\n",
    "$\\mu$ = a $p \\times 1$ vector of parameters representing the effect of $p$ origin attributes on flows\n",
    "\n",
    "$\\alpha$ = a $p \\times 1$ vector of parameters representing the effect of $p$ destination attributes on flows\n",
    "\n",
    "$\\beta$ = an exponential parameter representing the effect of movement costs on flows. \n",
    "\n",
    "When data for $T$, $V$,  $W$, and $d$ are available we can estimate the model parameters (also called calibration), which summarize the effect that each model component contributes towards explaining the system of known flows ($T$). In contrast, known parameters can be used to predict unknown flows when there are deviations in model components ($V$, $W$, and $d$) or the set of locations in the system are altered.\n",
    "\n",
    "Using an entropy-maximizing framework, Wilson derives a more informative and flexible \"family\" of four spatial interaction models (Wilson, 1971). This framework seeks to assign flows between a set of origins and destinations by finding the most probable configuration of flows out of all possible configurations, without making any additional assumptions. By using a common optimization problem and including information about the total inflows and outflows at each location (also called constraints), the following \"family\" of models can be obtained:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "&Unconstrained \\ (Gravity) \\\\\n",
    "&Tij = V_{i}^\\mu W_{j}^\\alpha  f(d_{ij}) \\quad & (2) \\\\\n",
    "\\\\\n",
    "&Production-constrained \\\\\n",
    "&T_{ij} = A_{i}O_{i}W_{j}^\\alpha f(d_{ij}) \\quad & (3) \\\\\n",
    "&A_{i} = \\sum_{j} W_{j}^\\alpha f(d_{ij}) \\quad & (3a) \\\\\n",
    "\\\\\n",
    "&Attraction-constrained \\\\\n",
    "&T_{ij} = B_{j}D_{j}V_{i}^\\mu f(d_{ij}) \\quad & (4) \\\\\n",
    "&B_{j} = \\sum_{i} V_{i}^\\mu f(d_{ij}) \\quad & (4a) \\\\\n",
    "\\\\\n",
    "&Doubly-constrained \\\\\n",
    "&T_{ij} = A_{i}B_{j}O_{i}D_{j}f(d_{ij}) \\quad & (5) \\\\\n",
    "&A_{i} = \\sum_{j} W_{j}^\\alpha B_{j} D_{j} f(d_{ij}) \\quad & (5a) \\\\\n",
    "&B_{j} = \\sum_{i} V_{i}^\\mu A_{i} O_{i} f(d_{ij}) \\quad & (5b)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "where \n",
    "\n",
    "$O_{i}$ = an $n \\times 1$ vector of the total number of flows emanating from origin $i$\n",
    "\n",
    "$D_{j}$ = an $m \\times 1$ vector of the total number of flows terminating at destination $j$\n",
    "\n",
    "$A_{i}$ = an $n \\times 1$ vector of thevorigin balancing factors that ensures the total out-flows are preserved in the predicted flows\n",
    "\n",
    "$B_{j}$ = an $m \\times 1$ vector of the destination balancing factors that ensures the total in-flows are preserved in the predicted flows\n",
    "\n",
    "$f(d_{ij})$ = a function of cost or distance, referred to as the distance-decay function. Most commonly this an exponential or power function given by,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "&Power\\\\\n",
    "&f(d_{ij}) = d_{ij}^\\beta \\quad & (6) \\\\\n",
    "\\\\\n",
    "&Exponential \\\\\n",
    "&f(d_{ij}) = exp(\\beta*d_{ij}) \\quad & (7) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $\\beta$ is expected to take a negative value. Different distance-decay functions assume different responses to higher costs associated with moving to more distant locations. Of note is that the unconstrained model with of power function distance-decay is equivalent to the basic gravity model in equation (2), except that the scaling factor, $k$, is not included. In fact, there is no scaling factor in any of the members of the family of maximum entropy models because there is a total trip constrained implied in their derivation and subsequently incorporated into their calibration (Fotheringham & O'Kelly , 1989). Another aside is this in the doubly-constrained maximum entoropy model the values for $A_{i}$ and $B_{j}$ are dependent upon each other and may need to be computed iteratively depending on calibration technqiue. It is also usually assumed that $n=m$ for doubly-constrained models.\n",
    "\n",
    "The family provides different model strucutre depedning on the available data or the research question at hand. The so-called unconstrained model does not conserve the total inflows or outflows during parameter estimation. The production-constrained and attraction-constrained models conserve either the number of total inflows or outflows, respectively, and are therefore useful for building models that allocate flows either to a set of origins or to a set of destinations. Finally, the doubly-constrained model conserves both the inflows and the outflows at each location during model calibration. The quantity of explanatory information provided by each model is given by the number of parameters it provides. As such, the unconstrained model provides the most information, followed by the two singly-constrained models, with the doubly-constrained model providing the least information. Conversely, the model's predictive power increases with higher quantities of built-in information (i.e. total in or out-flows) so that the doubly-constrained model usually provides the most accurate predictions, followed by the two singly-constrained models, and the unconstrained model supplying the weakest predictions (Fotheringham & O'Kelly, 1989). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Calibration framework"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spatial interaction models are often calibrated via linear programming, nonlinear optimization, or increasingly more often through linear regression. Given the flexibility and extendability of a regression framework it was chosen as the primary model calibration technqiue within the SpInt module. By taking the natural logarithm of both sides of a spatial interaction model, say the basic gravity model, is is possible to obtain the so-called log-linear or log-normal spatial interaction model:\n",
    "\n",
    "\n",
    "$$\\ln{T_{ij}} = k + \\mu \\ln{V_{i}} + \\alpha \\ln{W_{j}} - \\beta \\ln{d_{ij}} + \\epsilon \\quad (8)$$\n",
    "\n",
    "\n",
    "where $\\epsilon$ is a a normally distributed error term with a mean of 0. Constrained spatial interaction models can be achieved by including a fixed effects for the origins (production-constrained), a fixed effect for the destinations (attraction-constrained) or both (doubly-constrained). However, there are several limitations of the log-normal gravity model, which include:\n",
    "\n",
    "1. flows are often counts of people or objects and should be modeled as discrete entities; \n",
    "2. flows are often not normally distributed;\n",
    "3. Bias\n",
    "4. zero flows \n",
    "\n",
    "Therefore, Flowerdew (1984) proposes the Poisson log-linear regression specification for the family of spatial interaction models. This specification assumes that the number of flows between $i$ and $j$ is drawn from a Poisson distribution with mean, $\\lambda_{ij} = T_{ij}$, where $\\lambda_{ij}$ is assumed to be logarithmically linked to the linear combination of variables,\n",
    "\n",
    "$$\\ln{\\lambda_ij} = k + \\mu V_{i} + \\alpha W_{j} - \\beta d_{ij}) \\quad (9)$$\n",
    "\n",
    "and exponentiating both sides of the equation yields the following family of Poisson log-linear spatial interaction models:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "&Unconstrained \\\\\n",
    "&T_{ij} = \\exp(k + \\mu V_{i} + \\alpha W_{j} - \\beta d_{ij}) \\quad & (10) \\\\\n",
    "\\\\\n",
    "&Production-constrained \\\\\n",
    "&T_{ij} = \\exp(k + \\mu_{i} + \\alpha W_{j} - \\beta d_{ij}) \\quad & (11) \\\\\n",
    "\\\\\n",
    "&Attraction-constrained \\\\\n",
    "&T_{ij} = \\exp(k + \\mu V_{i} + \\alpha_{j} - \\beta d_{ij}) \\quad & (12) \\\\\n",
    "\\\\\n",
    "&Doubly-constrained \\\\\n",
    "&T_{ij} = \\exp(k + \\mu_{i} + \\alpha_{j} - \\beta d_{ij}) \\quad & (13) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $\\mu_{i}$ are origin fixed effects and $\\alpha_{i}$ are destination fixed effects that achieve the same constraints as the balancing factors in equations (2-5). Using Poisson regression automatically alleviates limiations (1) and (2) and since we no longer need to take the logarithm of $T_{ij}$, issue (4) is also absolved. Using fixed effects within Poisson regression to calibrate the doubly-constrained model also avoids the need for iterative computation of the balancing factors that exists in other calibration methods. \n",
    "\n",
    "Calibration of Poisson regression can be carried out within a generalized linear modeling framework (GLM) using iteratively weighted least sqaures (IWSL), which converges to the maximum likelihood estimates for the parameter estimates (Nelder & Wedderburn, 1972). To maintain computational efficiency with increasingly larger spatial interaction datasets, SpInt is built upon a custom GLM/IWLS routine that leverages sparse data structures for the production-constrained, attraction-constrained, and doubly-constrained models. As the number of locations in these models increases, so the does the number of binary indicator variables needed to construct the fixed effects that enforce the constraints. Therefore, larger spatial interaction datasets become increasingly sparse and the utilization of sparse data structures take advantage of this feature. As a metric, constrained models with $n = 3,000$ locations, which implies $n^2 = 9,000,000$ observations can still be calibrated within minutes on a standard macbook pro notebook.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 An Illustrative example: migration in Austria"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 The data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following example purposefully utilizes data that has was previously used to demonstrate spatial interaction modeling the R programming language for validation and consistency. The data is migration flows between Austrian NUTS level 2 municipalities in 2006. In order to use a regression-based calibration, the data has to be transformed from the matrices and vectors described in equations (1-5) to a table where each row represents a single origin-destination dyad, $(i,j)$ and any variables associated with either location. Details on how to do this are outlined further in (LeSage & Pace, 2008), though this has already been done in the example data. Let's have a look! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Origin</th>\n",
       "      <th>Destination</th>\n",
       "      <th>Data</th>\n",
       "      <th>Oi</th>\n",
       "      <th>Dj</th>\n",
       "      <th>Dij</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AT11</td>\n",
       "      <td>AT11</td>\n",
       "      <td>0</td>\n",
       "      <td>4016</td>\n",
       "      <td>5146</td>\n",
       "      <td>1.000000e-300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AT11</td>\n",
       "      <td>AT12</td>\n",
       "      <td>1131</td>\n",
       "      <td>4016</td>\n",
       "      <td>25741</td>\n",
       "      <td>1.030018e+02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AT11</td>\n",
       "      <td>AT13</td>\n",
       "      <td>1887</td>\n",
       "      <td>4016</td>\n",
       "      <td>26980</td>\n",
       "      <td>8.420467e+01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>AT11</td>\n",
       "      <td>AT21</td>\n",
       "      <td>69</td>\n",
       "      <td>4016</td>\n",
       "      <td>4117</td>\n",
       "      <td>2.208119e+02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>AT11</td>\n",
       "      <td>AT22</td>\n",
       "      <td>738</td>\n",
       "      <td>4016</td>\n",
       "      <td>8634</td>\n",
       "      <td>1.320075e+02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Origin Destination  Data    Oi     Dj            Dij\n",
       "0   AT11        AT11     0  4016   5146  1.000000e-300\n",
       "1   AT11        AT12  1131  4016  25741   1.030018e+02\n",
       "2   AT11        AT13  1887  4016  26980   8.420467e+01\n",
       "3   AT11        AT21    69  4016   4117   2.208119e+02\n",
       "4   AT11        AT22   738  4016   8634   1.320075e+02"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAACjCAYAAABmOs66AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnWV0VFcXhp8bdyNKEiBo0ODubsVdW+CjOLRQoIKWosWt\nWHGKFXeXUCyQAEECASIQd89kcs/3IymFEkoSYsB91mIxM/fIvqxhz7n77LNfSQiBgoKCgsLngUZ+\nG6CgoKCgkHcoTl9BQUHhM0Jx+goKCgqfEYrTV1BQUPiMUJy+goKCwmeE4vQVFBQUPiO0MttQkiQN\n4BbgL4ToIEmSC7AaMAR8gL5CiLgM+vkA0YAMpAghauaA3QoKCgoK2SArK/2xwP3X3q8DJgohXID9\nwMR39JOBxkKIKorDV1BQUMhfMuX0JUlyANoC61/7uLQQwjX99Rmg67u6Z3YeBQUFBYXcJbPOeDHw\nHfD68V1PSZI6pL/uATi8o68ATkuSdFOSpP9lz0wFBQUFhZzgvU5fkqR2QLAQwoO0VfvfDAZGSpJ0\nk7S4vuodQ9QTQlQl7UlhpCRJ9T/QZgUFBQWFbJKZjdx6QAdJktoC+oCxJElbhBADgFYAkiSVAtpl\n1FkIEZj+d6gkSfuBmoDrv9tJkqQUAVJQUFDIIkII6f2t/uG9K30hxA9CiCJCiOJAL+CcEGKAJElW\n8Cqr5yfgt3/3lSTJQJIko/TXhkBLwPM/5ipQf6ZNm5bvNig2fTo2FVS7FJs+Xpuyw4dssPaWJMkL\neAC8FEJsApAkyU6SpCPpbWwAV0mS3IFrwGEhxKkPmFNBQUFB4QPIdJ4+gBDiInAx/fUyYFkGbQKB\n9umvnwOVP9xMBQUFBYWcQEml/A8aN26c3ya8hWJT5iiINkHBtEuxKXMURJuyg5TduFBOI0mSKCi2\nKCgoKHwMSJKEyOmNXAUFBQWFTwfF6SsoKCh8RihOX0FBQeEzQnH6CgoKCp8RitNXUFBQ+IxQnL6C\ngoLCZ4Ti9BUUFBQ+IxSnr6CgoPAZkWmnL0mShiRJ7pIkHUp/7yJJ0l+SJN2RJOng34XVMujXWpKk\nR5IkPZYkaVJOGa6goKCgkHVyVS4xvQLnCtJKMJcnrUibc/bNVVBQUFD4EHJbLrEm8EQI4SuESAF2\nAh0/wF4FBYUcRpbl/DZBIQ/JbblEe8D/tfcv0j9TUFDIZ1JTU3G0t0dTU5PiRRzp1rkzV69e5fTp\n0wQFBb3RVghBUlJSPlmqkJO8t7Ty63KJkiQ1fu3SYGCZJElTgEO8Wy4x00yfPv3V68aNG38yVe0U\nFAoi169f50VAAGeG98c/KoaV190Z1KMbFgYG/PXYGz8/PxwdHQHYunUrAwcOZED//kQGB1GsREnG\njh9PiRIl8vkuPi8uXLjAhQsXPmiM91bZlCRpNtAPUJMulwjsE2lyiX+3KQVsFULU/lff2sB0IUTr\n9PeTASGEmJfBPEqVTQWFPGTMyBFYPnvI1FaN3vhclgWa42dSuXJlyhR34rGXF+73H+BsY0VXl7KU\ns7bkUVgEq67epl27dpQo44yjgwM9e/XCwMAgn+7m8yQ7VTazVFpZkqRGwHghRAdJkqxEmu6tBrAR\nOP+3etZr7TUBL6AZEAjcAHoLIR5mMLbi9BUU8pAfvv+ePb+v58n3I9+61mvHQR4Gh9GzkjMtSzlx\nPyiEvtUqoaX5T0Q4Ij6RLW538AyNwC8qhqexidz38kJPTy8vb+OzJq+d/hhgJGlx/n1CiB/S29gB\n64QQ7dPftwaWkrZ/sEEIMfcdYytOX0EhD7lx4wa1atVCLJ6WI+N1336AY55e1K9di+Zt2qKhqUnn\nzp0pXrx4joyv8Da57vRzE8XpKyjkLV8NHMi5o4fx/WlMjo0Zm5TMiUfeXPELICIhkZ237tGoXl1+\n37oNY2NjjI2NmTd7Nn4+PkiamlhbWfHjtGno6OjkmA2fE4rTV1BQyDSmJiZ0KVuCjb1zL4s6VZb5\n+Ywrv5y6hK6ONro6urgUtqFH+VIcvv+YYw+fANCwQQO2bttGkSJFcs2WTxHF6SsoKGSKqKgozM3N\n2dG/C72rVsz1+VJlGQkJv6ho7E2N0dbUBCAsLoFaS9bzLDwSAKdixXj2/PkHz5eUlERQUBDFihXD\nx8eHqKgoKleu/MHjFjSy4/Tfm7KpoKDwcRMZGcmDBw8oXLgwpqamWFhYsHfvXio62NG1Urk8sUFT\nI20DuJiF2RufWxoZ8PSnMajUqbRbt4Mzj5+xa+dOuvfogYZG9kqDybJM+bJleebjg76uLonJyVia\nmtJ/4EBatmmDjo4OlpaWVKpU6YPv62NEWekrKOQgkZGR9P+yP6oUFV07d2XIoCFopq9q8wtNTU1k\nWUZXW5vklBRaNGyASSFL/ty/n8WdWjGuUe33D5JHXPT2YcLxiwTGJ9C5cxe69OhBgwYN0NLK/Pp0\n3bp1zJv6Ex7jBqOWZRJUKcQlq1h88RoX/AIw1tPF0/8lU6dMZfKPP+bi3eQ+SnhHQSGf6dOvD96h\n3rT9ui0HFh+ABDi47yBFixbNF3uCg4OxtbUlcMZ4bE2MeBIazuOQcNxeBqGhITG2fk1M9HTzxbb/\n4mFwKPvvebHv0VN8wiMp5uBAXEICcfHxOBYujIODI1a2tgS+eIFzhQocPXaMGjVqkJKSwrHDhzgw\nsBv1i797f2Dd1VsM3X2EyRMnUrJUKfr07Yu+vv6r60IIJClLvjRfUJy+gkI+cuToEb5o/wVzT86l\nesvqCCHYt2QfGyZvoHPXzlSpXIWh/xuKubl5rtuSlJTEpo0bWbNyBSYpyZwc0gs97Y8zmvsyKoaX\n0bGoUlNxMDN59T4wJhY7E2P+8g/g1CNv4lQpVHWwY3brxpSztfrPMR+HhDPuyFlqFbZm1z0vEjQ0\nqVq1GqbGxqSmqvnz4CHq1qzJ9t27sba2zqM7zTqK01dQyEeePXv2qizBGXHm1edBPkFcP3qdNRPW\noEpS4ebmRtWqVXN1JdmtYwe877jzbb3q9KtWCQ2Ngr9qzS9SUlM59uAJbi8CMdLRRl9bG1tjI6ae\nvEjtlq3ZtGVLfpv4ThSnr6CQzzRs0pAqParQYXiHt66lpqYypvYYvNy8uHTpEg0aNMjS2EIIwsPD\nMTIy+s9Tr/fv36dK5cpE/zIRfR3tLN+DQhpf/nGQzTc8KMh+KTtOX1HOUlDIQQoVKkREYESG185s\nPYOXmxdOTk7UqFEjy2MfP34cKysripcszuHDh0lJSQHAzc2NAwcOvGrn5eWFvo624vA/kD5VKgCf\nXunpD1HOqixJ0tX0z25IklT9Hf180tW13CVJupFThiso5BWJiYmZahcfH8+zp8/Q0894Fd6wW0PK\n1ynP8+fPadGiBUeOHOHkyZPcvn2bp0+fkpCQkGE/lUqFj48PPj4+WDtY029WP/43/H8YGhpiZGRE\ni9Yt6DegH81bNcfd3Z24uDhK2RbcOPTHQmFTYwAWL1qUz5bkLFnZ2flbOcsk/f08YJoQ4pQkSW2A\nBUCTDPrJQGMhROQHWaqgkIcIIZBlmXXr1jF8+HCatWzGiqUriIyMpEqVKm+FV+Lj49m+fTsP7j/g\n10m/ZjimvpE+iy4vYs34Nexbug9XV1dKu5QmPiYeBESERFC7bm109XTp16sf7dq14/79+0ycOJG/\n/voLs0JmfP/H91RrUY1WX7YiPiaeG8dvUK9TPQKfBXJg+QGqVq0KQCFDAyLiEzDT18t2vvvnTgU7\nawbUqU5CJn/0PxYyFdNPV87aCPwCfJtecO048LsQYo8kSb2BdkKIfhn0fQ5UF0KEv2cOJaavUGDo\n1LUTB/cdBGDS5knERsSy6adNaOtqY25uTscOHZkxbQampqacOnWKVq1aATB191Qadm+YrTmTEpK4\neugqiXGJXNlzheunrgNQt11d2o9qT83WNf+zvzpFzbN7zziw7ACnNp969bm+jjYdypXm144tcTAz\n+Y8RFP7NTycucF2lwelz5/LblAzJtY1cSZL2kObwTfmnyqYzcBKQ0v/UFUL4Z9D3GRAFpAJrhRDr\n3jGH4vQVCgzde3Tngf8DZh2dhYlFmqNUJavQ0tbi4bWHHFl9BK8rXvTp3Yf9+/bj9ciLHX47sHbM\nubDKs3vPkCQJpwpOWe772O0xIf4h6Bno4X7uNnsX7iU1VcbRwowB1Sowq22zHLPzUyU8PoGWv23D\nNz6JsIiM92nym1xx+unKWW2EEKPSlbP+XukvJa2G/gFJkroBXwshWmTQ304IEShJkhVwGhj1mrbu\n6+0Up69QYNi2bRtzls5h2c1lGV4XQuDl5sWVP6/wxO0JY9eNxc7JLo+tzDxxUXFc3H2RJcOXIGTB\n77068FWtKvltVoElLC4BqykLgDSFsZo1//spK7/Irdo79YAOkiS1JV05S5KkrUB7IcRYACHEXkmS\nNmTUWQgRmP53qCRJ+0kTS3/L6YMil6hQMEhNTWXBogV0ntT5nW0kScK5hjPONZzz0LLsY2RmxJHf\nDiNkwYDqlehcsWx+m1RgCYqJo+GqLWhpaXHjxg2qVCk4P455Ipf4RuM3RVTuAyOEEBclSWoGzBVC\n1PhXewNAQwgRJ0mSIXAKmCGEOJXB2MpKX6FA4OPjQ/kK5fkz4k+0P6G0x60zt7J52mbali/N0SG9\n89ucAkuxWcswsbPn2vXrBV7+Ma/z9IcCCyVJcgdmpb9HkiQ7SZKOpLexAVzT21wDDmfk8BUUChJG\nRkZoaGgQEVQw47jZpf/U/iy8uJATD73R+HYmRzy98tukAom2pibt2rYt8A4/u2SpGIcQ4iJwMf31\nFeCt3Pz0cE779NfPgU+viLXCJ83ylctp0rsJNkVs8tuUHMeloQuF7CwIfRnGiD+P0b5Cmfw2qUAR\nFpeAd0gYRZ2yvnn+saAk8CoopBMZGcnEyRNZunQpbb9um9/m5BqtBrWmcMnCvIyOzW9TChyuz/0A\naNeuXT5bknsoTl9BAbh16xYuVV14EPaA1XdWU6pqqfw2Kdf4cuaX6BvqYW1siFr9aZUYiEtSERaX\n8cnmzNC0lBPFbawoUqQI69ete1Xq4lNCKbim8FmRmppKbGwssiwjyzI3btxgxeoVXL16lTG/jaFh\nt+wdrPrY6FesL0G+wVibGBE8Y3x+m5MjLLpwlfEH07YMS9ta8fC7YQBZPpEshGDDdXf+t+sw9+7d\no0KFCjlua06hVNlUUPgXISEh7N+/n4SEBO49uMfBAwdRqVRoamoiaUjYF7en9bDWNOnVBD2Dd1eu\n/BQ5s+0MSwcvJGHeD/ltSrbxDAxml/t9Zp2+DICeiTH9929nXbN/qpx2qFCGg4N7ZWlclTqV0vN/\nY8KUqfTq3RtLS8sctTunUJy+gkI6x44fY8r0KXg99KJeh3qYWJlg6WhJ/S71sS1mm9/mFQiWjVzG\noVWHEIun5bcpmSYqIYlivyzDysgQfU0N7gWGAGBb3hnTIo60WzgLm7Jl8Nx/mKPf/EBMUAjq5GRs\nTY0JnP5tlua6FxDM8EOnueL1lMTExP8sZ51fKMLoCp89siwzfeZ01q5fy6jVo6javCq6+gVPDvBD\nOLT6EBsmriM5SYW5pSnNvmzJkDlD/rOPLMv80HIyPp7PSVGpUatTiY9NwKXwx5Wh1Hz1FqITEolX\np6JWqShatxZDzx9GS0fnjXYVOn9Bhc5fALCtx5fc23MAWZazFOqpWNiGFV80p2lgKCqVqkA6/eyg\nOH2Fj5qUlBSmzZiG9zNvQkND8XrkReFShVnuthwLW4v8Ni/HkWWZFaOW82UNFzpXdObAPS82L9iN\nzz0fZh6amaFTUyWp2Lt4L25nbzOsbjUsDQ1ehUM29e6Y17fwQRQyTNOx7bl9LeU6tc+UYLqBmRnG\nerrZqjbqUtiG1s7FWfzrr0ybOTPL/QsiSnhH4aNFCMG3E77lgtsFWv+vNaZWptgWs8WxjGN+m5Zr\nRIZE0t2mOykLpqCllebEXJ/50fy3rega6tH5m66kqlM5/ftJtPW0CXsZRlKSCh0tTSY1rcvMNk0B\niElKosHyTezo14Xydh9P7X2TyXOITVYxJdQbo0zE2dVqNT9qp7V7NHkkZWyyHpt38wug2+6juN/z\nzBN946ygxPQVPhvOXzjP+O/GExgSyMLLCz/Jg1QZsWrcKg4uP4B64ZQ3Po9JSmLc/pPs9LhPoiqF\nRiWK4mBmQuMSxWhfvjS2Jkb5ZHHOovHNDASgqaPNZF9PtPTSQncvb7rj1KQBWlpaxIWFcXXlei7O\nWURKsuqN/tnZv1Cnygw/cJL9nl6M/eZbyjg706NHj5y4nQ9GcfoKnzxubm7UqFED80Lm9PyxJ+2/\nbv9ZZd1813QC7uc9ePzDKEpZFXrreloqKq+eAj41ZFnGaPIcElPUb3ye7vwwMDYiITYOgL7VK7Ks\nU2tW/+XG8ituJKWkEPXLpGzNe+zBE7Q0JA4+fMpvl6/z5759dOrU6YPv50PJVacvSZIGcAvwTy+4\nVhlYDegBKaQVX3PLoF9rYAlpB8E2CCHmvWN8xekrvJcBAwawdetWVtxY8dFUuMxJOpp2JD4mnmWd\nWzO6Ya38NidfeB4eyaWnvvStWomA2Fhmn77Eok6tCIyOY6eHJ/cCQ2hUvCjD62ddhzgjph4/z8+n\nLr31eVBQEDY2+fuEmdtO/xugGmCS7vRPAgtfk0ucKIRo8q8+GsBjoBkQANwEegkhHmUwvuL0Fd7L\nvXv3aNS0EePWjaNep3r5bU6eIssyLTVb0r5cKQ7/r09+m/PJsvD8X0w/eYlUIRjbsBb9qlag88bd\nPAkNB0nC1MqR6BA/EhIS0NfXz1dbc63KZrpcYltg/Wsfy6QpaQGYAS8z6FoTeCKE8BVCpAA7gY8r\nXUChQFGxYkUO/HmAzRM3M7vnbCJD/pFevnrkKnMHzMXjggdbZmxhUouJzO47Ox+tzVnU6rSQRhFz\n0/e0VMgO9wNDqPrrGiYcOo2tS1OM7Usz98xlqixci3dYOMZmlvx4IgFN7bT0UB8fn/w1OJtkNmVz\nMfAd/zh5gG+Ak5IkLSRdLjGDfvbA6xKKL0j7IVBQyDYNGzbE844nP037iWGVhjHj0Ay2z9zG1aPX\ncLQ15czWM+jpapGSkkqqLPB74MuKmyszTO+7euQqF3ZeYOxvYzEwKtildHV0dHBp5MLqi250qVSW\nZqWL57dJHzUnH3oz68xlIhOTCY6LJyx9L6BI+br0nXcMgFDfh+yY3JaoIB86fL8NLR09KjXrxYUt\ns6hRuy5aGhJ/7t1Ls2ZN8/NWskSuyiVKktQVaCWE+LvWfj+gphBiTAbziGnT/tlZV5SzFDLDzJ9n\nMm1q2vfm7KoBNK3u9MYhnI7j/+DQ5cev2tfpUIcWA1rw5697eH7vOQnxSQBoamnStHcT2g//gvJ1\nyuf9jWSBMbVH8+D6Q4pbWTCohgvfN6ufrRz0z50ai9bhERROUZdGWBevRNX2Q7B0zLjU9L8PdsWE\nvuTmvqW47lzA6t/WMuzr/+WJzf9WzpoxY0auaOTOBvoBatLlEoH9pMklmr/WLloIYfqvvrWB6UKI\n1unvJwMio81cJaavkB2io6OxtraiTkVHvh9Qh1Z1Sr7VRqVSc9btOe2/+QNZCLS1NSlia8bo7tUJ\nDo9n/WF3EpNSiEtQYWFtxu7gvflwJ5lHlmVun73N8fXHcd3nipYksaV3J3pUKdg/VgWFS94+DP/z\nOA+DQrAu4cKw9R5ZHuPwgiHcPraB58+fU6xYsZw3MpPkespmNuQSNQEv0jZyA4EbQG8hxMMMxlac\nvkK2uHfvHqdPnWL8hAmM7FGHFRNaZrqvWdN5oK9LfFQ86hQ1bQa3Yfz6j6fqpCzLLB6yiBObTiKE\nwERfj+IWZkQlJXPwq55Usv88zi+8zsPgUPps249nYAjq1FS0NDVxH/8/KtjZEBQTh920hegamFCr\n2zfU6f4tekYmWZ5jzeCKJIT5ER0dnQt3kHny2unXA5YCmkASaT8A7pIk2QHrhBDt0/u0Tm/3d8rm\n3HeMrTh9hQ/Cw8ODVi2asemn1rSpm7l6+MaN51KukQuzjswi8Hkgdk52H2WoRJZl4qLi2LdkH1cP\n/kXA0wBSklTc/W4YZW2s8tu8PEOWZTTH/4y2rgGtRy/FsogzZ9dM4sXDawyrU401V29RyKEkwzd/\nmFTkvp97oxX9lNtuN3LI8uyhHM5S+OxxdXWlW+cvuLN1MDaF/vsUqkqlRrf+L7Qf1p5xq8flkYV5\nR0fTDlhq6+Dz09j8NiXXCYiOYeml66x0dSNepWL0Nm8s7Eu8ur7126Y8cz9Pvd6TaD40w3Xne0mI\nDiMqyIeXj25yfOlIli9fwciRI3LqFrKF4vQVFIDvJ33Hw5snODCv63vbmjWdS3RcMmfEmTywLG84\ntfkk879cgCRJTG3ZkOmtG+e3SbnO4J2H+P26O9p6hrT7ZhUuLQe81SYmLAATy8LZnmP1QGdC/Lxo\n36EzTRs3ZNy4sUhSlvxtjpNrefoKCh8TcXFxBITGkplFRHH7tFyEg6sO8uzus9w2LU+4tDetgqZq\n/k856vC3u92jxeotOMxYTJk5K5l6/Hy+yy2mpqbNP7ZB2unklKT4DB0+8EEOH6DbzP2YWjtw985t\nRo4cke8OP7sopZUVPhlSU1NZsngxJ48d4NbGLzP1n3LLtE50mbybVWNWIssy5lZmaOtokxSfSHxc\nIjZFbNDR16HlwJZ0n1Awimy9Dwu7tJLSJ7ye0L58ximImWXy4TMce/iEao6F2XLzDqYGBtQuVYqI\nuDh+PnWJK8/9mNe+BdWLfJhDzQ6zTl1iyvHzr95r6+rz1TLXXJvPqmhZev58kLVfV2Pr1m0MHjwo\n1+bKTRSnr/DJMKB/P3b8sZNHe0ZibJg54ZQKJW14vHc0AEFhcYxbfAIjfR3cHgZwJzKOl08DAFg3\nef1H4fRVKhVeN9M2KR1Ms56VAhASG8eLqFgiExKZd+4K1ZycOPXEl8IWFqwcNAgTg7RDbDP27OGy\nlxe1l6ynQmEb7IwNWdalTYaF4HIDHU3NV6+LVmpAnznH0DF49z6OKjGeqCAfdPSNMLMtmq05NTS1\nMDA0olOnj7ewgBLTV8h1kpOTWbVqFb1798bWNnekCsPDw7G0tKR9w7Ic/vXDnbNF8wVU61iXTqM7\ncf6P89TpUIfKjSvngKU5jyzLqFVqrh29xsrRK4kOiUKdmgrA+REDWXXlJld9X9CslBPzv2iOtbER\nIbFxjNl3Ap+IKB6GhKOjpYGNkSGBMXHEJqtISe8PcHjSJIz+QzVq4/nzPA0O5opX2o+Nka4uqbKM\npobExZFfUtXRLtfuvc6SDTxO1mL0H8/f2SbipTdrBldCQ1MDIyNjwkKCqNb2KyydKmFqXYSStdqg\nrZv5Gjp/fNuQ1b9Op2nT/D+Fq2zkKhQohBBcvHiRqVOnc/nyRVauXMmIEbmT7TBm9EjWrVvLwrEt\nqVDcmvVH7qFWy4zrWY2a5e2zPF7twb9zw9MfbR1tDE0McCjjyMjloyhZ+e3DXwDu59zZMHk9eoZ6\nBD0PItgvBH1DPdoNa08hu0LUaFMDOyc7dPR0MuyfVcICwpjd6xeCngURHR5NclJa3fiStrYs7N+f\nK15erDhxggSVCktjY6ISEtDU0AAh+KqmC7/f8EBXW4fiVlZUdnLiZUQEd319AWjh4sKAhg1fzaWT\nCXUqAA8fH+76+pIqy9hbWLDvxg38w0JpWaYEVext0dbQ4PsWDXLk/v/GYcYSdIpXZ8DCjDfiY8MC\nOLZ0JI9cDxAZGYmZmRm+vr7s3r2H575+bNq4gWbDFlPti6FA2nf23pkdRAd602DA27X31apkdkxs\ngYOpBhfOn8v39F7F6SsUGO7cuUPlypUpXLg4+voWREQ859tvv+H69ZucPn2CTp064+DggKfnQ1JT\n1ezbtwdjY+O3xklJSSEyMhJPT0+srKwwNjamaNGib8XrAwIC+HnGVMLDQnnxwp+q1WpSxrkcM2dM\nZd7Ixgz6Iuur9IQkFct23eD3Q+488Y/ghx0/0LT326u7rTO3snnaZgx1dXGytsbM0JDi1tb4h4fz\n15PHaGpokJCUDIC5lRmLLi/+IHWvq4eusmLMCkJ8g+lQowZGenr0q1+f6MREbEzfLMamUqtfOe3n\nISEsPHyYsNhYyjk4MLVbt2zbkBkSVCp+PXSI697eJCSn3X/qwik56ignHT7N/HN/0eXHbdiXrf0q\nTfPArJ68fHKH+PCX1Kpdl769ejJkyNsxeFdXVzp06UHrb9YC4H//Kn/9MRdZlqnfZzKpyfEkhPkT\nE+pPRLA/8TGR6BsYokpKJCQkBBOT7IXQcgrF6SsUGGJjYylevAQNGvyPMmXqcujQL+jpGePkVIuE\nhGhiYgLR0THC0rIo+/bNYtGiRYwZM/pV/5iYGBrWr8Odew8wNjKgallHXgRHER2XSFhEDM2aNubM\n2fPvNiCd48ePM2b4IB7vHpqtbIshvxxiw0F3arepyaxjb1fs/L7FJG6euUU5BwdWDh78znFCY2JQ\np6Yyc9+fPA4IZHfwbswszV5dVyWpuHvxLiWrlXzjc0gL35zbcY4Xj19wbM1RIkIi0dXWZmizZnSp\nVfBr6qvUamRZpsOC+Sxo35yxjWrn2NiyLFN+/m88Cg5FS1uHH08lc3bNRFx3LsDUzIyTJ05Q6z3/\nRr8uWsKRYyfQ0tJCR1ub2T9Pw83NjeDgYAwMDHB0dHz1x8bGBs3X9hLyG8XpK+Q5QghiYmIwNX27\n3K+HhweNGzdjxIht6Oi8O2YaGurL9u3jad++La1bt8Td3YO/XC9RvjAs/aY5iclqLM3SNg8jYxKx\naD4fJycn3N3dM5z3dZKTk6lUoSxz/1eDzk3KZunewqISsGq5AFsnW7Z4b3lrhRrkG0S/Yv2Y0rUr\nTStUyNSYoTEx9Fi8mEOxh3A/68657WfxvuVNsF8wanUqeoZ6HIk78qq922k3lo9Yzkvvl+hoa1HM\n2pqfOnfBsVDebJbmJJ1//ZWo+HisjY0w1NHG0cyEk0P7oafz4fkkfpHRFJ25hC4/7cDQzIprO+cQ\n8eIJ+jrPNZoTAAAgAElEQVSa6OrpM37cGL78ciBCCDw9PSlTpsx7vzsfA9lx+p9t9o5KpUIIga5u\n5rI8FODhw4dcunSJW7fcuXXLnejoaGJjY4iICKN48ZKoVCpq166Nrq4ulpYWnDlzlujoCJ4+daNs\n2XfHcq2sijJs2CZOnVqKm9tSLCyKc+3GTb4Y3pSA0FhKFfnHwZmb6OP5x3DafbsTMzMznjx5QsmS\nGcfZAS5fukRoWDjaWllfnVmY6FGqaCGePA/i8G+H6TjizYyNuxfvoqOtlWmHDxAdHw9AB+MOaGhI\nWJmaUsnBkS7NK2FhbMzMvXvpZNoRAyMDhCaEvQijiJUlY9u0oVPNj7sq+crBg5m8fTs2pqbo6+jg\n7uOD/qRf2NS7IwNrZn+TXJZlGq3YlPZGyBSv1pzi1ZqTmqLC754ryQkxrPpjI+MnTEAIgaV9caKC\n/bh54zrOzp+f+tqHyCXuBEqnXzYHIoUQVTPo5wNEkya6kiKEyPCbm5cr/VWrfmPUqBFoamrx3Xff\nM3v2jByf48SJE/z22zouXjyPr6/Pq9ifEIKkpKR8V9x5HSEE8+cvYPLkSdSoURstLW1CQ0NRqZJR\nqZIxMjKhbFlnrlxxpUSJWpiY2OPoWB49PSM0NDTQ0dEnOjoESFu1p6QkoVIlEhj4CH19I5o2/Roj\nI4tM2yPLajZtHE1omC9JSYloaWpwekV/GlcrBqSt9ot3WUnvPv1ZuGjRO/8tL128SKPGafH8DVOy\nnmInyzI9f/iTo9ef8WfEPnR03tyE3TJjC1umb6FRuXKMbdMGc6PMiY8/DwlBQ5KwMzd/a5PUNzQU\n10ePCIiM5EVEOBO+6PBRruozw9Hbt/n18GHW9GjPgGou2V7xL7l4jW8OnKTtuJXU6PjuRIHk+Bg0\ntXXQ0tFjRhOJCRO+Y8GC+dk1v0CQp3KJ/7r2KxAlhJiVQb9nQDUhROS/r/2rXZ44/Xv37lGpUiWG\nDLmOn58rly79RFhYKO7u7lSqVAkTExNCQkKwtrb+oHm++24iv/66AEjbjNTS0iIpKYkmTZpx7dpf\nVK9eC3NzC6ytC7Fp08YMBT7yijVr1jJt2mwqV26PmZkt2tp6GBlZoKWljaamNklJcYSG+qKpqUW5\nco3y1LY9e6bz4MFFNk/rxIB2LgA88QundLcVbN+2jZKlSlG9evUMNwenTpnCz7NmEX1uMiZGWX+i\nm7vZle9XnmXYomF0++btTU+f+z5832oyMRGx6GlocWDChKzf4GfMhQcPmLFnDwAmenpEz8meaLn+\nxNk4VG5K3/knMt1n6/imPLt9Hm9vb0qUKPH+DgWUXHP66XKJG4FfSBdR+dd1P6CJEOJpBn2fA9WF\nEOHvmSNXnL4QAj8/P27fvs3Nm7dZvnwplSuPolmz2Tx4sJc9e7q/amtpaUtYWBAAwcHBH+T4Hz58\nSKVKLtjaFiU09AVmZhbo6+sRE5PA11//jp/fPVSqRM6fX0d8fCS2tvb873+DsLQsRFhYGCdOnMbX\n15evvhrI5MmTciwMFR0d/UYsUwhB+/Yd0dIqS5UqbXJkjpxkzpy2qFSJxF74HiODf1baF2758POm\n65y/4cXkSRMZP+E7Cv1rRWxf2I4apc04sKBntuYeMusQO84+5PBrMfaM+HvFf/LHHzOd3qiQhkqt\nJjg6mgErVtCwRFF6Vi6PkY4OVR1tqWCXVhZ6zunLlLIqRLfK5d7qL8sy2hNmYW5TlBHbn2Y6Mygh\nOpwDs/vx5MaJV4uyj5HcdPp7SHP4pqSXVn7tWgPSBNLfFbZ5BkQBqcBaIcS6d7TLMacvhGDcuG+5\nccOdBw/uoaWlR+HCVSlUqAq2tlUpXrw5Ojr/PIoHBaWJKCQlRaGtbcCuXV8gywlERISjqamZrd36\nb78dz9mzHnTq9BMpKUkkJESTmBiLmZkN+vpvpnklJsbw4sVDjh37FXNzW/T0zHFwqMDdu8cJDvbh\nwIEDdOz4YScAo6OjGTZsBDt37qBUqbKYmpoQEBBAbGwMhoZm9OgxB3Pz3DtE8yHMmtWCUo5mXF0/\nGDOTNw8JeT4NYczCk/iEpbB5y3bq1av36j9+v379SAq6y945XbI1r9uDAGp+tQ5rByt6/dCHdkPb\nZehUPC54MKHJBLaMHImjpWW25vrcOXr7Nvtv3OBlRAQCSE5JoYSlBaFx8cSkp7uKxW/nzQP8cOQs\nc866Mn5fMEbm71+oJcVFs3NyK4KfezJp4iSmT5uSk7eSp+SK089ALnG8EOKL166vIk38fPE7+tsJ\nIQIlSbICTgOjhBBvFcjIabnEUqWc8fb2olu3XZQvn7UTmklJUcybZ46engGpqWo8PNwpV+7tVcZ/\nsWHD74wcOZLWrUdRuXK7LPUFePTIlV27/vkyxsfHY2CQPQ3XyMhIrKysKFmyOk2bfo0spxIfH4WJ\niSWGhuYYGBTcLAYfnzu4um7n6dObOBez4uHut2O2siz4/bA7i3fdJiQygbDwSFq1bMnly5eYMbQx\nE/pmP60xLCqBFqO38eB5KLIsMDAxoHLTyjTo1hArBytunb7F1hlbqVGyJPP79v2QW1V4jbCYGAat\nXk1sUhL1nZ1xffSIexOHvVr9/83fefrFXBoxcMmFTI2tSoxjTltjli9fzqhRo3LB+twjP+US9wkh\nBqQrY70EqgohAt47mSRNA2KFEIsyuJaj4Z3k5GT69BnIvn270NMzon37tZQv3zvT/VNSEpBlNbdu\nLeP06SkMHDiETZsyfEh5J2vXrmXFiu106ZL1jWJZTiUgwAtdXUM2bBhOixYt2Lp1M2ZmZu/so1ar\nAV49qsbGxjJ16jQSEuJZuzbt8MmUKWfQ0Cg4ecbv45dfWqJWp1Dc3oI5I5rSo8W7JQGFENx5Eszv\nh9y5fDeA5d+2oH7lIjlihyzLtBy9jbM3n2Nrbk5EXCzqVBlZltHV1ubEDz/kyDwK/yDLMmpZJkGl\novOCBVgZGRDy83dAmuRhty37CI2NxdSyMOP2vMzUmKnqFDyOb+TIoq9xdXWlXr16uXkLuU6eKmel\nv28NTBJCNHlHewNAQwgRJ0mSIXAKmCGEOJVB21yJ6aekpKCjo0O/ficpUSLzMnqvc+PGSo4fH0VA\nQAB2dpkLgcTExNC9e09Onz7J1KnnsjXv36jVKrZuHY+BgczDh/ffuObt7c3x48dxdb3K0aOHiY+P\nw8WlCp07d2b69Kk4O9fG2roUenpGRES8pFWrkWhp5UwpgNzmypWdnDmzhjkjmzF5YP18tSUmLolC\nLRbQvmo1xrZtm6+2fE4kqVScuHMHt6dPueLlhb62Nno6OkQlJFCqZmsaD5qFXem3kgaJDQsgNjwQ\nq6JliQl9ga6hCadWj+fxlYOUcXZm7i8/06ZNwdvDyir5kaffE/jjX0a8LpdoA+yXJEmkz7U9I4ef\nm2hra9O8eVu8vPbg5NQsy6tctTqZ48dH0axZc2xtbYmNjWXdunXUrVsXXV1dqlSp8lYfIQRly5Yj\nIOAl48f/+cH3oKWlg7V1EdzcjtGzZx+EEAgB2tpa/PHHNqpUaYmDQyW6d/+FTZvGEhISzPz58+nd\nezalS9f54Pnzitu3j3L37iliY4KIT4hFnZJM3zaV8t3hA0TFJaNOlXHJRxHsz43nISEMWr361XsL\nOyfMCpekqEsjilVuRJGK7/5ebP22EaH+3q/eGxqbEB8bw9q1axk8eHC+18zJTz7pE7lCCDZu3Mzx\n4yc4duwEY8b4o6v7dn2X93HlylxevNiGqWkhwsICePrUG0nSQAgZlUqFtrb2q7bPnz+nQYNGvHzp\nzxdfTKBq1azH8zMiISEaHx8PkpPTDvdoaGghhIyxcSFKlEjTo09NVePpeQ5NTS3s7EpTqJBDjsyd\nF5w6tZqrV3dT0tGC8k5WpKSmMm9kcyqULDjC3pV6r0adpM2qIUPy25RPnpazZpGSmoqBsTlGlva0\nH78Wx/KZW8CkqlNYN7QqDWuUw8LcgmnTpmb6Cf1jQzmR+y8WLVrCkiWbqFBhMAMH/pQthw9Qrdow\nvL0PU7bs/7CxSaF9+ybs2tWWXr3avOHwHz16RMeOnSlatDaDBm3M0di5gYHpe3PkNTW1cHHJXggr\nL1CrVdy9expHxwpYWf1Tzzw83P/V601TO1LPJWfi8DnN+L51GPTzISLj4jJ9EEshe1iamBAYGcng\n1Tff0LrNDMeWjCT4uScVv+rBlCkfb2ZObvHJrvSvX79O69bt6dv3MpaWOXvU+uXLG6xfX4vAwEBs\nbW1JSUnhypUrdO/eg7JlW9KkyeCPVkrtQ5FlmUePXHF2rk94uD/x8VFcvryNwoVLc8fjOLFxaWf0\njAxNMLdwJDT0GUlJiQBYWxjxZO9ITIzeXbs9v5FqzmBAw4Z81STDbSyFHEKlVtN2zhwKFSnL8I2e\nWeorp6Zy/8JuTi8fie/zp5ibm+eSlfmPstJPZ/v2HYwePZ7mzVfmuMOPiwvi4MFeFC5chLt377J0\n6TLWrVuPvr4pZcs2o1GjgZ+tw09KimPt2sFERoagoaGBLKfplxob6BLw8g7F7c05v3c8sgy/bLzE\nvachtKrqzPLxbZGRMcihWvO5iSTB5UePqF6iBCb6+tiZm3PCw4OStrb8ccUV10de1Hd25uee2TsQ\nppCGjpYWFR0d8fC5jyzLWYrBa2hqUrFZbzxPbWLZ8hVMm6qs9l/nk1vpb9u2jf79+9Oq1RJq1x6b\nA5a9yeXLszl37sdX72vW7ICzcyOcnN7OIPjcOHBgLnfunGTdD18QFp3AiK7VC/SqPTvM2+zK3C1/\nEROXhJz+fdXU0CBVflMgvF+DBgwuAMpKHzNxSUl0WrAA50bd6TZ1Z6b7qRLjuHdmB67bZ9G9UzvW\n/rb6/Z0+UpTSyqRt3i5fvoKxY8cwcWI4+vqZL/SVGdTqJNavr0NwsAcVKzanS5cf39/pMyEpKY41\na4aiLcUQcvLTr0PzxC+cH1ad48dB9dHV1sLJzpxktZrWY7ZxzfMlG4YNo7hNwdmI/hjptWQJEUkq\nfjiRkOk+j1wPsGtKZ0aP/YafZ0z7JEoov4vsOP1PLm9JkiQaNkwr43v79socH19LS48+fY4C4Oz8\ncR/syGlkWU1UVCBLv22V36bkCaWKFGLP3O5ULm1HWScr9PS0MDXSY2B6YbiVJ08Sk5B5Z6XwNoXN\nzUlJTsxSn1sHV+JY1IllSxZ90g4/u3xyTh+gcuXK1K1bH2/vk7kyvptb2uNisWIFUyg7v3j27DYA\nvVtVzGdL8pdhXWtwfvVAbj9/TscFC3jw4kV+m/TRMjG95tTNA6sy3efFwxv4+z5nx44/3t/4M+ST\ndPqQVobBx+cKBw58SWJiRI6OXb36cABu3TqWo+N+7JQunSaD5+ru9962N++/pOOEnfT4fg+9f9xL\nrx/3Yt5sAQ2GbiQhXeT7Y6ZxtWIcmN8TYwNdRv2+gf4rVuDh45PfZn102JqZUd7RkWu7F2S6z7jd\nLyjfqCujx4xm7dq1HD58OBct/Pj4JLN3ANzcbuDt7U2pUqUwMrKiefPMf2n+ixcvrrNtW1r44r8k\nAD9HdHQMcLAvRfPRW0m49P07My5kWab+0E0YG1uiVquQNDSRU1Owsa3EzQcemDebT4papqSDBWN6\n1sTe2pjVe93wDY4BBF57Po4iWTo6mliZGxKfqOJFeDjfbN7M+WkZV4pUeDflHRx46nEv0+11DYwp\nXqUpCWH+fP311wDvVVj7nMgL5azWwBLSnio2CCHmvWP8XKm9U6VKdRISbOjd++gHj/XkyTF27GiH\ngYEZ9er1om5dJS3v36jVKubOacvQzlUY1b0mL0JiqFvJESMDHTYcvM3CHdd57BuWplo28fBbdYBS\nUpK5fn0fpqbW3Lp1CF/fu6+u1apTlOtXfSlTzJLihc1QpaQyf3QLbAsZYlvIKMtH62/ef8nqfW5U\nKW3LnSfBeD4LoVZ5e5aO//CaLHM2XuaH1eeo6uSES9GiePr707xiRVq6uHzw2J8TkXFxdFm4kNrd\nxtFqZIaFfDMkOSGWue3SSpgvXryYESNGvKV89imQ23n6Y4H7gAmAEKLXaxP/SlrN/H8bpAGsAJoB\nAcBNSZIOCiEeZcXI7PLixQs8PG5RqVIPgoLuUKhQKbS1s1eeGODQoa+wtS3J119nrdrm50BiYiyR\nkQHcuLGLVDmVNftvccDVhxLFnXAdsw1dbU1UahlLS3s6dvqB8uUboaHx9tdPW1uX+vXTqqFWrNgM\nf39PDh1aQFiYH1Y2xvwwowV33F/idT8UIQTVBqx91VdDQ2Jo52qsnvT+0hdP/MKp+dV67B3M2HX2\nPkZGejiXt2bZrhv8ef4RwRFxSEhYmOpT0tGCjg1KM7xrDXS0NAiKiKeI7X9vENpbmyBJUNbengGN\nMj5JLcsy/uHhFDI2Zuafe1GnprJowMD32v45cfjWLQCaDP4lS/1SU9JChPUbNKBly5afpMPPLply\n+unKWW1JV87KoEkPIKMjijVJq7Xvmz7OTqAjkCdOPzg4GAAvr8OEhV0nIMCXChV6ULx4K/T1Ld4S\nU3kflpbl8PG5gEqVgI5O9n88PnbUahX+/p68fPmIwEBPNDW18fW9h46ONj169GD69DE4OTlRvnx5\n3NzcqF+/AZrahnzV/xccHd9dGjkjHB0rMHLkZvz9Pdm4cTQR4fE8exrO4XNDKVXGitCQODS1NIiN\nTeLo/vv88M0RgsJiiU9KYdAXVejVMk20XJZlNh+9w3XPFzz0DeeJf9o+z6mrw7F3+Kdc9aljD1m3\n4io/9qyEubkBVy4/x+0vX2b8fomJy8+8amdnZUzz6sXo3KQsnRuXfcvuCiWs0ZCkV7n8/yZBpeKb\nTZt4HBgIgKWlEZGRCWw8f57udepgpPfh5xvu+Phw+dEjQmNjuOn9lL7169O3wbsF6gsiLV1c2HTh\nAsv7ODF+X3Cm+gghePHgGjYlKuF6+TK9+/TjjsftXLb04yFXlbMkSeoKtBJCDE1/3w+oKYQYk0Hb\nXAnvhIamrQitra25desWR48e4+pVN8LCwvD2fkqFCgPw9j6Ks3MXGjX6+Z3jxMWFsHChDbq6Rkyc\nePCjrtIXEvKclJQkbG1LoqmpTWpqCgEBXqjVKTg5vV019G+SkxPw9DyLu/shTEz0KVzYlh49ulGo\nUCFcXFxexUzv37/P6NHjSEhI4MGD+9St24fatbMmZJMRf5daBth99Etatn3b2U6ffIyt691IVqmJ\ni03GxsKI0Kh4ZFmgqSFR2MEMe0dTSjpbUatOEfp+lbG+bkY8fxqOnr422joadGq2Hs+7gZgZ6xF5\n9m1tV8OGcyhn78CCfv3f+NwrIIDgqChmHziArr4WI76pR+eelShT1oY+HTdz7NADgFc/GNpamhjp\n6bH8q0HYW1iw4sQJDt26RWpqKt1q10ZDkihsYYGpvj7HPTyISUwkLjkJv9AwAKysjTHQ1yYuPpnw\nsHgqFHFkXt9+GHxEK9/1586x/fJl2o5bRY2Ow9+45n//KiaW9ugamhAT+gKrYuWJCfFnSa+02k6N\nGjfm0MGDmJiYZDT0R0+BU84qCE7/v3j48CGbNm3l3LnTuLm5MXbsc8zMimXYVpZlVq8uS1jYYwYM\nWJjhCVx392P4+3sQHR1I06YjsLd/2ynlBbGx4fj63iEqKoC4uAgMDc2pVasbOjr63L59iMOHF1Oq\nlDNxccmUKdOIW7cO4OjoiJ+fD4ULl0ZTUwtjYyuaNx+Fv78njx9fIjVVzcOHrjRp0piBA/vTtWvX\nd5abaNfuC27cuEvz5kOxsyuNiYlVjtzXvXtn2LfvF7bs7UeHru9PC/3L9Tkrf72Eta0xXw2rRcnS\nVhgY5JyzW7HoMlO/O0pJx0KYGOpwdkV/TIz0mLbmPDM3XOLo99+/4Vz3XrvGypMn0dCQaNikJPtO\nDXrjB0eWZeLiVERHJRL4Mhora2NuXPXl6/67sDAyIiIuDoAiRc1JTRXExSajq6tFSEgsAOUq2GFl\nY4iVtRHx8cksWdMFG9t/nN3p416M/HIPYWFxWJuaUrdUaUZ/BDXlZVmm2c9pC7JhG+5hUzzt6S0l\nKYHZbQzfaCtJEgZGJsTHRr/6bPHiJYwbl/On8wsCuRXTrwd0kCSpLenKWZIkbXlNOasL8K4aBC+B\n10smOqR/liHTp09/9fpD5RIzQ9myZZk3bzb799ekT58B6Oq+O06roaGBtnbaF8zY+G0n9vjxVW7d\n2sWPP05GT0+PQYMG0afPXEqWrEHa1kbuERj4hCtXthATE4osq4mICKRJk2YUK2aLiUkxnjx5ypo1\nX2JhYUdcXChXrlyhbt267Nixg7/+usbMmYdo0KABSUlJnDp1Ci0tLX77bR3r1g0mOTmeqVN/QldX\nlwYNZuOSiY3I+vXrcuzYEUqUqJGjgi1FilREkiRcLz7LlNOvW9+JuvWdcmz+fzPq2wbY2Bmxf9dd\nzhz3ot7QjYzqVoOZGy7Ru169Vw7/6O3b7LtxnWfBIVSobIer+7gMx9PQ0MDERA8TEz0ci6QVCXMq\nUYht69144BnE8CH1KFnaisHD3ywxHBWVQMCLaMpV+O/ywS3alOFR4A+cPPKIIwfu88fmm1ibmtKz\nbt0c+NfIPX7cuRMNTU1K1GiNsWVhAGLDAzk4dwB29o543HYjMjISR0dHVCoVL1++5MWLF5ibm2Nn\nZ4e9vX0+30HO8W+5xOyQ28pZmoAXaRu5gcANoLcQ4mEGbfN8pf83s2bNYsqUKYwc+fCdBdrOnv2R\nK1fm0LXrFMqXT7vdwMAnpKamoFar2Lz5G/bt20fnzp2RZZnhw0dy9uw5EhJUtG37HQ4OWdPYzQwq\nVSLPnt3iypU/iIl5wYYNG3BycsLZ2RlDw39WQLIsc+vWLaKjo6lXrx76+u9PNU1NTeX+/fvY29tT\nqFChLNn18uVLmjdvgUqlTf/+S7N8X//FhQubuHhxM5bWRrh5jcfMrGDsrSxfeJEpE9LObXStVYtR\nrVtz8OZN9t64zouwcGxsjRn9XSOGjqqDjk7ByJSe8t1Rflt8hcUDB1Le0TG/zcmQxUeOcOjWLfrM\nOUqp2m2Jiwji+K+DeHD1OJUquTBv3lxat26d32bmG/khl7gRuCqEWPtam9eVs/7+YVjKPymbc98x\ndr45fW9vbwYN+h+PHj3G2NgRXV0jChdugbNzRywtnVGp4pgzx5iqVdvxxRdpNWXSavz0Ji4ukpQU\nFZaWVgQFBaKp+U8NfSEEVapUo2jRFlSpkvOP0du2fUtEhC9t27bj559n4OSUe6varOLi4sLdu3dz\nRb3r6VM3du78HrVaTcMmpTh0rmCImvRov5HTxx6hr6NDcoqaVFmmanUHOnavyP9G1c3R0FJOIMsy\nXVr+zsVz3qz/+usCVSdILcuM37KFu75+lKjenH4LThEd4s+p5aMoY6NDvz696dChwys96M8VpeDa\nB6JSqbh48SLPnj1jx459eHjcpFy5njRvvpRly0qgoaGmb9/5WFo68uLFAzZsGElKSgpCiDfEVP7m\n2rVrNGzYiO++O4S2tu5b14OCvLG0LJKtEIhKlcicOW25ePEiDRs2zNb95iYnTpzg3LlzLFiwgBIl\nqiNJKdjbV6Jx40E5Nsfx48u5cWMfURkf/cgX7tx+ybpVVylT1opO3Su9CtMUVO56vKRhlWWsGDSo\nwKz21bJM32XLCYuLo8uP2ynXuDu+dy5yZG5/vmjXmiWLflVq6qSjOP0cJjo6mhIlSpGUpKZSpf5c\nvboMc/PC2NkVxc/vPnp6eoSGvjuN7MSJE7Rp04Y+feagp2eMEDKRkS958uQywcHPSUqKIy4ulv79\nf6V48WpZsu3UqeUUKqTmwIF9bzxdFDTOnj3LunUb6Ny5I7169aJVq5HUrt0tR8ZWqRKYM6cdg0fU\nZtqcNpiYfFplnPOC//X9g6N/PuDo5O/z25RXnL9/n5///JNRW5/w4sFVLv3+I9qagp+nT2PIkMH5\nbV6BQhFRyWFMTU25f/8ely5dYsyYb3FyKo2BgR4DBnRm0KBDGBj8dzy5Tp06DBkyhNu39+Hv70to\naAiVK1dj2rRJVKtWjWLFijF27LesWjWRYsXKUadOf0qWfCvzNUOuXt3H/PnzC7TDB2jWrBnNmjUD\nIDExka+++oqqVdvlSAkLHR0DGjYcwI6Nf7J7mwe/7+xDizZlPnjczwVZlnnhF0VicsGqdTRz714A\nlvcriZ29Iwtmz6Rv3z7KAascQlnp5yHJycno6r4d5jl27Bjt2rWjcuWWNGgwAAuL/8422L37Bx4+\nvMquXbvp0aN7bpmb4zRv3gI/vxj69MnZcIxarWLjxpEEBj6lkKURqzd3V5x/Jvh9zTW+G3GQb9q1\no13VgiMC1GTGDEaMGMHQoUMpV65chqFThTSUlX4BJyOHD3D9+k0AHjy4yIMHFylbtj46OvoUK1ad\n2NhwYmOD0NIywMqqGE+eXMPf35Po6OiP7sCJr68f3t6P2bHjGxwcKtOgwYAckZaMigpGT8+E+vX7\ncPnydlYuuqw4/Uywa0vaKdWLDx5QuVgx7C1yVnAou5QpWpRevXplKj1YIesoK/0CQmJiIjo6Onh5\neXH8+HF8fX05dOgoVapUpl69OoSHR+DhcZe4uFimTZtC8+bN89vkLBMXF0d0dDReXl4MGjSEpk2/\noUiRzNfeDwjw4tatI/j7uxMTE4YsBFqa2iQmxmPvaEZYaDwlS1tx/ubIApMWWZCJiEhgzTJX5s04\nS62SJZnbt29+m4QQgk4LF3L7zh1KlCiR3+YUeJSNXIWPhgULFrBx40E6d/4Jbe3MbcD+8ktLQKbN\nF2VxrmCDn08Exsa6TJ/XBqNPTIs3L2nTYDVXXX1yreyzEOKNJzpZCNyePuWevz+psoypvj4VixSh\nrL09kiQxctMm6rdsyZy5c7EoIE8fBRUlvKPw0fD1119z7dpNVq8eQMWKLXB0dKFoURe0tHQyDPls\n2gxEnEYAAAyTSURBVDQatTqFwSNqs/D/7d19dFT1ncfx95dMnhMyCQ8JKxCSQIClHrBSgRaUqkUq\nZ0Vxj1sfuoKnW6mVZekq6rZnqa5nK0hdu7u6ttuWtbtSRJCIK7Ho4oHoNjyYIAFJICEoCIEkkOeE\nPMx3/5gLO4SQBEzmDsz3dU5OZu7MnXwyM/nOze/e3/2+dJcLia9e8YnRxET17bh5+cmTnGpo4I38\nfLYfPMjo4cPJTE2lpqmJkzU1JKWkcMdddxEbF8cXR47wQm4u3shI7p48mQenT2f5mjWUlZby/pYt\nfZrL2Ja+cdnevXt58cV/Ji8vjwMHivF6BzFt2n1MnHgbp04do66ukrS0LF588TuMTE8m/9Mfhdwk\npytdWuyPaWlpZ+V3v8s1KSlEeTzUNDYyNCmpy7N9dvh85753+HzEBhxVc/D4cZ7NyeGIc6LD9BEj\nePmVV/B6vZSWlpKamkpSUhJTpkw578O9o6ODl196iQ3r1rFn717GpqWRX1JCfX19j0fJhTMb3jFX\nLFWlrq6Ow4cPs2TJY+TlbSUzczR1dbXU1JzG4/Hg0zbuuX8SZQdqGZoWz4KF1/OVicPweq2D2aVq\nbm6juqqRxIHRXJf1PKeqG/FERBAfF0djUxPD0tI4VlHB1HHjGOH1UtnURFlFBZU1NTQ2+xuVR0VG\n0trWxnemT+eeqVNJjo/nX3JzeXPHDtra2i57tmxDQwOLFy2itLSUrXl5fflrX3Ws6Jurhs/nO3cG\nyrNjwuXl5fzsuWe59ZbbOF5xnCeWLsWbHMszz88mJjaS2ppm0jNSyBw9iIFJMXg8EcTEeGhoOENS\nkn0wVFU2sPa1QrZ/9Dn/8+4Bzpxpp62tnQfnP8CK5T+nsrKSYcOGER8fT3R0NLW1teTm5lJWVsao\nUaPIysoiLS2NoUOHEhERQXR0NNu3b+fuefM4XV3N+PR0Cg4eZOljj/HcihVu/7phoV+LvtMFaxdw\nNODcO4uAR4B24B1VfbKL9Q4DtYAPaOvqvPvO/azom0ty9OhRHl74kHMct5KQkEhBQQHVVadoamqh\nra0Nn0+JjPQwZuwwvjIxleiYASxcPI3scUMveLz29g72FVWQNWYwCQldH157paqvP8O0Cb9g/PhJ\n7Nq5i+997/uM/9NxvPHG67y9cdOXmuRXVlbG22+/TVpaGnPmzCExMbEPk5vu9HfRXwJcDwx0euR+\nE3gKuF1V20VksKpWdbHeIeB6VT3dw+Nb0Td9qra2ltjYWESEnTt3sm/fPj76KI+cnA3ExUdx/4JJ\n1Na0kpAYhc+n/GLFB+fWvX/BZLLGDOZQaQ3XDB/IN2aOYviIJDKyBvXJ3IJge+gv1hIfM5ZX/+O/\nrsj8pmv9VvSddomrcNolOkX/deCXqtrt7nURKQcmq2p1D/ezom+CoqOjg61bt7Ju/VoyMrKor6/j\n+PGjLJj/V0yYMIEtW7ZQXFJMbe1pfD7Y+NabJHkT2bF9N7Gx0XiTE5g1Zyz3PHAtSd5YEhKjGZVx\n/qGFra3tVFU2MjApxvX/Gv43r5zbb3yFqqqqSz5Ntglt/Vn0L2iXKCKFwFvAbKAZeFxVd3Wx7iH8\nTdM7gF+papddxa3om1C3f/9+Bg0aRH19Pb9fs5oNG9bS2trKyRNVNDY2kZrmZfLUEbyTU4THE4HH\nE0lzcwvXTU5n4d/cwG1zxhEZGeFv3RjRf411mppa2fxOMXsKj5HkjSF3437GjZnBqlW/67efadwR\nrHaJZ7f0i4AtqrpYRL4GvK6qmV2sP0xVj4vIEOA94FFV/bCL++mygMkhweicZUxf6OjooLm5mfz8\nfIqLi5k7dy5RUVGkpqbS0tLChpwNLF/+LPX1p0hOieOTws8YMXIIN96SwX3zJzHl66P6LEt7eweP\nzF/H2tcKGDMmi+kzvsHgwUNY+PAjZGZe8OdprjCdO2c9/fTT/VL0/xF4AP/O2lggEXgTGAwsV9Wt\nzv1KgSndDeOIyDKgXlVf6OI229I3V7Vt27ZRU1PDrFmzKC8vZ/2b63nh5yvYWvhD4uIiqa87w4h0\nL7U1LTQ2tjIyPRlV5URFPaUHqqiva+Fb3x6Lx3PhTte2tg6KPz3BjEn+TmVnW2Kaq1tQO2eJyMPA\nn6jqMhHJBt5T1fRO948DBqhqg4jEA5uBp1V1cxePbUXfhJ1n/mEZK59fyYABA1AFj0dob/cRGRlJ\nckoslSfriIqKInvsaD76cCfXTR7Okz+9lZ35R7lh2ghm3jqa22/8DYW7PiMj4xpKSz+npKSE7Oxs\nt381EwTBLvqRwG+BScAZZ/nWwHaJIpIBbAAU/ykfXgvFdonGuKm6upqoqCgSEhI4dOgQ6enptLS0\nUFRURHZ29rmdrx9//DG//NVLbHpnE9+8+WZyN+Xy1DMzeeyHObS3t4d8bwXT92xyljFh5GxntszM\ndMrKDrsdx7jAir4xYURVWb9+PTfddBNDhgxxO45xgRV9Y4wJI5dT9PvvYGFjjDEhx4q+McaEESv6\nxhgTRqzoG2NMGLGib4wxYcSKvjHGhBEr+sYYE0Z6XfRFZICIFIjIxoBli0Rkv4gUicjFTq8wW0SK\nReSAiDzRF6GNMcZcnkvZ0l8MfHr2itM568+Aa1X1WmBl5xWcFov/CtwGTADuFZFxXypxEAWewjRU\nWKbeCcVMEJq5LFPvhGKmy9Grou90zrod+HXA4oXAc6raDtBVq0TgBuCgqn6mqm3AGmDul4scPKH4\nIlum3gnFTBCauSxT74RipsvR2y39fwIex3+2zLOygRtFJF9EPhCRyV2sdw1wJOD6UWeZMcYYF/RY\n9J3OWSdUdTcQeI4HD5CsqlOBpcDa/olojDGmr/Rr5ywRmQr8VFVnO9efBFRVl3fxc+xsa8YYc4lC\nrXNWBFAC3AIcB3YA96rq/ksJaYwxpm98meP0fwtkOg3SVwN/Cf5G6CLy3wCq2gE8ir9N4j5gjRV8\nY4xxT8icT98YY0z/c31GrogsdiZ3FYnIX7uY4zcickJE9gQsSxaRzSJSIiJ/EJGkEMj05yKyV0Q6\nROSrwczTTaYVziS93SKyXkQGhkCmZ0TkExEpFJF3RSTN7UwBt/2tiPhEJCWYmS6WS0SWichRZ/Jl\ngYjMdjuTs7zHyZ/BzCQiawKeo3IRKQiBTBNF5I/O+3zHRY6iPJ+quvaFf8LWHiAaiMA/DJTpUpbp\n+Ju87wlYthxY6lx+Av+8BLczjQXGAFuAr4bI83QrMMC5/BzwsxDIlBBweRHwb25ncpYPB94FyoGU\nEHn9lgE/CnaWHjLNdOqBx7k+2O1MnW5fCfzE7UzAH4BZzuVvAx/09Dhub+mPB7ar6hn1j/9vA+a5\nEURVPwROd1o8F3jVufwqcKfbmVS1RFUPcv7hs25nel9Vfc7VfPyFze1MDQFX4wEfQXSR9xP8/5wX\nV3STy5X3E1w00w/oefJnsDMFugf4fZDiABfN5APOjkB4gS96ehy3i/5eYIYzjBKHf9bvCJczBRqq\nqicAVLUCGOpynivBQ0Cu2yEARORZEfkcuA/4+xDIcwdwRFWL3M7ShUed4blfB3sY8yJ6M/nTFSIy\nA6hQ1TK3swBLgJXO+3wF8FRPK7ha9FW1GP8QynvAJqAQ6HAzUw9sr3c3ROTHQJuqrnY7C4Cq/kRV\nRwKv4R/icY2IxAJ/h38o5dxil+J09jL+YdVJQAXwgst5ILQnf95LkLfyu/EDYLHzPl+C/6jKbrm9\npY+qrlLVyao6E6gBDrgcKdAJEUkFcHYEnnQ5T8gSkfn4/1O7z+UoXVkN3O1yhixgFPCJiJTjHwL7\nWERc/+9RVSvVGRQG/h34mpt5HEfwTwJFVXcCPhEZ5G6kc3OP5gGvu53F8aCq5gCo6jr85zvrlutF\nX0SGON9HAnfh/wN1LQ7nb31tBOY7lx8E3gp2IC7M1Pk2N5yXyTna43HgDlU9EyKZRgfcdifgxvyQ\nc5lUda+qpqlqpqpm4D8P1XWq6saGROfnKvDIpnn4h12DrfP7PAe4GcCZ/BmpAbP9XcoE8C1gv6oe\nC3KWszpn+sKZNIuI3EJvNpqDuff5Inukt+F/kxUCM13MsRo4BpwBPgcWAMnA+/hnFW8GvCGQ6U78\nW0HN+Gc554ZApoPAZ0CB8/VyCGRaBxQBu/F/WA9zO1On2w/hztE7XT1Xv8N/FN1u/MU2NQQyeYD/\ndF7DXcBNbmdylq8Cvh/s162b5+nrzvNTCPwR/4ZEt49jk7OMMSaMuD68Y4wxJnis6BtjTBixom+M\nMWHEir4xxoQRK/rGGBNGrOgbY0wYsaJvjDFhxIq+McaEkf8DIqDwGKJMiLsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1196e4e90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import geopandas as gp\n",
    "%pylab inline\n",
    "austria_shp = gp.read_file('austria.shp')\n",
    "austria_shp.plot()\n",
    "austria = pd.read_csv('http://dl.dropbox.com/u/8649795/AT_Austria.csv')\n",
    "austria.drop(['Oi2007', 'Dj2007', 'OrigAT11', 'OrigAT12', 'OrigAT13', 'OrigAT21', 'OrigAT22', 'OrigAT31', 'OrigAT32', 'OrigAT33', 'OrigAT34', 'DestAT11', 'DestAT12', 'DestAT13', 'DestAT21', 'DestAT22', 'DestAT31', 'DestAT32', 'DestAT33', 'DestAT34', 'beta', 'Offset'], axis=1, inplace=True)\n",
    "austria.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **Origin** and **Destination** columns refer to the origin, $i$, and destination, $j$, location labels, the **Data** column is the number of flows, the **Oi** and **Dj** columns are the number of total out-flows and total in-flows, respectively, and the **Dij** column is the distance between $i$ and $j$. In this case we use the total out-flow and total in-flow as variables to describe how emissive an origin is and how attractive a destination is. If we want a more informative model we can replace these with application specific variables that pertain to different hypotheses. Next, lets format the data into arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "austria = austria[austria['Origin'] != austria['Destination']]\n",
    "flows = austria['Data'].values\n",
    "Oi = austria['Oi'].values\n",
    "Dj = austria['Dj'].values\n",
    "Dij = austria['Dij'].values\n",
    "Origin = austria['Origin'].values\n",
    "Destination = austria['Destination'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Oi and Dj vectors need not be $n^2 \\times 1$ vectors. In fact, they can be $n^2 \\times k$ where $k$ is the number of variables that are being used to describe origin or desitnation factors associated with flows. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Calibrating the models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, lets import the main SpInt functions and calibrate some models. The main SpInt functions are found within the gravity namespace of the SpInt module ans the estimated parameters can be accessed via the **params** attribute of a successfully instantiated spatial interaction model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pysal.contrib.spint.gravity import Gravity, Production, Attraction, Doubly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unconstrained (basic gravity) model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.44314667  0.51739961 -0.00979932]\n"
     ]
    }
   ],
   "source": [
    "gravity = Gravity(flows, Oi, Dj, Dij, 'exp')\n",
    "print gravity.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Production-constrained model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.90285448 -0.0072617 ]\n"
     ]
    }
   ],
   "source": [
    "production = Production(flows, Origin, Dj, Dij, 'exp')\n",
    "print production.params[-2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Attraction-constrained model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.90037216 -0.00695034]\n"
     ]
    }
   ],
   "source": [
    "attraction = Attraction(flows, Destination, Oi, Dij, 'exp')\n",
    "print attraction.params[-2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doubly-constrained model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.00791533]\n"
     ]
    }
   ],
   "source": [
    "doubly = Doubly(flows, Origin, Destination, Dij, 'exp')\n",
    "print doubly.params[-1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that for the constrained models we have limited the params attribute to print only those associated variables (i.e., not fixed effects), though it is still possible to access the fixed effect parameters too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.16851884  0.52128801  0.98284063 -0.56934181 -0.28515686  0.0381801\n",
      " -0.47906115 -0.0141766  -0.1583821   0.90285448 -0.0072617 ]\n"
     ]
    }
   ],
   "source": [
    "print production.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first parameter is always the overall intercept with the subsequent 8 parameters representing the fixed effects in this case. You might ask, \"why not 9 fixed effects for the 9 different municipalities?\". Due to the coding scheme used in SpInt, and many popular statistical programming languages, you would use $n - 1$ binary indicator variables in the design matrix to include the fixed effects for all 9 municipalities in the model. While the non-zero entries in these columns of the design matrix indicate which rows are associated with which miunicipality, where a row has all zero entries then refers to the $nth$ municipality that has been left out. In Spint this is always the first origin or destination for the production-constrained and attraction-constrained models. For the doubly-cosntrained model, both the first origin and the first destination are left out. In terms of interpetting the parameters, these dropped locations are assumed to be 0. Since the fixed effects parameters are interpretted as deviations from the overall intercept, this essentially means the intercept acts as the fixed effect for the first location. Therefore, we could also say first 9 parameters are the origin fixed effect parameters and that the last two parameters are for destination attractiveness and distance, respectively. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Interpretting the parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Assessing model fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.812858632041\n",
      "0.910155702595\n",
      "0.909354566583\n",
      "0.943539912198\n",
      "\n",
      "0.834893355219\n",
      "0.913167247331\n",
      "0.912363460094\n",
      "0.946661920896\n",
      "\n",
      "0.68757357636\n",
      "0.740913740348\n",
      "0.752155260541\n",
      "0.811852110904\n",
      "\n",
      "1.01757072219\n",
      "0.464519773826\n",
      "0.58404798182\n",
      "0.37928618533\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print gravity.pseudoR2\n",
    "print production.pseudoR2\n",
    "print attraction.pseudoR2\n",
    "print doubly.pseudoR2\n",
    "print ''\n",
    "print gravity.D2\n",
    "print production.D2\n",
    "print attraction.D2\n",
    "print doubly.D2\n",
    "print \"\"\n",
    "print gravity.SSI\n",
    "print production.SSI\n",
    "print attraction.SSI\n",
    "print doubly.SSI\n",
    "print \"\"\n",
    "print gravity.SRMSE\n",
    "print production.SRMSE\n",
    "print attraction.SRMSE\n",
    "print doubly.SRMSE\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 Local models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.6 Testing for overdispersion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Additional functionality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Existing features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Future additions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, there are several paradigms for incorporating spatial effects into spatial interaction models (competing destinations, spatial autoregressive, eigenvector spatial filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python (spint)",
   "language": "python",
   "name": "spint"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}