examples/17_advanced_discretionary_locations.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e225bed3",
   "metadata": {},
   "source": [
    "# Applying a location model to discretionary activities\n",
    "\n",
    "This notebook applies a time-prism approach for choosing the location of non-mandatory activities. With this approach, agents try to sensibly combine the location of their discretionary activities around the locations of “anchor points” (defined as the locations of mandatory activities – home, employment, education).\n",
    "\n",
    "The `pam.planner.choice_location.DiscretionaryTrips` class allows the user to apply the time-space prism approach on a PAM plan, given accessibility and attraction data.\n",
    "\n",
    "The typical workflow goes as follows:\n",
    "\n",
    "``` python\n",
    "model = DiscretionaryTrips(population, od) # initialize the model and point to the data objects \n",
    "model.update_plan() # apply the model and update the location of each non-mandatory activity in the plan\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "82e63bf7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:27.464569Z",
     "iopub.status.busy": "2024-04-05T15:22:27.464373Z",
     "iopub.status.idle": "2024-04-05T15:22:28.529618Z",
     "shell.execute_reply": "2024-04-05T15:22:28.529206Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/6n/0h9tynqn581fxsytcc863h94tm217b/T/ipykernel_95220/2386317155.py:6: DeprecationWarning: \n",
      "Pyarrow will become a required dependency of pandas in the next major release of pandas (pandas 3.0),\n",
      "(to allow more performant data types, such as the Arrow string type, and better interoperability with other libraries)\n",
      "but was not found to be installed on your system.\n",
      "If this would cause problems for you,\n",
      "please provide us feedback at https://github.com/pandas-dev/pandas/issues/54466\n",
      "        \n",
      "  import pandas as pd\n"
     ]
    }
   ],
   "source": [
    "import logging\n",
    "import random\n",
    "from copy import deepcopy\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from prettytable import PrettyTable\n",
    "from shapely.geometry import Point\n",
    "\n",
    "from pam.activity import Activity, Leg, Plan\n",
    "from pam.location import Location\n",
    "from pam.planner.choice_location import DiscretionaryTripOD, DiscretionaryTrips\n",
    "from pam.planner.od import ODFactory, ODMatrix\n",
    "from pam.planner.utils_planner import get_trip_chains_either_anchor\n",
    "from pam.utils import minutes_to_datetime as mtdt\n",
    "from pam.variables import END_OF_DAY\n",
    "\n",
    "logging.getLogger(\"pam\").setLevel(logging.DEBUG)\n",
    "random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f8b4b6b-24a2-4f12-a821-7e0966393d3d",
   "metadata": {},
   "source": [
    "## Methodology"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a4a3e18-afa5-4775-8bda-d466b1ace936",
   "metadata": {},
   "source": [
    "We start by identifying fixed mandatory locations (ie home, work, education), splitting the PAM plan into chains of activities."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e889c39-de05-4efa-9ecf-1066f7df7160",
   "metadata": {},
   "source": [
    "We then distinguish between 4 possible cases of trip chains:"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "9923e0b2-e543-4f78-a8d4-a05a3edb815f",
   "metadata": {},
   "source": [
    "![trip chain cases](../resources/example-discretionary-cases.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "535cc558-a480-43b3-869f-8301bcb53122",
   "metadata": {},
   "source": [
    "**Case 1** is a simple tour between mandatory activities, and may be approached with a discrete choice model (see 15_advanced_choice_modelling.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec8a8900-760c-4bcd-858f-d5d076626657",
   "metadata": {},
   "source": [
    "**Case 2** selects a possible destination using origin-destination probabilities from the anchor location (saved under the `od_probs` variable of the `od` object)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "748c08f9-0009-41c3-b8fc-63f3310649cd",
   "metadata": {},
   "source": [
    "**Case 3** weights candidate locations as a function of three elements:\n",
    "* leg ratio, which captures how likely it is for the discretionary location to be closer to the previous or to the next fixed location (between 0 and 1)​\n",
    "* diversion factor, which captures the deviation from the straight path the individual is likely to accept in order to reach a discretionary activity location\n",
    "* attraction (OD probabilities), which captures the number of opportunities in a destination\n",
    "\n",
    "The final destination probabilities are calculated by multiplying the three elements: $p_j = p\\_leg_{ij} * p\\_diversion_{ij} * p\\_attraction_{ij}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14d163b9-b914-4576-aef9-db4087b38126",
   "metadata": {},
   "source": [
    "**Case 4** solves the first location as case 2, and then the rest as case 3."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d863cc2c",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6d13dd57",
   "metadata": {},
   "source": [
    "Let's create a plan with the following activity sequence: home->work->shop->work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a02cfd1e-33d8-4760-a368-f59d10727b38",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.532087Z",
     "iopub.status.busy": "2024-04-05T15:22:28.531837Z",
     "iopub.status.idle": "2024-04-05T15:22:28.691053Z",
     "shell.execute_reply": "2024-04-05T15:22:28.690570Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAADTCAYAAADJTpGsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMaUlEQVR4nO3dd3QUZd/G8WvTCwkhhJAECL2EGkLoSC+CdKVIF6yANEUEUVBQEEQQeEB6EXnAgqh0lCTCC0LoIAhIC0jo0kJP5v0jD6sxoSc7bPb7OSfn7M7O7l5z7z2be+a3M2MxDMMQAAAAAAAAAADAE87J7AAAAAAAAAAAAAAPgqIGAAAAAAAAAACwCxQ1AAAAAAAAAACAXaCoAQAAAAAAAAAA7AJFDQAAAAAAAAAAYBcoagAAAAAAAAAAALtAUQMAAAAAAAAAANgFihoAAAAAAAAAAMAuUNQAAAAAAAAAAAB2gaIGAAAAANxFvnz5NG7cOLNjAAAAAPgfihoAAABABunSpYssFossFotcXV1VoEABvfnmm0pISDA72kPr0qWLmjdvnuL+P5ctZ86cqlevnmbOnKmkpCTzggIAAADI1ChqAAAAABno6aefVnx8vA4dOqThw4dr0qRJevPNNx/ptQzD0O3bt9M54aO7s2xHjhzR8uXLVatWLfXu3VuNGzd+onICAAAAyDwoagAAAAAZyN3dXUFBQcqTJ4/atWun9u3ba/HixZKSixSjRo1SgQIF5OnpqTJlyuibb76xPjc6OloWi0UrV65UZGSk3N3dtXbtWu3YsUO1atWSj4+PfH19Va5cOW3evNn6vG+//VYlSpSQu7u78uXLpzFjxqTIlC9fPn300Ufq2rWrfHx8FBoaqqlTpz7ysuXKlUsREREaNGiQvv/+ey1fvlyzZ89+qNcaOnSoQkND5e7urpCQEPXq1cv62Lx58xQZGSkfHx8FBQWpXbt2On36dJrtVLZsWXl6eqp27do6ffq0li9frrCwMPn6+ur555/X1atXrc+rWbOmevbsqZ49e8rPz0/Zs2fX4MGDZRjGXXNevHhRL7/8sgIDA+Xr66vatWtrx44dD7WsAAAAAB4dRQ0AAADAhjw9PXXr1i1J0uDBgzVr1ixNnjxZv/32m/r27asOHTooJiYmxXPeeustjRgxQnv37lXp0qXVvn175c6dW7GxsdqyZYvefvttubq6SpK2bNmi1q1bq23bttq1a5eGDh2qd999N1WRYcyYMYqMjNS2bdvUvXt3vfbaa/r9998fe/lq166tMmXKaNGiRdZpXbp0Uc2aNe/6nG+++UZjx47VlClTdODAAS1evFilSpWyPn7z5k0NGzZMO3bs0OLFi3X48GF16dIl1esMHTpUEydO1Pr163Xs2DG1bt1a48aN0/z587V06VKtXr1aEyZMSPGcOXPmyMXFRRs3btT48eM1duxYTZ8+Pc2chmHomWee0cmTJ7Vs2TJt2bJFERERqlOnjs6fP/9wDQUAAADgkbiYHQAAAABwFJs2bdL8+fNVp04dJSQk6NNPP9WaNWtUuXJlSVKBAgW0bt06TZkyRTVq1LA+74MPPlC9evWs9+Pi4tS/f38VK1ZMklS4cGHrY59++qnq1Kmjd999V5JUpEgR7dmzR6NHj05RCGjUqJG6d+8uSRowYIDGjh2r6Oho62s+jmLFimnnzp3W+8HBwfe8zkZcXJyCgoJUt25dubq6KjQ0VBUqVLA+3rVrV+vtAgUKaPz48apQoYKuXLmiLFmyWB8bPny4qlatKknq1q2bBg4cqIMHD6pAgQKSpOeee05RUVEaMGCA9Tl58uTR2LFjZbFYVLRoUe3atUtjx47VSy+9lCpnVFSUdu3apdOnT8vd3V2S9Mknn2jx4sX65ptv9PLLLz9sUwEAAAB4SBypAQAAAGSgJUuWKEuWLPLw8FDlypVVvXp1TZgwQXv27NH169dVr149ZcmSxfo3d+5cHTx4MMVrREZGprjfr18/vfjii6pbt65GjhyZYv69e/dad+zfUbVqVR04cECJiYnWaaVLl7betlgsCgoKSnFKp8dhGIYsFov1/ogRIzR37ty7zt+qVStdu3ZNBQoU0EsvvaTvvvsuxTU5tm3bpmbNmilv3rzy8fGxHvURFxeX4nX+uUw5c+aUl5eXtaBxZ9q/l7FSpUopslauXDlVW92xZcsWXblyRdmzZ0/xmR0+fDjVZwYAAAAgY3CkBgAAAJCBatWqpcmTJ8vV1VUhISHW00QdPnxYkrR06VLlypUrxXPuHAVwh7e3d4r7Q4cOVbt27bR06VItX75cQ4YM0YIFC9SiRYtUBQVJaV4j4k6OOywWyz2PpngYe/fuVf78+R94/jx58mjfvn1avXq1fvrpJ3Xv3l2jR49WTEyMbt68qfr166t+/fqaN2+ecuTIobi4ODVo0EA3b9686zJZLJZ0X8akpCQFBwcrOjo61WN+fn6P/LoAAAAAHhxFDQAAACADeXt7q1ChQqmmFy9eXO7u7oqLi0txqqkHVaRIERUpUkR9+/bV888/r1mzZqlFixYqXry41q1bl2Le9evXq0iRInJ2dn7k5XhQa9as0a5du9S3b9+Hep6np6eaNm2qpk2bqkePHipWrJh27dolwzB09uxZjRw5Unny5JGkFBdFf1y//vprqvuFCxdOs60iIiJ08uRJubi4KF++fOmWAQAAAMCDo6gBAAAAmMDHx0dvvvmm+vbtq6SkJFWrVk2XLl3S+vXrlSVLFnXu3DnN5127dk39+/fXc889p/z58+v48eOKjY3Vs88+K0l64403VL58eQ0bNkxt2rTRhg0bNHHiRE2aNCndl+HGjRs6efKkEhMTderUKa1YsUIjRoxQ48aN1alTJ+t8AwcO1J9//nnXU1DNnj1biYmJqlixory8vPTFF1/I09NTefPmVVJSktzc3DRhwgS9+uqr2r17t4YNG5Zuy3Ds2DH169dPr7zyirZu3aoJEyZozJgxac5bt25dVa5cWc2bN9fHH3+sokWL6sSJE1q2bJmaN2+e6jRhAAAAANIfRQ0AAADAJMOGDVNgYKBGjBihQ4cOyc/PTxERERo0aNBdn+Ps7Kxz586pU6dOOnXqlAICAtSyZUu9//77kpKPJvjqq6/03nvvadiwYQoODtYHH3yQ4iLh6WXFihUKDg6Wi4uLsmXLpjJlymj8+PHq3LmznJz+vnxffHx8qutf/JOfn59Gjhypfv36KTExUaVKldKPP/6o7NmzS0ouegwaNEjjx49XRESEPvnkEzVt2jRdlqFTp066du2aKlSoIGdnZ73++ut3veC3xWLRsmXL9M4776hr1646c+aMgoKCVL16deXMmTNd8gAAAAC4N4uR1gl2AQAAACCTq1mzpsLDwzVu3DizowAAAAB4QE73nwUAAAAAAAAAAMB8FDUAAAAAAAAAAIBd4PRTAAAAAAAAAADALnCkBgAAAAAAAAAAsAsUNQAAAAAAAAAAgF2gqAEAAAAAAAAAAOwCRQ0AAAAAAAAAAGAXKGoAAAAAAAAAAAC7QFEDAAAAAAAAAADYBYoaAAAAAAAAAADALrg86hNz586dnjkAwG6dPn1aiYmJcnZ2VmBgoNlxTEEbAAAAACkxRqYNJNpAog1AH5BogzuOHz+eLq9jMQzDeKQnWizpEgAAAAAAAAAAAGRuj1iKSOWRj9SQJCcni4JzeqdLEACwV3/GX7He9vDwMDGJea5fvy5JcrI4yT+Lv8lpAAAAAPOdvXzWetvhtxOcnJQjRw6T05jj1KlT1tv0AycFBzjuL9Qd2Z+nT1pvsx447v70+FMJ6fZaj1XUCM7preO7u6VXFgCwS845xispyZCHh4fq1atndhxT/Pjjj5Ik/yz+WvDGQpPTAAAAAOar/349JRlJbCdIypEjh3755ReT05gjLCxMSUn0A0kKDgjU0eVbTE4DM7iVz8N6cGc9cOD96blLzki31+JC4QAAAAAAAAAAwC5Q1AAAAAAAAAAAAHaBogYAAAAAAAAAALALFDUAAAAAAAAAAIBdoKgBAAAAAAAAAADsAkUNAAAAAAAAAABgFyhqAAAAAAAAAAAAu0BRAwAAAAAAAAAA2AWKGgAAAAAAAAAAwC5Q1AAAAAAAAAAAAHaBogYAAAAAAAAAALALFDUAAAAAAAAAAIBdoKgBAAAAAAAAAADsAkUNAAAAAAAAAABgFyhqAAAAAAAAAAAAu0BRAwAAAAAAAAAA2AWKGgAAAAAAAAAAwC5Q1AAAAAAAAAAAAHaBogYAAAAAAAAAALALFDUAAAAAAAAAAIBdoKgBAAAAAAAAAADsAkUNAAAAAAAAAABgFyhqAAAAAAAAAAAAu0BRAwAAAAAAAAAA2AWKGgAAAAAAAAAAwC64mB0AAHB/FSpUeOB5N23alIFJAAAAADxJ2FYAfQBIxrrgOChqAIAdGDRo0APNZxiGWrRokcFpAAAAADwp2FYAfQBIxrrgOChqAIAdaN68udkRAAAAADyB2FYAfQBIxrpwdzdvJurw0UsqmD+rXFzs/4oU9r8EAODAXF1dzY4AAAAA4AnEtkJqhmGYHcGm6ANAMkdeF65evaVuvVbLK/d/VKLqF4o7flmS1OvtaI0cF2tyukdHUQMA7IyTk5PatGmjWbNm6auvvlLOnDklSe3bt1e9evVMTgcAAADALGwrSFOmTElzemJiot544w0bp7E9+gCQjHUh2cBh67Vj91lF//CsPDz+PmlT3RqhWrj4gInJHg9FDQCwM61atVLt2rU1e/Zs3b592zr9yJEjDvWPGQAAAEBKbCtIc+bM0cKFC1NMS0xMVN++fbV3716TUtkOfQBIxrqQbPGyg5r4cU1Vq5RLln9ML17UXwcPXzQt1+OiqAEAdqZ27dr6z3/+o5iYGCUmJlqnHzlyRLlz5zYxGQAAAAAzsa0gTZ06VZ988omWLVsmSbp165Z69+6tP/74Q3PnzjU5XcajDwDJWBeSnTl3TYE5vFJNT7h6SxZLGk+wE1woHADsjL+/v+Lj41NNd3JykrOzswmJAAAAADwJ2FaQSpYsqYkTJ6p79+5yc3PTN998o2PHjmnu3LkKCAgwO16Gow8AyVgXkpUvm1NLVx3W6y+HS5K1kDFt7m5VLh9sXrDHxJEaAGBn4uLiVKJEiVTTq1atqkOHDpmQCAAAAMCTgG2FZBUrVtSoUaPUu3dv/fnnn/riiy8coqAh0QeAO1gXko14t4re+XCDXntjjW4nJumzKdtVr+UizV6wVx++U9nseI+MIzUAwM4sWLBAffv2lb+/v5ycnFS5cmXlypVLtWvX1rBhw8yOBwAAAMAkjrqt0LNnzzSnZ8uWTb6+vnrvvfes0yZOnGirWKZw1D4A/BvrQrIqFUL0f8tb6ZOJW1Uwn59WRcUponQObVjRWqWK22+xl6IGANiZ2NhYjR49Wq1atZJhGGrfvr0OHjyoYcOGafv27WbHAwAAAGASR91W8PHxSXP6U089ZeMk5nPUPgD8G+vC30oVD9CcSfXNjpGuKGoAgB3atm2btm3bZnYMAAAAAE8YR9xWGDFihCTJMAydOHFC/v7+8vT0NDmVeRyxDwBpYV342+kzV3X67FUlJRkpppcukcOkRI+Ha2oAgJ3p1auXSpcubXYMAMATYE7UHL0y+eV0nxcAYJ8cfVvBMAw1aNBAp06dMjuKaRy9D5jJpVwufR+1wuwY+B/WhWRbtp9SyarzFFx8uko/9aXCa8y3/pWt+V+z4z0yjtSAKbr0WKULF29o8bwmKaZHrzuuWs2+1V+HXpVfVneT0gFPNh8fH7333nu6fPmy1q5dq+joaIe6yBUAOIrB89/RjVs3Nbrz6FSP7Tn2m3rN6KVPu4xV84rNH+j1WldpnWLeUd99rCvXr+iD5x3nnMIAkNk5+raCk5OT8ubNqwsXLpgdxTSO3gcy0unzZzVk8iit+L8onTp/Vtl8s6p04eJ675V+qlw60ux4+BfWhWQvvL5aRQr6acZndZUz0EsWi9mJ0gdFDQCwMx9++KG8vb1VrVo1Va9eXU2aNNGff/6pmJgYxcTE6PTp02ZHBACkg4ZlG2rowqE6deGUcvrlTPHY8m0rVDCokErne/Bfn3m6e8pTjnsqDgBwBGwrSP3799eoUaM0dOhQFSlSxOw4NkcfyDit+r+kW7dvaeb741Qgd16dOndGazat0/mLF8yOhjSwLiQ7fPSSFs1prEIF/MyOkq4oauCJ9u0PB/TeyF/1x+GLCs7ppddfCtcbPSKsj+cLn6kXO5TU/oN/adHSg8qezUPjR9RQlQrBerH3z/p57THlD/XVrAn1FFn2750B6zed0Nsf/J9it51SgL+nWjxTUCPerSpvb1czFhN4aAkJCVq5cqVWrlyp7Nmzq3r16qpbt67atWunFi1amB0PAJAOKhWpLD9vP63cvlKdanayTr9+87pidkera51umhM1R+t//z9NeW2qJGn74e2atnqqjp45KmcnZ+ULzKdBz76jnH45U8w7J2qOVu1YJUmqO7SOJOmTzmMUnj/c5ssJAEhfjr6t0L9/f127dk3NmjWTq6urPDw8Ujy+adMmk5LZjqP3gYxw4fJF/d/2Tfp56jeqUa6yJClvcG5VKFk2xXxnL5zXs29006oN0coVGKTRfYeoSY2/L9Acs2WDBnw2XDv375G/r586Nm6lYd3fkotL8i7a2i8/p5IFi0qSvly2SM7OTnrluU764LW3ZMksP7G3IdYFqU71PNqx+wxFDcBWtmw/pdbdlmvogIpq07yI1m+KV/e3opQ9m4e6tCtunW/s59v00eAqevfNCho7eZs6vrZKVSsGq2v7Ehr9fjUNeP//1Kn7Kv22voMsFot27TmrBs8t1rBBlTXjs7o6c+6aeg6IVs8BUZo1sf49EgFPHmdnZxUqVEhFihRRYGCgQx9mDQCZjbOzs+qVqadV21eqY42O1g3ZX/bE6FbiLdUpXUeLfl1knT8xMVFDFrynRuWe0TvPDdatxFva9+fvab526yqtFXc2TldvJKh/s7ckST6ePhm/UAAAm3HUbYVBgwaZHeGJ4ah9ICNk8fRWFi9vfR+9QpVKRcjdLe1Tpg+b9qlG9hqsj/sM1n8WzFLHwT11aMlG+WfNpj9Px6tJr47q3KS1Zr//mfYd+UOvDO8vD3d3DXnlDetrzF3ytbo2a6v1c37Ulj079eqHbylvUG692LK9rRY303HkdWH6Z3XVuccq7f79nEoWyy5XV+cUjzdtWMCkZI+HogZMs2TVYWUJnZRiWmJikvX2p5O3qU71PHr3zYqSpCKFsmnP/vMaPXFLiqJGo7r59EqXUpKk9/pX1ORZu1S+bE61alZYkjSgVzlVfvornTp9VUE5vTV6wha1e66o+ryaXE0vXDCbxo+ooRpNvtXkT2rLw4PVAk++UqVKqXr16qpSpYqcnJy0YcMGDR8+XDt37jQ7GgAgHT1dtqG+Wv+VdhzZrvD8yWOXFdtWqFrYU6mKEAk3EpRwI0GVilRSiH+IJClvjrxpvq6nu6fcXdx06/ZN+fv4Z+xCAABsytG3FRzl19f34uh9ICO4uLho5tCxemX4W5r67TyVLVZS1SMqqU2DZipd+O99VJ0at1bbp5tLkob3fFsTF87Upt+26+kqtTT56znKkzNE4wd8KIvFomL5C+nEmZMaOOEjvftSXzk5OUmS8uQM0Zg33pfFYlHRfIW064/f9dn8aRQ1HgHrgrR+U7zW/XpCy386kuoxi8WixDO9bB8qHbD3FqapVS23Jn9SO8W0jVtOqsOrKyVJe/efV7N/VQurVgjWuM+3KTExSc7OyV/2pUsEWB/PGeglSSpVPPW002eTixpbdpzWH4cv6stv9lnnMQxDSUmGDh+9pLCibNjjyTZz5kz5+Pho27ZtmjRpkjZt2qRbt26ZHQsAkAFCc4SqRJ4SWrFthcLzl9WJ8ye06+gujez4cap5fb181SC8gd7+YoDKFSyniAIRqlGiprL7ZDchOQDADGwrJEtMTNRPP/2kgwcPymKxqFChQqpdu7acnZ3v/2Q7Rx/IOC3rPKNG1epo7bZN+nXnFq3cEKVP5k7W1MGj1blpG0lS6cJh1vm9Pb3k45VFZ86flST9fvgPVSpdLsVppKqEl9eVqwk6fipeocG5JEkVS0WkmKdy6XIaO2+KEhMTHaIPpxfWhWS9BkarY+tievfNCsoZ6G12nHRDUQOm8fZyTXU+t+MnrlhvG4ZSnS/QMFK/jqurk/X2nfldXVJPS/rfQSBJSYZe6VxSvV4OT/Vaobk57QKefAsXLtS6deuUkJBgdhQAgA08XbahJi6boNcb9dKKbSuU0y+nIgpEpDlv/+ZvqXnFFor9I1bRu6M1a80sfdxxlIrnKZ7m/ACAzIVtBeno0aN6+eWXderUKeXPn1+GYWjKlCkKCgrS1KlTFRoaanbEDEUfyFge7h6qV6m66lWqrndf7quXP3hT708ZYy1quLqkvFarxWJR0v92SBmGkcZ+LuN/89kgvINhXUh27vx19X2tbKYqaEiS0/1nAcxRvKi/1v16IsW09bHxKlIwm/UojUcRUSZQv/1+ToUK+KX6c3Oj4o0n38qVKx3+nzIAOJKaJWrKyclJa3at0eodq9QgvME9LxRZOLiw2j3VTuNfnKB8gfm0ZteaNOdzcXZVkpGU5mMAAPvEtoI0fPhw5cmTR9HR0fruu++0ePFiRUVFKXfu3Bo+fLjZ8TIcfcC2wgoUVsK1qw8874Ydm62FDEnasGOzfLyzKFdgsHXaxl1bUzzv111bVTg0P0dpPCTWhWQtGxdS1LrjZsdIdxypgSfWG90jVL7uAg37ZKPaNC+iDbHxmjh9hyaNqvVYrzugVzlVavCVevSP0kudSsrby0V79/+l1dFxmvBxzfQJD2SwQoUKqVq1agoICJCra8pfgowYMcKkVACAjODp7qmaJWpq5s8zlHAjQfXDG6Q5X/xf8Vq6ZamqFK2s7D4BOnb2mI6fO656ZeqnOX+QX05tPhirY2ePydfTV94e3nJxZvMAAOydo28rxMbGauHChfLz87NOy5Ytm9588009//zz5gWzIUfvAxnh3IXzajPgFb3QrK1KFQ6Tj1cWbdmzQ5/MnaymNdMem/3ba606a/z86eo9arC6t35B+48e1PtTxqhP+5et19OQpGOnTuiNT4fq5ZYdtO333frPwpka3fe9jFq0TI11QSpS0E8Dh63Xul9PqFTxgBRnt5GkXq+EmxPsMbHVgidWRJlAfTWjod4b+auGfbJJwTm99cHblVJcJPxRlC6RQzE/PKt3Ptygp575Woakgvmyqk3zwukTHMhgTz31lPr06aPt27crPDxc27ZtU0hIiPz8/PTrr7+aHQ8AkAEaRjTU8m3LVa5gpHL65UxzHg9Xdx07G6f3t6/SpWuX5J/FX80rNFfjco3TnL9RuWe048gOdZ/6mq7dvKZPOo9ReP7wDFwKAEBGY1tBcnNzS/PX2QkJCal2amZG9IGMkcXLWxVKRuizL6fp4PGjunX7lvLkDFG35u00sOvrD/QauQKD9eP4LzTgs+Ga/nw9+fv66YVmz+udbr1TzNfxmed0/fp1Ve7UWM7OzurRpqteatkhIxYrU2NdSDZ93m/K4u2qmPV/Kmb9nykes1jst6hhMYy0rlLwAE+0WJQrOIuO7+6W3pkAwK445xivpCRDHh4eqlevXoa/3/jx47VixQotW7ZMCxYsUO/evXXq1Cn16NFD58+f13//+98Mz/BvP/74oyQpwCdAC95YaPP3BwAAAJ409d+vpyQjyWbbCdKTt61wZzshZ86c+uWXX2zynm+99Zb27NmjDz/8UKVLl5Yk7dixQ++++65KlCihkSNH2iTHHWFhYUpKsl0/eNL6gPR3P8gVGKSjy7fY/P3tSe2Xn1N4keL69M0PzI6SrtzK57HpeiA9eeuCdT1w4P3puUvO0PETl9PltbimBgDYmaCgIG3evFmSdPv2bXl4eEiSvv/+ezVo8GCHvQIAAADIfBx5W+Ho0aOSpMGDBytPnjxq06aNSpUqpVKlSqlt27YKDQ3VO++8Y3LKjOfIfQD4J9aF1AzD0CMe3/DE4fRTAGBnLl++LE9PT0nSuXPnlDdvXh09elTe3t5yd3c3OR0AAAAAszjytkKDBg2UM2dOVaxYUXXr1tVbb72lw4cPyzAMFSpUSHnz5jU7ok04ch8A/ol14W9zF+zV6IlbdODQBUnJ19no37OcOrYJMzfYY6CoAQB2Zs+ePQoPD9fRo0e1bt06vfjiiypVqpTCw8O1c+dOs+MBAAAAMIkjbyvMmzdPmzZt0qZNmzRs2DDduHFDwcHBqlSpkq5cuSIPDw/lzJn2dakyE0fuA5nBmqnfmB0h02BdSPbppK16d8QG9exWRlUrhsgwDP3fpni9+uYanT1/TX1fizA74iOhqAEAdmbKlClyc3OTJH3zzTdKTExUWFiYNmzYoIULuZ4FAAAA4KgceVshMjJSkZGR6t69u27duqXt27dr48aN2rRpk5YuXaqbN28qNDRUK1euNDtqhnLkPgD8E+tCsgnTdmjy6Nrq1PbvozKaNSqoEkX9NXTURooaAICM5+TkpAoVKmjr1q2Sks+HuGjRIpNTAQAAADAb2wp/c3V1Vfny5VWqVCmVLVtW69at01dffaW4uDizo2Uo+gCQjHXhb/GnElSlQnCq6VUqBCv+VIIJidIHFwoHADuSlJSk1157Ta6urmZHAQAAAPAEYVtBunHjhjZs2KBx48apXbt2Kl++vIYPH66EhAQNHTpU0dHRZkfMUPQBIBnrwt8K5ffTV4v3p5q+8LsDKlzAz/aB0glHagCAndm3b58KFCigM2fOmB0FAAAAwBPEkbcVOnTooF27dik0NFSRkZHq0KGDKlSooICAALOj2ZQj9wHgn1gXkr0/oKLavLhcv2w4oaoVg2WRRes2ntDPvxzTVzMamh3vkVHUAAA7s3z5cnXt2lUBAQH6448/dOPGjRSPHzlyxJxgAAAAAEzlyNsK27ZtU44cOVSxYkVVqFBBkZGR8vf3NzuWzTlyHwD+iXUh2bNNC2vjKl+NnbxNi5cdkmEYKl7EX5tWt1HZ0oFmx3tkFDUAwE706tVL06ZNU//+/SVJL730kvUxwzBksVhkGIZatGhhVkQAAAAAJmBbQYqNjdXmzZu1adMmTZs2Tf369VP+/PlVvnx5VahQQRUqVMjURQ76AJCMdSHZpUt/F3EKF/DTpNG10pzH19fdlrHSDUUNALATtWrV0pw5c1L8QwYAAAAAthUkLy8vVa9eXdWrV5ckXblyRVu2bNHGjRs1ffp0vfnmm8qXL5+WLFlictKMQR8AkrEuJPMr8LksFst950s808sGadIfRQ0AsBN3/hk5+vkgAQAAAKTEtkJqXl5e8vPzk5+fn7JmzSoXFxcdPHjQ7FgZhj4AJGNdSBb1/bPW24YhNWr7vaaPq6NcwVlMTJV+KGoAgB0xDMPsCAAAAACeQI6+rZCUlKTdu3dr48aN2rhxo7Zu3apr164pZ86cqlixot577z1VrFjR7JgZytH7AHAH64JUo2ruFPednSyqFBmsAvmympQofVHUAAA78vnnn9/3n3P79u1tlAYAAADAk8LRtxUiIyN17do15ciRQxUqVNDAgQNVsWJFhYaGmh3NZhy9DwB3sC5kfhQ1AMCOzJ8/X1evXjU7BgAAAIAnjKNvK7z11luqWLGi8ufPb3YU0zh6HwDuYF3I/ChqAIAdWbt2rS5evGh2DAAAAABPGEffVmjbtq3ZEUzn6H0AuIN1IW0PcN1wu0FRAwDsBOeEBAAAAJAWthVAHwCSsS4ka9lpSYr7128k6tU31sjbyzXF9EVzG9syVrqhqAEAdsKSmUrqAAAAANIN2wqgDwDJWBeSZfV1S3G/Q6uiJiXJGBQ1AMBONG/e3OwIAAAAAJ5AbCuAPgAkY11INmtifbMjZCgnswMAAAAAAAAAAAA8CIoaAAAAAAAAAADALlDUAAAAAAAAAAAAdoGiBgAAAAAAAAAAsAsUNQAAAAAAAAAAgF2gqAEAAAAAAAAAAOwCRQ0AAAAAAAAAAGAXKGoAAAAAAAAAAAC7QFEDAAAAAAAAAADYBYoaAAAAAAAAAADALlDUAAAAAAAAAAAAdoGiBgAAAAAAAAAAsAsUNQAAAAAAAAAAgF2gqAEAAAAAAAAAAOwCRQ0AAAAAAAAAAGAXKGoAAAAAAAAAAAC7QFEDAAAAAAAAAADYBYoaAAAAAAAAAADALlDUAAAAAAAAAAAAdoGiBgAAAAAAAAAAsAsUNQAAAAAAAAAAgF2gqAEAAAAAAAAAAOyCxTAM45GeaLHIycmi4Jze6Z0JAOzKn/FXrLc9PDxMTGKe69evS5KcLE7yz+JvchoAAADAfGcvn7XedvjtBCcn5ciRw+Q05jh16pT1Nv3AScEBgSangRn+PH3Sepv1wHH3p8efSlBiYlK6vNZjFTUAAAAAAAAAAADu5xFLEam4POoTc+XKlS4BAMDenT59WomJiXJ2dlZgoGP+6oQ2AAAAAFJijEwbSLSBRBuAPiDRBuntkY/UAAAAAAAAAAAAsCUuFA4AAAAAAAAAAOwCRQ0AAAAAAAAAAGAXHqmocePGDQ0dOlQ3btxI7zx2gzagDSTaQKINJNpAog0cffkl2kCiDSTaQKINJNpAog0k2sDRl1+iDSTaQKINJNpAog0k2sDRl1+iDaT0bYNHuqbGpUuXlDVrVl28eFG+vr6PHcIe0Qa0gUQbSLSBRBtItIGjL79EG0i0gUQbSLSBRBtItIFEGzj68ku0gUQbSLSBRBtItIFEGzj68ku0gZS+bcDppwAAAAAAAAAAgF2gqAEAAAAAAAAAAOwCRQ0AAAAAAAAAAGAXHqmo4e7uriFDhsjd3T2989gN2oA2kGgDiTaQaAOJNnD05ZdoA4k2kGgDiTaQaAOJNpBoA0dffok2kGgDiTaQaAOJNpBoA0dffok2kNK3DR7pQuEAAAAAAAAAAAC2xumnAAAAAAAAAACAXaCoAQAAAAAAAAAA7AJFDQAAAAAAAAAAYBcoagAAAAAAAAAAALvwSEWNSZMmKX/+/PLw8FC5cuW0du3a9M71xPrll1/UpEkThYSEyGKxaPHixWZHsrkRI0aofPny8vHxUWBgoJo3b659+/aZHcumJk+erNKlS8vX11e+vr6qXLmyli9fbnYs04wYMUIWi0V9+vQxO4rNDB06VBaLJcVfUFCQ2bFs7s8//1SHDh2UPXt2eXl5KTw8XFu2bDE7ls3ky5cvVT+wWCzq0aOH2dFs5vbt2xo8eLDy588vT09PFShQQB988IGSkpLMjmZTly9fVp8+fZQ3b155enqqSpUqio2NNTtWhrnfeMgwDA0dOlQhISHy9PRUzZo19dtvv5kTNoPcrw0WLVqkBg0aKCAgQBaLRdu3bzclZ0a51/LfunVLAwYMUKlSpeTt7a2QkBB16tRJJ06cMC9wBrhfHxg6dKiKFSsmb29vZcuWTXXr1tXGjRvNCZtBHmbb6JVXXpHFYtG4ceNsls8W7tcGXbp0STVOqFSpkjlhM8iD9IO9e/eqadOmypo1q3x8fFSpUiXFxcXZPmwGuV8bpDVetFgsGj16tDmBM8D92uDKlSvq2bOncufOLU9PT4WFhWny5MnmhM0g92uDU6dOqUuXLgoJCZGXl5eefvppHThwwJywGeBB9hVl9jHig7RBZh8j3q8NMvs48UH6QGYfIz7sfuNHHSM+dFFj4cKF6tOnj9555x1t27ZNTz31lBo2bJipBiT3kpCQoDJlymjixIlmRzFNTEyMevTooV9//VWrV6/W7du3Vb9+fSUkJJgdzWZy586tkSNHavPmzdq8ebNq166tZs2aZap/xg8qNjZWU6dOVenSpc2OYnMlSpRQfHy89W/Xrl1mR7Kpv/76S1WrVpWrq6uWL1+uPXv2aMyYMfLz8zM7ms3Exsam6AOrV6+WJLVq1crkZLbz8ccf6/PPP9fEiRO1d+9ejRo1SqNHj9aECRPMjmZTL774olavXq0vvvhCu3btUv369VW3bl39+eefZkfLEPcbD40aNUqffvqpJk6cqNjYWAUFBalevXq6fPmyjZNmnPu1QUJCgqpWraqRI0faOJlt3Gv5r169qq1bt+rdd9/V1q1btWjRIu3fv19NmzY1IWnGuV8fKFKkiCZOnKhdu3Zp3bp1ypcvn+rXr68zZ87YOGnGedBto8WLF2vjxo0KCQmxUTLbeZA2ePrpp1OMF5YtW2bDhBnvfm1w8OBBVatWTcWKFVN0dLR27Nihd999Vx4eHjZOmnHu1wb//Pzj4+M1c+ZMWSwWPfvsszZOmnHu1wZ9+/bVihUrNG/ePO3du1d9+/bV66+/ru+//97GSTPOvdrAMAw1b95chw4d0vfff69t27Ypb968qlu3bqbZl/Ig+4oy+xjxQdogs48R79cGmX2c+CB9ILOPER9mv/FjjRGNh1ShQgXj1VdfTTGtWLFixttvv/2wL2X3JBnfffed2TFMd/r0aUOSERMTY3YUU2XLls2YPn262TFs6vLly0bhwoWN1atXGzVq1DB69+5tdiSbGTJkiFGmTBmzY5hqwIABRrVq1cyO8UTp3bu3UbBgQSMpKcnsKDbzzDPPGF27dk0xrWXLlkaHDh1MSmR7V69eNZydnY0lS5akmF6mTBnjnXfeMSmV7fx7PJSUlGQEBQUZI0eOtE67fv26kTVrVuPzzz83IWHGu9eY8PDhw4YkY9u2bTbNZEsPMibetGmTIck4evSobULZ2IO0wcWLFw1Jxk8//WSbUDZ2tzY4fvy4kStXLmP37t1G3rx5jbFjx9o8m62k1QadO3c2mjVrZkoeM6TVBm3atHGoccGDfB80a9bMqF27tm0CmSCtNihRooTxwQcfpJgWERFhDB482IbJbOffbbBv3z5DkrF7927rtNu3bxv+/v7GtGnTTEiY8f69r8gRx4j32l/mCGNEw3iwfYaZeZz4IMuf2ceId2uDxx0jPtSRGjdv3tSWLVtUv379FNPr16+v9evXP3xFBZnCxYsXJUn+/v4mJzFHYmKiFixYoISEBFWuXNnsODbVo0cPPfPMM6pbt67ZUUxx4MABhYSEKH/+/Grbtq0OHTpkdiSb+uGHHxQZGalWrVopMDBQZcuW1bRp08yOZZqbN29q3rx56tq1qywWi9lxbKZatWr6+eeftX//fknSjh07tG7dOjVq1MjkZLZz+/ZtJSYmpvrFqaenp9atW2dSKvMcPnxYJ0+eTDFedHd3V40aNRgvOrCLFy/KYrE41NF8/3Tz5k1NnTpVWbNmVZkyZcyOYzNJSUnq2LGj+vfvrxIlSpgdxzTR0dEKDAxUkSJF9NJLL+n06dNmR7KZpKQkLV26VEWKFFGDBg0UGBioihUrOuRpnO84deqUli5dqm7dupkdxaaqVaumH374QX/++acMw1BUVJT279+vBg0amB3NJm7cuCFJKcaLzs7OcnNzy7TjxX/vK3LEMaKj7y+THqwNMvM48X7L7whjxLTaID3GiA9V1Dh79qwSExOVM2fOFNNz5sypkydPPlIA2DfDMNSvXz9Vq1ZNJUuWNDuOTe3atUtZsmSRu7u7Xn31VX333XcqXry42bFsZsGCBdq6datGjBhhdhRTVKxYUXPnztXKlSs1bdo0nTx5UlWqVNG5c+fMjmYzhw4d0uTJk1W4cGGtXLlSr776qnr16qW5c+eaHc0Uixcv1oULF9SlSxezo9jUgAED9Pzzz6tYsWJydXVV2bJl1adPHz3//PNmR7MZHx8fVa5cWcOGDdOJEyeUmJioefPmaePGjYqPjzc7ns3dGRMyXsQd169f19tvv6127drJ19fX7Dg2tWTJEmXJkkUeHh4aO3asVq9erYCAALNj2czHH38sFxcX9erVy+wopmnYsKG+/PJLrVmzRmPGjFFsbKxq165t3cGZ2Z0+fVpXrlzRyJEj9fTTT2vVqlVq0aKFWrZsqZiYGLPjmWLOnDny8fFRy5YtzY5iU+PHj1fx4sWVO3duubm56emnn9akSZNUrVo1s6PZRLFixZQ3b14NHDhQf/31l27evKmRI0fq5MmTmXK8mNa+IkcbIzry/rI7HqQNMvM48V7L7yhjxLu1QXqMEV0e5Un//gWqYRgO9atU/K1nz57auXNnpv1lwb0ULVpU27dv14ULF/Ttt9+qc+fOiomJcYjCxrFjx9S7d2+tWrUqU50L92E0bNjQertUqVKqXLmyChYsqDlz5qhfv34mJrOdpKQkRUZG6qOPPpIklS1bVr/99psmT56sTp06mZzO9mbMmKGGDRtmyvOF38vChQs1b948zZ8/XyVKlND27dvVp08fhYSEqHPnzmbHs5kvvvhCXbt2Va5cueTs7KyIiAi1a9dOW7duNTuaaRgvQkq+GGTbtm2VlJSkSZMmmR3H5mrVqqXt27fr7NmzmjZtmlq3bq2NGzcqMDDQ7GgZbsuWLfrss8+0detWh17327RpY71dsmRJRUZGKm/evFq6dKlD7NROSkqSJDVr1kx9+/aVJIWHh2v9+vX6/PPPVaNGDTPjmWLmzJlq3769w21HjR8/Xr/++qt++OEH5c2bV7/88ou6d++u4OBghzjy39XVVd9++626desmf39/OTs7q27duim2KzOTe+0rcpQxoiPvL7vjfm2Q2ceJ91p+RxkjptUG6TVGfKgjNQICAuTs7Jyqgnr69OlUlVZkfq+//rp++OEHRUVFKXfu3GbHsTk3NzcVKlRIkZGRGjFihMqUKaPPPvvM7Fg2sWXLFp0+fVrlypWTi4uLXFxcFBMTo/Hjx8vFxUWJiYlmR7Q5b29vlSpVSgcOHDA7is0EBwenKuKFhYUpLi7OpETmOXr0qH766Se9+OKLZkexuf79++vtt99W27ZtVapUKXXs2FF9+/Z1uKO4ChYsqJiYGF25ckXHjh3Tpk2bdOvWLeXPn9/saDYXFBQkSYwXoVu3bql169Y6fPiwVq9enel+ffcgvL29VahQIVWqVEkzZsyQi4uLZsyYYXYsm1i7dq1Onz6t0NBQ63jx6NGjeuONN5QvXz6z45kmODhYefPmdZgxY0BAgFxcXBgz/s/atWu1b98+hxszXrt2TYMGDdKnn36qJk2aqHTp0urZs6fatGmjTz75xOx4NlOuXDnrDyPj4+O1YsUKnTt3LtONF++2r8iRxoiOvr9Mun8bZPZx4v2W3xHGiHdrg/QaIz5UUcPNzU3lypXT6tWrU0xfvXq1qlSp8jAvBTtmGIZ69uypRYsWac2aNZnuH/CjMgzDYQ4jr1Onjnbt2qXt27db/yIjI9W+fXtt375dzs7OZke0uRs3bmjv3r0KDg42O4rNVK1aVfv27Usxbf/+/cqbN69Jicwza9YsBQYG6plnnjE7is1dvXpVTk4phxPOzs7WX2Y6Gm9vbwUHB+uvv/7SypUr1axZM7Mj2Vz+/PkVFBSUYrx48+ZNxcTEMF50IHc2VA8cOKCffvpJ2bNnNzvSE8GRxosdO3bUzp07U4wXQ0JC1L9/f61cudLseKY5d+6cjh075jBjRjc3N5UvX54x4//MmDFD5cqVy7TnTb+bW7du6datW4wZ/ydr1qzKkSOHDhw4oM2bN2ea8eL99hU5whiR/WUP1gaZeZz4qH0gM40R79cG6TVGfOjTT/Xr108dO3ZUZGSkKleurKlTpyouLk6vvvrqw76UXbpy5Yr++OMP6/3Dhw9r+/bt8vf3V2hoqInJbKdHjx6aP3++vv/+e/n4+Fir7FmzZpWnp6fJ6Wxj0KBBatiwofLkyaPLly9rwYIFio6O1ooVK8yOZhM+Pj6pzgfo7e2t7NmzO8y5It988001adJEoaGhOn36tIYPH65Lly451Ol2+vbtqypVquijjz5S69attWnTJk2dOlVTp041O5pNJSUladasWercubNcXB7prI52rUmTJvrwww8VGhqqEiVKaNu2bfr000/VtWtXs6PZ1MqVK2UYhooWLao//vhD/fv3V9GiRfXCCy+YHS1D3G881KdPH3300UcqXLiwChcurI8++kheXl5q166dianT1/3a4Pz584qLi9OJEyckybpDLygoyPpLRXt2r+UPCQnRc889p61bt2rJkiVKTEy0jhf9/f3l5uZmVux0da82yJ49uz788EM1bdpUwcHBOnfunCZNmqTjx4+rVatWJqZOX/dbD/69k8LV1VVBQUEqWrSoraNmmHu1gb+/v4YOHapnn31WwcHBOnLkiAYNGqSAgAC1aNHCxNTp6379oH///mrTpo2qV6+uWrVqacWKFfrxxx8VHR1tXuh09iD7CS5duqSvv/5aY8aMMStmhrpfG9SoUUP9+/eXp6en8ubNq5iYGM2dO1effvqpianT1/3a4Ouvv1aOHDkUGhqqXbt2qXfv3mrevHmKC2fbs/vtK7JYLJl+jPgg+8sy+xjxfm1w+/btTD1OvN/yJyQkZPox4v3aIHv27OkzRjQewX/+8x8jb968hpubmxEREWHExMQ8ysvYpaioKENSqr/OnTubHc1m0lp+ScasWbPMjmYzXbt2ta4DOXLkMOrUqWOsWrXK7FimqlGjhtG7d2+zY9hMmzZtjODgYMPV1dUICQkxWrZsafz2229mx7K5H3/80ShZsqTh7u5uFCtWzJg6darZkWxu5cqVhiRj3759ZkcxxaVLl4zevXsboaGhhoeHh1GgQAHjnXfeMW7cuGF2NJtauHChUaBAAcPNzc0ICgoyevToYVy4cMHsWBnmfuOhpKQkY8iQIUZQUJDh7u5uVK9e3di1a5e5odPZ/dpg1qxZaT4+ZMgQU3Onl3st/+HDh+86XoyKijI7erq5Vxtcu3bNaNGihRESEmK4ubkZwcHBRtOmTY1NmzaZHTtdPey2Ud68eY2xY8faNGNGu1cbXL161ahfv76RI0cOw9XV1QgNDTU6d+5sxMXFmR07XT1IP5gxY4ZRqFAhw8PDwyhTpoyxePFi8wJngAdpgylTphienp6ZdnxwvzaIj483unTpYoSEhBgeHh5G0aJFjTFjxhhJSUnmBk9H92uDzz77zMidO7f1+2Dw4MGZasz8IPuKMvsY8UHaILOPEe/XBpl9nHi/5XeEMeKj7Dd+lDGi5X9vBgAAAAAAAAAA8ER7qGtqAAAAAAAAAAAAmIWiBgAAAAAAAAAAsAsUNQAAAAAAAAAAgF2gqAEAAAAAAAAAAOwCRQ0AAAAAAAAAAGAXKGoAAAAAAAAAAAC7QFEDAAAAAAAAAADYBYoaAAAAAAAAAADALlDUAAAAAAAAAAAAdoGiBgAAAAAAAAAAsAsUNQAAAAAAAAAAgF2gqAEAAAAAAAAAAOwCRQ0AAAAAAAAAAGAXKGoAAAAAAAAAAAC7QFEDAAAAAAAAAADYBYoaAAAAAAAAAADALlDUAAAAAAAAAAAAdoGiBgAAAAAAAAAAsAsUNQAAAAAAAAAAgF2gqAEAAAAAAAAAAOwCRQ0AAAAAAAAAAGAXKGoAAAAAAAAAAAC7QFEDAAAAAAAAAADYBYoaAAAAAAAAAADALriYHQAAAKSfuLg4nT171uwYCggIUGhoqNkxHBJ9ANKT0Q/oA+mHz9MxPAmfs8RnDQAAnnwUNQAAyCTi4uIUFlZMV69eMzuKvLw8tXfv74+0U2T27Nl64YUXFBsbq8jIyFSPN27cWLt379aRI0fSIWnmEhcXp6JFi+r69etmR5GHh4f27dv3UH3AYrE80HxRUVGqWbPmIyZLf3f67OHDh5UvXz6z4yR/FxQL09VrV03N4eXppb2/733o74EWLVpoxYoVio+Pl5+fX5rztG/fXl999ZWGDx+ut99++6HbPl++fKpZs6Zmz54tSTpx4oSmTp2q5s2bKzw8/KHyZrTk7/YwXb1q8ufp5aW9ex/+8/zmm2/UqlUrLViwQG3atEnxWJkyZbRz506tWLFCDRo0SPFYwYIFlTVrVm3duvWxs1ssFvXo0UMTJ0587NfKKE/K5yw9+mctSRs3btTIkSO1ZcsWnTp1Sn5+fipQoICqVKmiMWPGSEpe/0qWLKklS5akd3QAAOAgKGoAAJBJnD17VlevXtO8zxsorIi/aTn27j+vDq+u1NmzZ/mlp42dPXtW169fV9myZZUlSxbTcly5ckXbtm176D6wYcOGFPeHDRumqKgorVmzJsX04sWLp0vOzOrs2bO6eu2qBrYcqNAAc9bBuLNxGrFoxCN9D3Tr1k2LFy/W/Pnz1b1791SPX7x4Ud99950aN26srl27qkaNGgoODn6o9/juu+/k6+trvX/ixAm9//77ypcv3xNX1Ej+br+q0aNHq2DBgqZkOHjwoPr37/9In2fNmjVlsVgUFRWVoqhx/vx57dq1S97e3oqKikpR1Dh+/LgOHTqkfv36pdsyPOnufM5zh01QsfyFTcvx++ED6vTu64/0WS9dulRNmzZVzZo1NWrUKAUHBys+Pl6bN2/WggULrEUNAACAx0VRAwCATCasiL8iygSaHQMmypIly11/4f4kq1SpUor7OXLkkJOTU6rp/3b16lV5eXllZDS7FBoQqsIhRcyO8dAaNmyokJAQzZw5M82ixn//+19du3ZN3bp1U44cOZQjR46Hfo+yZcumR1SbKliwoEqUKGF2jIcWEBCgkiVLKjo6OsX0mJgYubi4qFu3boqKikrx2J37tWrVeqz3vnbtmjw9PR/rNWytWP7CiggrZXaMRzJq1Cjlz59fK1eulIvL37sa2rZtq1GjRpmYDAAAZDZcKBwAANi169eva+DAgcqfP7/c3NyUK1cu9ejRQxcuXEgxX758+dS4cWMtWbJEZcuWlaenp8LCwqynv5g9e7bCwsLk7e2tChUqaPPmzanea/PmzWratKn8/f3l4eGhsmXL6quvvrLFYuJ/atasqZIlS+qXX35RlSpV5OXlpa5du0qSFi5cqPr16ys4ONj6+b799ttKSEiwPn/cuHGyWCz6448/Ur32gAED5ObmluKc9j/99JPq1KkjX19feXl5qWrVqvr5558zfkEdmLOzszp37qwtW7Zo165dqR6fNWuWgoOD1bBhQ82ePVsWiyXF6ei2bdumxo0bKzAwUO7u7goJCdEzzzyj48ePW+fJly+funTpIkmKjo5W+fLlJUkvvPCCLBaLLBaLhg4dmpGL6VBq1aqlffv2KT4+3jrtTrs3atRIW7Zs0eXLl1M85uzsrKeeeuqhv+MXLVqksmXLysPDQ++//36aeQzD0KBBg+Tq6qpp06ZlyDI7onPnzikgICBFQeMOJ6fUux5WrFihiIgIeXp6qlixYpo5c2aqeXbv3q1mzZopW7Zs8vDwUHh4uObMmZNinujoaFksFs2bN0/9+vVTUFCQPD09VaNGDW3bti39FhAAADwxKGoAAIAnUmJiom7fvp3qzzAM6zyGYah58+b65JNP1LFjRy1dulT9+vXTnDlzVLt2bd24cSPFa+7YsUMDBw7UgAEDtGjRImXNmlUtW7bUkCFDNH36dH300Uf68ssvdfHiRTVu3FjXrv19fZKoqChVrVpVFy5c0Oeff67vv/9e4eHhatOmjfW8/LCN+Ph4dejQQe3atdOyZcusv+Y/cOCAGjVqpBkzZmjFihXq06ePvvrqKzVp0sT63A4dOsjNzS3VZ5aYmKh58+apSZMmCggIkCTNmzdP9evXl6+vr+bMmaOvvvpK/v7+atCgAYWNDNa1a1dZLJZUOzn37NmjTZs2qXPnznJ2dk71vISEBNWrV0+nTp3Sf/7zH61evVrjxo1TaGhoip3m/xQREaFZs2ZJkgYPHqwNGzZow4YNevHFF9N/wRzUnSMu/nm0RlRUlGrUqKGqVavKYrFo7dq1KR6LiIiQr6/vQ33Hb926Vf3791evXr20YsUKPfvss6my3LhxQ+3atdPEiRP1448/6qWXXsqYhXZAlStX1saNG9WrVy9t3LhRt27duuu8O3bs0BtvvKG+ffvq+++/V+nSpdWtWzf98ssv1nn27dunKlWq6LffftP48eO1aNEiFS9eXF26dEnzyI9Bgwbp0KFDmj59uqZPn64TJ06oZs2aOnToUIYsLwAAMA+nnwIAAE+ke51yKG/evJKkVatWaeXKlRo1apT69+8vSapXr57y5MmjNm3aaO7cuSl2WJ07d06//vqrcuXKJUkKCQlReHi4pk2bpj/++MN6CiOLxaLmzZvrp59+su4Q7969u0qUKKE1a9ZYf4XaoEEDnT17VoMGDVKnTp3S/CUq0t/58+f19ddfq3bt2immDx482HrbMAxVrVpVYWFhqlGjhnbu3KnSpUsrICBAjRs31pw5c/TBBx9YP7NVq1bpxIkTeuGFFyQln9Kqd+/eaty4sb777jvr6zZq1EgREREaNGiQNm7caIOldUyFChVS9erVNW/ePI0aNUqurq6SZC1y3Dk6599+//13nTt3TjNmzFCzZs2s01u3bn3X9/L19VXJkiUlJZ/i6X6nO8PDq1GjhpycnBQdHa3nn39e586d0+7duzV69GhlyZJFERERioqKUqNGjXTs2DEdPnxYrVq1eujv+NOnT2vPnj0qUiTt066dP39ezZo10+HDh7V27VqVKVPGJsvvKEaOHKnff/9dEyZM0IQJE+Tq6qry5curSZMm6tmzZ4prPZ09e1b/93//Z71uR/Xq1fXzzz9r/vz5ql69uiRp6NChunnzpqKiopQnTx5Jyd/BFy5c0Pvvv69XXnlFWbNmtb5mjhw59N1338lisUiSqlWrpsKFC2vEiBEckQMAQCbDljcAAHgizZ07V7Gxsan+qlWrZp3nzgWk75xG5o5WrVrJ29s71a/pw8PDrQUNSQoLC5OUfEqjf16T4c70o0ePSpL++OMP/f7772rfvr0kpThypFGjRoqPj9e+ffvSaclxP9myZUtV0JCkQ4cOqV27dgoKCpKzs7NcXV1Vo0YNSdLevXut873wwgs6fvy4fvrpJ+u0WbNmKSgoSA0bNpQkrV+/XufPn1fnzp1TfN5JSUl6+umnFRsbm+K0Vkh/3bp109mzZ/XDDz9ISl7v5s2bp6eeekqFC6d9IeVChQopW7ZsGjBggD7//HPt2bPHlpFxF9myZVOZMmWsR2rExMTI2dlZVatWlZRc9LhzHY1/Xk/jYb/jS5cufdeCxuHDh1W5cmVdunRJv/76KwWNDJA9e3atXbtWsbGxGjlypJo1a6b9+/dr4MCBKlWqVIpT+4WHh6e4ELmHh4eKFCli/b8rJf+Pr1OnjrWgcUeXLl109epVbdiwIcX0du3aWQsaUvIPIKpUqZLqmi0AAMD+UdQAAABPpLCwMEVGRqb6++evMs+dOycXF5dUFwq2WCwKCgrSuXPnUkz39/dPcd/Nze2e069fvy5JOnXqlCTpzTfflKura4q/O6c++ufOGmSs4ODgVNOuXLmip556Shs3btTw4cMVHR2t2NhYLVq0SJJSnEqsYcOGCg4Otp5y6K+//tIPP/ygTp06WU9pdOczf+6551J95h9//LEMw9D58+czelEd2nPPPaesWbNaP6dly5bp1KlT6tat212fkzVrVsXExCg8PFyDBg1SiRIlFBISoiFDhtzzVDjIeLVq1dL+/ft14sQJRUVFqVy5ctZf7t+59sHFixcVFRUlFxcXVatW7aG/49P6brhj06ZN2r9/v9q0aaPcuXOn/wLCKjIyUgMGDNDXX3+tEydOqG/fvjpy5EiKU0Zlz5491fPc3d1TfFefO3cuzc80JCTE+vg/BQUFpZo3rX4CAADsH6efAgAAdit79uy6ffu2zpw5k2Knl2EYOnnypPXiv4/rzjUWBg4cqJYtW6Y5T9GiRdPlvXB///wl7h1r1qzRiRMnFB0dbT06Q1KqiwlLyRei7tixo8aPH68LFy5o/vz5unHjhvXUU9Lfn/mECRPuejqinDlzPuaS4F48PT31/PPPa9q0aYqPj9fMmTPl4+OjVq1a3fN5pUqV0oIFC2QYhnbu3KnZs2frgw8+kKenp95++20bpce/1apVS59++qmio6MVHR2tRo0aWR+7cwTeL7/8Yr2AeJYsWR76Oz6t74Y72rRpo6CgIL3zzjtKSkpKcbo6ZBxXV1cNGTJEY8eO1e7dux/qudmzZ09xcfk7Tpw4Ienv7+k7Tp48mWrekydPpllAAQAA9o0jNQAAgN2qU6eOpOQLOv/Tt99+q4SEBOvjj6to0aIqXLiwduzYkebRI5GRkfLx8UmX98KjubMz093dPcX0KVOmpDn/Cy+8oOvXr+u///2vZs+ercqVK6tYsWLWx6tWrSo/Pz/t2bPnrp/5nSN6kHG6deumxMREjR49WsuWLVPbtm1TnCruXiwWi8qUKaOxY8fKz89PW7duveu8d/rNP38ljvRVvXp1OTs765tvvtFvv/2mmjVrWh/LmjWrwsPDNWfOHB05csR6YfH0/o4fPHiwxo0bp/fee08DBw58vAVCKmkVIKS/T/935wiLB1WnTh1rwfqf5s6dKy8vr1QF5//+978yDMN6/+jRo1q/fn2KvgYAADIHjtQAAAB2q169emrQoIEGDBigS5cuqWrVqtq5c6eGDBmismXLqmPHjun2XlOmTFHDhg3VoEEDdenSRbly5dL58+e1d+9ebd26VV9//XW6vRceXpUqVZQtWza9+uqrGjJkiFxdXfXll19qx44dac5frFgxVa5cWSNGjNCxY8c0derUFI9nyZJFEyZMUOfOnXX+/Hk999xzCgwM1JkzZ7Rjxw6dOXNGkydPtsWiObTIyEiVLl1a48aNk2EY9zz1lCQtWbJEkyZNUvPmzVWgQAEZhqFFixbpwoULqlev3l2fV7BgQXl6eurLL79UWFiYsmTJopCQkIfeCYu78/X1VUREhBYvXiwnJyfr9TTuqFGjhsaNGydJ1qJGRnzH9+7dW1myZNHLL7+sK1euaPz48fc8wgMPrkGDBsqdO7eaNGmiYsWKKSkpSdu3b9eYMWOUJUsW9e7d+6Feb8iQIVqyZIlq1aql9957T/7+/vryyy+1dOlSjRo1KsXpKKXkC8W3aNFCL730ki5evKghQ4bIw8ODAhYAAJkQRQ0AADKZvfvNPc+/Ld/fYrFo8eLFGjp0qGbNmqUPP/xQAQEB6tixoz766KNUv9p/HLVq1dKmTZv04Ycfqk+fPvrrr7+UPXt2FS9eXK1bt06390kPV65ccbj3z549u5YuXao33nhDHTp0kLe3t5o1a6aFCxcqIiIizee88MILevnll+Xp6ak2bdqkerxDhw4KDQ3VqFGj9Morr+jy5csKDAxUeHh4qgsXP4nizsZlivfu1q2bevfureLFi6tixYr3nLdw4cLy8/PTqFGjdOLECbm5ualo0aKaPXu2OnfufNfneXl5aebMmXr//fdVv3593bp1S0OGDNHQoUPTbTke18GDB+3+vWvVqqXY2FiVLVtWvr6+KR6rUaOGxo4dKzc3N1WpUkVSxn3Hd+vWTd7e3urYsaMSEhI0ffp0OTk9GScx+P3wAbt9/8GDB+v777/X2LFjFR8frxs3big4OFh169bVwIEDFRYW9lCvV7RoUa1fv16DBg1Sjx49dO3aNYWFhWnWrFlpfgd/9NFHio2N1QsvvKBLly6pQoUKWrBggQoWLPjIywQAAJ5MFuOfx2cCAAC7FRcXp7CwYrp61fzTp3h5eWrv3t8VGhpqdhSHEhcXp6JFi1ovcG4mDw8P7du3jz5ggri4OIUVC9PVa1dNzeHl6aW9v++lDzym5O/2MF29avLn6eWlvXv5PDPKk/I5S/b3WUdHR6tWrVr6+uuv9dxzz5kdBwAA2ABHagAAkEmEhoZq797fdfbsWbOjKCAgwG52hmQmoaGh2rdvH33AwYWGhmrv73tN7wf0gfSR/N3O55nZPSmfs8RnDQAAnnwcqQEAAAAAAOwSR2oAAOB4KGoAAAAAAAAAAAC78GRcDQ0AAAAAAAAAAOA+KGoAAAAAAAAAAAC7QFEDAAAAAAAAAADYBYoaAAAAAAAAAADALlDUAAAAAAAAAAAAdoGiBgAAAAAAAAAAsAsUNQAAAAAAAAAAgF2gqAEAAAAAAAAAAOwCRQ0AAAAAAAAAAGAXKGoAAAAAAAAAAAC7QFEDAAAAAAAAAADYBYoaAAAAAAAAAADALlDUAAAAAAAAAAAAdoGiBgAAAAAAAAAAsAsUNQAAAAAAAAAAgF2gqAEAAAAAAAAAAOwCRQ0AAAAAAAAAAGAX/h/fTg28xz7D0QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1600x175 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plan = Plan(home_area=\"h\", home_location=Location(loc=Point(0, 0), area=\"h\"))\n",
    "plan.add(\n",
    "    Activity(\n",
    "        seq=1, act=\"home\", area=\"h\", loc=Point(0, 0), start_time=mtdt(0), end_time=mtdt(5 * 60)\n",
    "    )\n",
    ")\n",
    "plan.add(\n",
    "    Leg(\n",
    "        seq=1,\n",
    "        mode=\"car\",\n",
    "        start_area=\"h\",\n",
    "        end_area=\"w\",\n",
    "        start_loc=Point(0, 0),\n",
    "        end_loc=Point(1, 2),\n",
    "        start_time=mtdt(5 * 60),\n",
    "        end_time=mtdt(6 * 60),\n",
    "    )\n",
    ")\n",
    "plan.add(\n",
    "    Activity(\n",
    "        seq=2,\n",
    "        act=\"visit\",\n",
    "        area=\"w\",\n",
    "        loc=Point(1, 2),\n",
    "        start_time=mtdt(6 * 60),\n",
    "        end_time=mtdt(17 * 60),\n",
    "    )\n",
    ")\n",
    "plan.add(\n",
    "    Leg(\n",
    "        seq=2,\n",
    "        mode=\"car\",\n",
    "        start_area=\"w\",\n",
    "        end_area=\"w\",\n",
    "        start_loc=Point(1, 2),\n",
    "        end_loc=Point(0, 0),\n",
    "        start_time=mtdt(17 * 60),\n",
    "        end_time=mtdt(18 * 60),\n",
    "    )\n",
    ")\n",
    "plan.add(\n",
    "    Activity(\n",
    "        seq=3,\n",
    "        act=\"work\",\n",
    "        area=\"w\",\n",
    "        loc=Point(0, 0),\n",
    "        start_time=mtdt(18 * 60),\n",
    "        end_time=mtdt(19 * 60),\n",
    "    )\n",
    ")\n",
    "plan.add(\n",
    "    Leg(\n",
    "        seq=3,\n",
    "        mode=\"car\",\n",
    "        start_area=\"w\",\n",
    "        end_area=\"na\",\n",
    "        start_loc=Point(0, 0),\n",
    "        end_loc=Point(1, 2),\n",
    "        start_time=mtdt(19 * 60),\n",
    "        end_time=mtdt(20 * 60),\n",
    "    )\n",
    ")\n",
    "plan.add(\n",
    "    Activity(\n",
    "        seq=4,\n",
    "        act=\"shop\",\n",
    "        area=\"na\",\n",
    "        loc=Point(2, 2),\n",
    "        start_time=mtdt(20 * 60),\n",
    "        end_time=mtdt(22 * 60),\n",
    "    )\n",
    ")\n",
    "plan.add(\n",
    "    Leg(\n",
    "        seq=4,\n",
    "        mode=\"car\",\n",
    "        start_area=\"na\",\n",
    "        end_area=\"h\",\n",
    "        start_loc=Point(2, 2),\n",
    "        end_loc=Point(0, 0),\n",
    "        start_time=mtdt(22 * 60),\n",
    "        end_time=mtdt(23 * 60),\n",
    "    )\n",
    ")\n",
    "plan.add(\n",
    "    Activity(\n",
    "        seq=5, act=\"home\", area=\"h\", loc=Point(0, 0), start_time=mtdt(23 * 60), end_time=END_OF_DAY\n",
    "    )\n",
    ")\n",
    "plan.plot()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1f6c39ea",
   "metadata": {},
   "source": [
    "The `od` object holds origin-destination data, for example travel time and travel distance between each origin and destination, for each travel mode.\n",
    "\n",
    "For the discretionary trips location choice, the `od` object needs to include the variable `od_probs`, which provides provides origin-destination probabilities to be considered -amongst others- as part of the location choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "15027351",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.693275Z",
     "iopub.status.busy": "2024-04-05T15:22:28.693091Z",
     "iopub.status.idle": "2024-04-05T15:22:28.700460Z",
     "shell.execute_reply": "2024-04-05T15:22:28.700132Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Origin-destination dataset \n",
       "--------------------------------------------------\n",
       "Labels(vars=['time', 'distance', 'od_probs'], origin_zones=('h', 'b', 'w'), destination_zones=('h', 'b', 'w'), mode=['car', 'bus'])\n",
       "--------------------------------------------------\n",
       "time - car:\n",
       "[[20. 30. 40.]\n",
       " [30. 10. 30.]\n",
       " [40. 30. 20.]]\n",
       "--------------------------------------------------\n",
       "time - bus:\n",
       "[[30. 40. 45.]\n",
       " [40. 10. 40.]\n",
       " [45. 40. 30.]]\n",
       "--------------------------------------------------\n",
       "distance - car:\n",
       "[[5. 6. 8.]\n",
       " [6. 2. 6.]\n",
       " [8. 6. 5.]]\n",
       "--------------------------------------------------\n",
       "distance - bus:\n",
       "[[5. 7. 9.]\n",
       " [7. 2. 7.]\n",
       " [9. 7. 5.]]\n",
       "--------------------------------------------------\n",
       "od_probs - car:\n",
       "[[ 500.  333.  250.]\n",
       " [ 333. 1000.  333.]\n",
       " [ 250.  333.  500.]]\n",
       "--------------------------------------------------\n",
       "od_probs - bus:\n",
       "[[ 333.  250.  222.]\n",
       " [ 250. 1000.  250.]\n",
       " [ 222.  250.  333.]]\n",
       "--------------------------------------------------"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zone_labels = (\"h\", \"b\", \"w\")\n",
    "od = ODFactory.from_matrices(\n",
    "    [\n",
    "        ODMatrix(\n",
    "            \"time\",\n",
    "            \"car\",\n",
    "            zone_labels,\n",
    "            zone_labels,\n",
    "            np.array([[20, 30, 40], [30, 10, 30], [40, 30, 20]]),\n",
    "        ),\n",
    "        ODMatrix(\n",
    "            \"time\",\n",
    "            \"bus\",\n",
    "            zone_labels,\n",
    "            zone_labels,\n",
    "            np.array([[30, 40, 45], [40, 10, 40], [45, 40, 30]]),\n",
    "        ),\n",
    "        ODMatrix(\n",
    "            \"distance\", \"car\", zone_labels, zone_labels, np.array([[5, 6, 8], [6, 2, 6], [8, 6, 5]])\n",
    "        ),\n",
    "        ODMatrix(\n",
    "            \"distance\", \"bus\", zone_labels, zone_labels, np.array([[5, 7, 9], [7, 2, 7], [9, 7, 5]])\n",
    "        ),\n",
    "        ODMatrix(\n",
    "            \"od_probs\",\n",
    "            \"car\",\n",
    "            zone_labels,\n",
    "            zone_labels,\n",
    "            np.array([[500, 333, 250], [333, 1000, 333], [250, 333, 500]]),\n",
    "        ),\n",
    "        ODMatrix(\n",
    "            \"od_probs\",\n",
    "            \"bus\",\n",
    "            zone_labels,\n",
    "            zone_labels,\n",
    "            np.array([[333, 250, 222], [250, 1000, 250], [222, 250, 333]]),\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "od"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d679888-5933-4ac3-9dba-babdb324e67c",
   "metadata": {},
   "source": [
    "## Model application"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08ebbb60-8756-4789-8ed4-11c12feaf9fc",
   "metadata": {},
   "source": [
    "Initially, the `shop` activity was allocated to zone `na`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e1163184",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.703428Z",
     "iopub.status.busy": "2024-04-05T15:22:28.703234Z",
     "iopub.status.idle": "2024-04-05T15:22:28.709407Z",
     "shell.execute_reply": "2024-04-05T15:22:28.709060Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Work locations and travel modes:\n",
      "+-----+---------+----------+\n",
      "| seq | purpose | location |\n",
      "+-----+---------+----------+\n",
      "|  0  |   home  |    h     |\n",
      "|  1  |  visit  |    w     |\n",
      "|  2  |   work  |    w     |\n",
      "|  3  |   shop  |    na    |\n",
      "|  4  |   home  |    h     |\n",
      "+-----+---------+----------+\n"
     ]
    }
   ],
   "source": [
    "def print_activity_locs(plan):\n",
    "    summary = PrettyTable([\"seq\", \"purpose\", \"location\"])\n",
    "    for seq, act in enumerate(plan.activities):\n",
    "        summary.add_row([seq, act.act, act.location.area])\n",
    "    print(summary)\n",
    "\n",
    "\n",
    "print(\"Work locations and travel modes:\")\n",
    "print_activity_locs(plan)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8af752f7-2d40-4a74-bc0e-a400dd9ed5a5",
   "metadata": {},
   "source": [
    "We apply the discretionary location choice model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0c8f720d-d4f1-4a54-8222-bd482fa6a2b4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.711268Z",
     "iopub.status.busy": "2024-04-05T15:22:28.711107Z",
     "iopub.status.idle": "2024-04-05T15:22:28.714717Z",
     "shell.execute_reply": "2024-04-05T15:22:28.714354Z"
    }
   },
   "outputs": [],
   "source": [
    "plan_new = deepcopy(plan)\n",
    "planner = DiscretionaryTrips(plan=plan_new, od=od)\n",
    "planner.update_plan()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31387271-af1b-4dd1-81f3-e014a6a83106",
   "metadata": {},
   "source": [
    "The location of the shop trip has now been updated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "67375c1e-9164-4341-b3e2-491058f05283",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.716764Z",
     "iopub.status.busy": "2024-04-05T15:22:28.716600Z",
     "iopub.status.idle": "2024-04-05T15:22:28.719623Z",
     "shell.execute_reply": "2024-04-05T15:22:28.719177Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----+---------+----------+\n",
      "| seq | purpose | location |\n",
      "+-----+---------+----------+\n",
      "|  0  |   home  |    h     |\n",
      "|  1  |  visit  |    w     |\n",
      "|  2  |   work  |    w     |\n",
      "|  3  |   shop  |    w     |\n",
      "|  4  |   home  |    h     |\n",
      "+-----+---------+----------+\n"
     ]
    }
   ],
   "source": [
    "print_activity_locs(plan_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f393b4d4-71de-4d6a-9244-c6d9097db7b0",
   "metadata": {},
   "source": [
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "960e5d3f-3786-4e60-bf4e-5f182293981e",
   "metadata": {},
   "source": [
    "## Breaking down the calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8bb1ad7-f073-4c34-92cd-77c71cae080a",
   "metadata": {},
   "source": [
    "This section goes through some of the calculations behind the model application."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "134e97af-c173-46c1-b738-d08a5008e2a0",
   "metadata": {},
   "source": [
    "We select the work->shop->home chain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6703633e-4912-4b53-badc-1df09d07c1a9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.721680Z",
     "iopub.status.busy": "2024-04-05T15:22:28.721517Z",
     "iopub.status.idle": "2024-04-05T15:22:28.724732Z",
     "shell.execute_reply": "2024-04-05T15:22:28.724323Z"
    }
   },
   "outputs": [],
   "source": [
    "plan_new = deepcopy(plan)\n",
    "chains = get_trip_chains_either_anchor(plan_new)\n",
    "chain = chains[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7c50299c-d018-4a95-a345-4d60275af486",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.727043Z",
     "iopub.status.busy": "2024-04-05T15:22:28.726840Z",
     "iopub.status.idle": "2024-04-05T15:22:28.730287Z",
     "shell.execute_reply": "2024-04-05T15:22:28.729941Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['work', 'travel', 'shop', 'travel', 'home']"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[act.act for act in chain]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8295978c-7a9e-4f94-90ea-967822e22b0c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.732096Z",
     "iopub.status.busy": "2024-04-05T15:22:28.731939Z",
     "iopub.status.idle": "2024-04-05T15:22:28.734503Z",
     "shell.execute_reply": "2024-04-05T15:22:28.734160Z"
    }
   },
   "outputs": [],
   "source": [
    "model = DiscretionaryTripOD(chain, od)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e27e9d6-bc1b-43d4-bade-82659d557d69",
   "metadata": {},
   "source": [
    "The durations of the two trips are one hour each:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "652ff1ad-8371-48fb-ae57-4f061b33f6e4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.736361Z",
     "iopub.status.busy": "2024-04-05T15:22:28.736117Z",
     "iopub.status.idle": "2024-04-05T15:22:28.739330Z",
     "shell.execute_reply": "2024-04-05T15:22:28.738947Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[datetime.timedelta(seconds=3600), datetime.timedelta(seconds=3600)]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[x.duration for x in chain if isinstance(x, Leg)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2034f85c-aef6-4db4-a0b4-a456e19b3e08",
   "metadata": {},
   "source": [
    "... and therefore, the observed leg ratio is 0.5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "285a6048-e7ed-408b-86e7-09770c3ae3c4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.741226Z",
     "iopub.status.busy": "2024-04-05T15:22:28.741072Z",
     "iopub.status.idle": "2024-04-05T15:22:28.748128Z",
     "shell.execute_reply": "2024-04-05T15:22:28.747802Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.observed_leg_ratio"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2cacab7-ed83-4050-8fbd-5f93784160fb",
   "metadata": {},
   "source": [
    "Of the possible destination zones (h/b/w), zone b has the closest-to-observed leg ratio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a960cc2f-ac94-4c1e-85e3-bb9a5543ee4c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.750145Z",
     "iopub.status.busy": "2024-04-05T15:22:28.749975Z",
     "iopub.status.idle": "2024-04-05T15:22:28.753481Z",
     "shell.execute_reply": "2024-04-05T15:22:28.752965Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.66666667, 0.5       , 0.33333333])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.leg_ratios"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1adfdcba-b973-4f5b-a597-459cebd44692",
   "metadata": {},
   "source": [
    "... and therefore, offers the highest \"leg_ratio\" probability:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "639c6ee3-de4f-4cb9-8cc0-49de505aa3bc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.755408Z",
     "iopub.status.busy": "2024-04-05T15:22:28.755257Z",
     "iopub.status.idle": "2024-04-05T15:22:28.758608Z",
     "shell.execute_reply": "2024-04-05T15:22:28.758288Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.66666667, 1.        , 0.66666667])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.leg_ratio_p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d30aeda-9925-4c7f-b4b9-90cdcf2a7198",
   "metadata": {},
   "source": [
    "The diversion factors are the same each destination:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "27449617-345d-4ded-b425-30684326347f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.760614Z",
     "iopub.status.busy": "2024-04-05T15:22:28.760443Z",
     "iopub.status.idle": "2024-04-05T15:22:28.763832Z",
     "shell.execute_reply": "2024-04-05T15:22:28.763467Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.5, 1.5, 1.5])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.diversion_factors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e215364-9a93-4b84-9154-e1238c2f9d19",
   "metadata": {},
   "source": [
    "... and therefore the diversion probabilities are equal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5a2d3f46-a141-4d89-a06c-600640523b91",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.765775Z",
     "iopub.status.busy": "2024-04-05T15:22:28.765619Z",
     "iopub.status.idle": "2024-04-05T15:22:28.768662Z",
     "shell.execute_reply": "2024-04-05T15:22:28.768367Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.54545455, 0.54545455, 0.54545455])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.diversion_p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e0b6ac6-64d5-4c0e-a8be-5236b9ced6f6",
   "metadata": {},
   "source": [
    "Attraction probabilities reflect the data in the `od` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b3cdb221-80b3-45ce-a5bc-e0e1255bc2da",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.770452Z",
     "iopub.status.busy": "2024-04-05T15:22:28.770311Z",
     "iopub.status.idle": "2024-04-05T15:22:28.774299Z",
     "shell.execute_reply": "2024-04-05T15:22:28.773958Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([250., 333., 500.])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.od[\"od_probs\", \"w\", :, \"car\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5156166-446c-42b1-ad06-0f4a9425ca54",
   "metadata": {},
   "source": [
    ".. the more opportunities it offers, the higher the probability of a destination:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "94d6b7be-8278-43a5-962f-cd9dd44d0725",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.776411Z",
     "iopub.status.busy": "2024-04-05T15:22:28.776245Z",
     "iopub.status.idle": "2024-04-05T15:22:28.779588Z",
     "shell.execute_reply": "2024-04-05T15:22:28.779155Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.23084026, 0.30747922, 0.46168052])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.attraction_p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ccdf3d1-0946-489d-95a3-cc0647d60a7d",
   "metadata": {},
   "source": [
    "The final probabilities are defined as the product of the three elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5fb51c28-7109-4a43-83da-e7e9865f25a6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.781661Z",
     "iopub.status.busy": "2024-04-05T15:22:28.781488Z",
     "iopub.status.idle": "2024-04-05T15:22:28.784974Z",
     "shell.execute_reply": "2024-04-05T15:22:28.784566Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.20008003, 0.3997599 , 0.40016006])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probs = model.leg_ratio_p * model.diversion_p * model.attraction_p\n",
    "probs /= probs.sum()\n",
    "probs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d37ac59a-056e-442c-8983-68228fdaefb5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.786900Z",
     "iopub.status.busy": "2024-04-05T15:22:28.786744Z",
     "iopub.status.idle": "2024-04-05T15:22:28.792874Z",
     "shell.execute_reply": "2024-04-05T15:22:28.792560Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.20008003, 0.3997599 , 0.40016006])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.destination_p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8aba630a-281d-4f24-b3d6-e1b843344380",
   "metadata": {},
   "source": [
    "Running the model multiple times gives us the same probabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "074ab041-2c0a-40bb-bb21-12c1ff83516c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-04-05T15:22:28.794810Z",
     "iopub.status.busy": "2024-04-05T15:22:28.794656Z",
     "iopub.status.idle": "2024-04-05T15:22:28.882662Z",
     "shell.execute_reply": "2024-04-05T15:22:28.882270Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "h    0.200\n",
       "b    0.403\n",
       "w    0.397\n",
       "Name: proportion, dtype: float64"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = []\n",
    "for i in range(10**3):\n",
    "    planner = DiscretionaryTrips(plan=plan_new, od=od)\n",
    "    planner.update_plan()\n",
    "    res.append(list(plan_new.activities)[3].location.area)\n",
    "res = pd.Series(res)\n",
    "res.value_counts(normalize=True).loc[list(od.labels.destination_zones)]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}