{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Route Choice with sub-area analysis\n\nIn this example, we show how to perform sub-area analysis using route choice assignment, \nfor a city in La Serena Metropolitan Area in Chile.\n\n.. admonition:: References\n \n   * :doc:`../../route_choice`\n\n.. seealso::\n    Several functions, methods, classes and modules are used in this example:\n\n    * :func:`aequilibrae.paths.Graph`\n    * :func:`aequilibrae.paths.RouteChoice`\n    * :func:`aequilibrae.paths.SubAreaAnalysis`\n    * :func:`aequilibrae.matrix.AequilibraeMatrix`\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Imports\nfrom uuid import uuid4\nfrom tempfile import gettempdir\nfrom os.path import join\nimport itertools\n\nimport pandas as pd\nimport numpy as np\nimport folium\n\nfrom aequilibrae.utils.create_example import create_example"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# We create the example project inside our temp folder\nfldr = join(gettempdir(), uuid4().hex)\n\nproject = create_example(fldr, \"coquimbo\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import logging\nimport sys"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# We the project opens, we can tell the logger to direct all messages to the terminal as well\nlogger = project.logger\nstdout_handler = logging.StreamHandler(sys.stdout)\nformatter = logging.Formatter(\"%(asctime)s;%(levelname)s ; %(message)s\")\nstdout_handler.setFormatter(formatter)\nlogger.addHandler(stdout_handler)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Model parameters\nWe'll set the parameters for our route choice model. These are the parameters that will be\nused to calculate the utility of each path. In our example, the utility is equal to\n$distance * theta$, and the path overlap factor (PSL) is equal to $beta$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "theta = 0.011  # Distance factor\n\nbeta = 1.1  # PSL parameter"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's build all graphs\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "project.network.build_graphs()\n# We get warnings that several fields in the project are filled with NaNs.\n# This is true, but we won't use those fields."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We grab the graph for cars\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "graph = project.network.graphs[\"c\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We also see what graphs are available\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "project.network.graphs.keys()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's say that utility is just a function of distance.\nSo we build our *utility* field as the $distance * theta$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "graph.network = graph.network.assign(utility=graph.network.distance * theta)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Prepare the graph with all nodes of interest as centroids\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "graph.prepare_graph(graph.centroids)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And set the cost of the graph the as the utility field just created\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "graph.set_graph(\"utility\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Mock demand matrix\n\nWe'll create a mock demand matrix with demand `10` for every zone and prepare it for computation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from aequilibrae.matrix import AequilibraeMatrix\n\nnames_list = [\"demand\"]\n\nmat = AequilibraeMatrix()\nmat.create_empty(zones=graph.num_zones, matrix_names=names_list, memory_only=True)\nmat.index = graph.centroids[:]\nmat.matrices[:, :, 0] = np.full((graph.num_zones, graph.num_zones), 10.0)\nmat.computational_view()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Sub-area preparation\n\nWe need to define some polygon for out sub-area analysis, here we'll use a section of zones and\ncreate out polygon as the union of their geometry. It's best to choose a polygon that avoids\nany unnecessary intersections with links as the resource requirements of this approach grow\nquadratically with the number of links cut.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "zones_of_interest = [29, 30, 31, 32, 33, 34, 37, 38, 39, 40, 49, 50, 51, 52, 57, 58, 59, 60]\nzones = project.zoning.data.set_index(\"zone_id\")\nzones = zones.loc[zones_of_interest]\nzones.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Sub-area analysis\n\nFrom here there are two main paths to conduct a sub-area analysis, manual or automated.\nAequilibraE ships with a small class that handle most of the details regarding the implementation\nand extract of the relevant data. It also exposes all the tools necessary to conduct this analysis\nyourself if you need fine grained control.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Automated sub-area analysis\n\nWe first construct out ``SubAreaAnalysis`` object from the graph, zones, and matrix we previously constructed, then\nconfigure the route choice assignment and execute it. From there the ``post_process`` method is able to use the route\nchoice assignment results to construct the desired demand matrix as a DataFrame. If we were interested in the original\norigin and destination IDs for each entry we could use `subarea.post_process(keep_original_ods=True)` instead. This\nwill attach the true ODs from the select link OD matrix as part of the index. However, this will create a\nsignificantly larger, but more flexible matrix.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from aequilibrae.paths import SubAreaAnalysis\n\nsubarea = SubAreaAnalysis(graph, zones, mat)\nsubarea.rc.set_choice_set_generation(\"lp\", max_routes=3, penalty=1.02, store_results=False)\nsubarea.rc.execute(perform_assignment=True)\ndemand = subarea.post_process()\ndemand"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll re-prepare our graph but with our new \"external\" ODs.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "new_centroids = np.unique(demand.reset_index()[[\"origin id\", \"destination id\"]].to_numpy().reshape(-1))\ngraph.prepare_graph(new_centroids)\ngraph.set_graph(\"utility\")\nnew_centroids"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can then perform an assignment using our new demand matrix on the limited graph\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from aequilibrae.paths import RouteChoice\n\nrc = RouteChoice(graph)\nrc.add_demand(demand)\nrc.set_choice_set_generation(\"lp\", max_routes=3, penalty=1.02, store_results=False, seed=123)\nrc.execute(perform_assignment=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's take the union of the zones GeoDataFrame as a polygon\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "poly = zones.union_all()\npoly"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And prepare the sub-area to plot.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "subarea_zone = folium.Polygon(\n    locations=[(x[1], x[0]) for x in poly.boundary.coords],\n    fill_color=\"blue\",\n    fill_opacity=0.1,\n    fill=True,\n    weight=1,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create a function to plot out link loads data more easily\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def plot_results(link_loads):\n    link_loads = link_loads[link_loads[\"demand_tot\"] > 0]\n    max_load = link_loads[[\"demand_tot\"]].max()\n    links = project.network.links.data\n    loaded_links = links.merge(link_loads, on=\"link_id\", how=\"inner\")\n    factor = 10 / max_load\n\n    return loaded_links.explore(\n        color=\"red\",\n        style_kwds={\n            \"style_function\": lambda x: {\n                \"weight\": x[\"properties\"][\"demand_tot\"] * factor,\n            }\n        },\n    )"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And plot our data!\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "map = plot_results(rc.get_load_results())\nsubarea_zone.add_to(map)\nmap"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Sub-area further preparation\n\nIt's useful later on to know which links from the network cross our polygon.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "links = project.network.links.data\ninner_links = links[links.crosses(poly.boundary)].sort_index()\ninner_links.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As well as which nodes are interior.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nodes = project.network.nodes.data.set_index(\"node_id\")\ninside_nodes = nodes.sjoin(zones, how=\"inner\").sort_index()\ninside_nodes.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's filter those network links to graph links, dropping any dead ends and creating a `link_id`,\n`dir` multi-index.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = (\n    graph.graph.set_index(\"link_id\")\n    .loc[inner_links.link_id]\n    .drop(graph.dead_end_links, errors=\"ignore\")\n    .reset_index()\n    .set_index([\"link_id\", \"direction\"])\n)\ng.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here we'll quickly visualise what our sub-area is looking like.\nWe'll plot the polygon from our zoning system and the links that it cuts.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "map = inner_links.explore(color=\"red\", style_kwds={\"weight\": 4})\nsubarea_zone.add_to(map)\nmap"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Manual sub-area analysis\n\nHere we'll construct and use the Route Choice class to generate our route sets,\n\nIn order to perform out analysis we need to know what OD pairs have flow that enters and/or exists\nour polygon. To do so we perform a select link analysis on all links and pairs of links that cross\nthe boundary. We create them as tuples of tuples to make represent the select link AND sets.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "edge_pairs = {x: (x,) for x in itertools.permutations(g.index, r=2)}\nsingle_edges = {x: ((x,),) for x in g.index}\nf\"Created: {len(edge_pairs)} edge pairs from {len(single_edges)} edges\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's prepare our graph once again\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "project.network.build_graphs()\ngraph = project.network.graphs[\"c\"]\ngraph.network = graph.network.assign(utility=graph.network.distance * theta)\ngraph.prepare_graph(graph.centroids)\ngraph.set_graph(\"utility\")\ngraph.set_blocked_centroid_flows(False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This object construction might take a minute depending on the size of the graph due to the\nconstruction of the compressed link to network link mapping that's required. This is a one\ntime operation per graph and is cached. We need to supply a Graph and an AequilibraeMatrix\nor DataFrame via the ``add_demand`` method, if demand is not provided link loading cannot\nbe preformed.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rc = RouteChoice(graph)\nrc.add_demand(mat)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here we add the union of edges as select link sets.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rc.set_select_links(single_edges | edge_pairs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the sake of demonstration we limit out demand matrix to a few OD pairs. This filter is also\npossible with the automated approach, just edit the ``subarea.rc.demand.df`` DataFrame, however\nmake sure the index remains intact.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ods_pairs_of_interest = [\n    (4, 39),\n    (92, 37),\n    (31, 58),\n    (4, 19),\n    (39, 34),\n]\nods_pairs_of_interest = ods_pairs_of_interest + [(x[1], x[0]) for x in ods_pairs_of_interest]\nrc.demand.df = rc.demand.df.loc[ods_pairs_of_interest].sort_index().astype(np.float32)\nrc.demand.df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Perform the assignment\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rc.set_choice_set_generation(\"lp\", max_routes=3, penalty=1.02, store_results=False, seed=123)\nrc.execute(perform_assignment=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can visualise the current links loads\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "map = plot_results(rc.get_load_results())\nsubarea_zone.add_to(map)\nmap"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll pull out just OD matrix results as well we need it for the post-processing, we'll also\nconvert the sparse matrices to SciPy COO matrices.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sl_od = rc.get_select_link_od_matrix_results()\nedge_totals = {k: sl_od[k][\"demand\"].to_scipy() for k in single_edges}\nedge_pair_values = {k: sl_od[k][\"demand\"].to_scipy() for k in edge_pairs}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the post processing, we are interested in the demand of OD pairs that enter or exit the\nsub-area, or do both. For the single enters and exists we can extract that information from\nthe single link select link results. We also need to map the links that cross the boundary to\nthe origin/destination node and the node that appears on the outside of the sub-area.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from collections import defaultdict\n\nentered = defaultdict(float)\nexited = defaultdict(float)\nfor (link_id, dir), v in edge_totals.items():\n    link = g.loc[link_id, dir]\n    for (o, d), load in v.todok().items():\n        o = graph.all_nodes[o]\n        d = graph.all_nodes[d]\n\n        o_inside = o in inside_nodes.index\n        d_inside = d in inside_nodes.index\n\n        if o_inside and not d_inside:\n            exited[o, graph.all_nodes[link.b_node]] += load\n        elif not o_inside and d_inside:\n            entered[graph.all_nodes[link.a_node], d] += load\n        elif not o_inside and not d_inside:\n            pass"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here he have the load that entered the sub-area\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "entered"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "and the load that exited the sub-area\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "exited"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To find the load that both entered and exited we can look at the edge pair select link results.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "through = defaultdict(float)\nfor (l1, l2), v in edge_pair_values.items():\n    link1 = g.loc[l1]\n    link2 = g.loc[l2]\n\n    for (o, d), load in v.todok().items():\n        o_inside = o in inside_nodes.index\n        d_inside = d in inside_nodes.index\n\n        if not o_inside and not d_inside:\n            through[graph.all_nodes[link1.a_node], graph.all_nodes[link2.b_node]] += load\n\nthrough"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "With these results we can construct a new demand matrix. Usually this would be now transplanted\nonto another network, however for demonstration purposes we'll reuse the same network.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "demand = pd.DataFrame(\n    list(entered.values()) + list(exited.values()) + list(through.values()),\n    index=pd.MultiIndex.from_tuples(\n        list(entered.keys()) + list(exited.keys()) + list(through.keys()), names=[\"origin id\", \"destination id\"]\n    ),\n    columns=[\"demand\"],\n).sort_index()\ndemand.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll re-prepare our graph but with our new \"external\" ODs.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "new_centroids = np.unique(demand.reset_index()[[\"origin id\", \"destination id\"]].to_numpy().reshape(-1))\ngraph.prepare_graph(new_centroids)\ngraph.set_graph(\"utility\")\nnew_centroids"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Re-perform our assignment\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rc = RouteChoice(graph)\nrc.add_demand(demand)\nrc.set_choice_set_generation(\"lp\", max_routes=3, penalty=1.02, store_results=False, seed=123)\nrc.execute(perform_assignment=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And plot the link loads for easy viewing\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "map = plot_results(rc.get_load_results())\nsubarea_zone.add_to(map)\nmap"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "project.close()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.18"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}