{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AequilibraE Routing\n",
    "\n",
    "Inputs:  demand, network\n",
    "\n",
    "Outputs:  shortest path skims, routing results\n",
    "\n",
    "## Major steps\n",
    "1. Set up Aequilibrae environment\n",
    "2. Obtain the shortest path skim from the network\n",
    "3. Run routing \n",
    "4. Generate summary statistics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aequilibrae environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#needs scipy, openmatrix (pip install)\n",
    "import sys\n",
    "from os.path import join\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import openmatrix as omx\n",
    "from math import log10, floor\n",
    "import matplotlib.pyplot as plt\n",
    "from aequilibrae.distribution import GravityCalibration, Ipf, GravityApplication, SyntheticGravityModel\n",
    "from aequilibrae import Parameters\n",
    "from aequilibrae.project import Project\n",
    "from aequilibrae.paths import PathResults\n",
    "from aequilibrae.paths import SkimResults #as skmr\n",
    "from aequilibrae.paths import Graph\n",
    "from aequilibrae.paths import NetworkSkimming\n",
    "from aequilibrae.matrix import AequilibraeData, AequilibraeMatrix\n",
    "from aequilibrae import logger\n",
    "from aequilibrae.paths import TrafficAssignment, TrafficClass\n",
    "import logging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fldr = 'C:/Users/Scott.Smith/GMNS/Lima'  #was aeqRepro\n",
    "fldr = 'C:/Users/Scott/Documents/Work/AE/Lima'  #was aeqRepro\n",
    "proj_name = 'Lima.sqlite'   #the network comes from this sqlite database\n",
    "dt_fldr = '0_tntp_data'\n",
    "prj_fldr = '1_project'\n",
    "skm_fldr = '2_skim_results'\n",
    "assg_fldr = '4_assignment_results'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = Parameters()\n",
    "p.parameters['system']['logging_directory'] = fldr\n",
    "p.write_back()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Because assignment takes a long time, we want the log to be shown here\n",
    "stdout_handler = logging.StreamHandler(sys.stdout)\n",
    "formatter = logging.Formatter(\"%(asctime)s;%(name)s;%(levelname)s ; %(message)s\")\n",
    "stdout_handler.setFormatter(formatter)\n",
    "logger.addHandler(stdout_handler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shortest path skim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2021-03-07 13:11:34,418;aequilibrae;WARNING ; AequilibraE might not work as intended without spatialite. ('The specified module could not be found.\\r\\n',)\n"
     ]
    }
   ],
   "source": [
    "project = Project()\n",
    "project.load(join(fldr, prj_fldr, proj_name))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Scott\\Anaconda3\\envs\\py37\\lib\\site-packages\\aequilibrae\\project\\network\\network.py:279: UserWarning: Fields were removed form Graph for being non-numeric: modes,link_type\n",
      "  warn(f'Fields were removed form Graph for being non-numeric: {\",\".join(removed_fields)}')\n"
     ]
    }
   ],
   "source": [
    "# we build all graphs\n",
    "project.network.build_graphs()\n",
    "# We get warnings that several fields in the project are filled with NaNs.  Which is true, but we won't\n",
    "# use those fields\n",
    "\n",
    "# we grab the graph for cars\n",
    "graph = project.network.graphs['c']\n",
    "\n",
    "# let's say we want to minimize free_flow_time #distance\n",
    "graph.set_graph('free_flow_time')\n",
    "\n",
    "# And will skim time and distance while we are at it\n",
    "graph.set_skimming(['free_flow_time', 'distance'])\n",
    "\n",
    "# And we will allow paths to be compute going through other centroids/centroid connectors\n",
    "# required for the Sioux Falls network, as all nodes are centroids\n",
    "graph.set_blocked_centroid_flows(True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "########## SKIMMING ###################\n",
    "\n",
    "# setup the object result\n",
    "res = SkimResults()\n",
    "res.prepare(graph)\n",
    "\n",
    "# And run the skimming\n",
    "res.compute_skims()\n",
    "\n",
    "# The result is an AequilibraEMatrix object\n",
    "skims = res.skims\n",
    "\n",
    "# We can export to OMX\n",
    "skims.export(join(fldr, skm_fldr, 'sp_skim.omx'))   #change for each run  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Routing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Base Skim Shape: (449, 449) Size= 449\n",
      "Number of tables 1\n",
      "Table names: ['matrix']\n",
      "attributes: []\n"
     ]
    }
   ],
   "source": [
    "#### Open the matrix to get its size ####\n",
    "f_demand = omx.open_file(join(fldr, dt_fldr, 'demand.omx'))\n",
    "matrix_shape = f_demand.shape()\n",
    "matrix_size = matrix_shape[1]\n",
    "print('Base Skim Shape:',f_demand.shape(),  \"Size=\",matrix_size)\n",
    "print('Number of tables',len(f_demand))\n",
    "print('Table names:',f_demand.list_matrices())\n",
    "print('attributes:',f_demand.list_all_attributes())\n",
    "\n",
    "f_demand.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### LOAD DEMAND MATRIX #####\n",
    "demand = AequilibraeMatrix()\n",
    "demand.load(join(fldr, dt_fldr, 'demand.omx'))\n",
    "demand.computational_view(['matrix'])  # We will only assign one user class stored as 'matrix' inside the OMX file\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2021-03-07 13:11:36,636;aequilibrae;INFO ; bfw Assignment STATS\n",
      "2021-03-07 13:11:36,636;aequilibrae;INFO ; Iteration, RelativeGap, stepsize\n",
      "2021-03-07 13:11:36,894;aequilibrae;INFO ; 1,inf,1.0\n",
      "2021-03-07 13:11:37,030;aequilibrae;WARNING ; # Alert: Adding 0.5 to stepsize to make it non-zero\n",
      "2021-03-07 13:11:37,049;aequilibrae;INFO ; 2,2.37924774254931e-06,0.5\n",
      "2021-03-07 13:11:37,175;aequilibrae;WARNING ; # Alert: Adding 0.3333333333333333 to stepsize to make it non-zero\n",
      "2021-03-07 13:11:37,180;aequilibrae;INFO ; 3,1.225687113873547e-06,0.3333333333333333\n",
      "2021-03-07 13:11:37,180;aequilibrae;INFO ; bfw Assignment finished. 3 iterations and 1.225687113873547e-06 final gap\n"
     ]
    }
   ],
   "source": [
    "######### TRAFFIC ASSIGNMENT WITH SKIMMING\n",
    "\n",
    "\n",
    "assig = TrafficAssignment()\n",
    "\n",
    "# Creates the assignment class\n",
    "assigclass = TrafficClass(graph, demand)\n",
    "\n",
    "# The first thing to do is to add at list of traffic classes to be assigned\n",
    "assig.set_classes([assigclass])\n",
    "\n",
    "assig.set_vdf(\"BPR\")  # This is not case-sensitive # Then we set the volume delay function\n",
    "\n",
    "assig.set_vdf_parameters({\"alpha\": \"b\", \"beta\": \"power\"})  # Get parameters from link file\n",
    "#assig.set_vdf_parameters({\"alpha\": 0.15, \"beta\": 4}) \n",
    "\n",
    "assig.set_capacity_field(\"capacity\")  # The capacity and free flow travel times as they exist in the graph\n",
    "assig.set_time_field(\"free_flow_time\")\n",
    "\n",
    "# And the algorithm we want to use to assign\n",
    "assig.set_algorithm('bfw')\n",
    "#assig.set_algorithm('msa')   #all-or-nothing\n",
    "\n",
    "# since I haven't checked the parameters file, let's make sure convergence criteria is good\n",
    "assig.max_iter = 100      #was 1000  or 100\n",
    "assig.rgap_target = 0.001   #was 0.00001,  or 0.01\n",
    "\n",
    "assig.execute()  # we then execute the assignment\n",
    "\n",
    "# The link flows are easy to export.\n",
    "# we do so for csv and AequilibraEData\n",
    "assigclass.results.save_to_disk(join(fldr, assg_fldr, 'linkflow.csv'), output=\"loads\") #change for each run\n",
    "#assigclass.results.save_to_disk(join(fldr, assg_fldr, 'link_flows_c1.aed'), output=\"loads\")\n",
    "\n",
    "# the skims are easy to get.\n",
    "\n",
    "# The blended one are here\n",
    "avg_skims = assigclass.results.skims\n",
    "\n",
    "# The ones for the last iteration are here\n",
    "last_skims = assigclass._aon_results.skims\n",
    "\n",
    "# Assembling a single final skim file can be done like this\n",
    "# We will want only the time for the last iteration and the distance averaged out for all iterations  \n",
    "kwargs = {'file_name': join(fldr, assg_fldr, 'rt_skim'+'.aem'),   #change\n",
    "          'zones': graph.num_zones,\n",
    "          'matrix_names': ['time_final', 'distance_blended']}\n",
    "\n",
    "# Create the matrix file\n",
    "out_skims = AequilibraeMatrix()\n",
    "out_skims.create_empty(**kwargs)\n",
    "out_skims.index[:] = avg_skims.index[:]\n",
    "\n",
    "# Transfer the data\n",
    "# The names of the skims are the name of the fields\n",
    "out_skims.matrix['time_final'][:, :] = last_skims.matrix['free_flow_time'][:, :]\n",
    "# It is CRITICAL to assign the matrix values using the [:,:]\n",
    "out_skims.matrix['distance_blended'][:, :] = avg_skims.matrix['distance'][:, :]\n",
    "\n",
    "out_skims.matrices.flush()  # Make sure that all data went to the disk\n",
    "\n",
    "# Export to OMX as well\n",
    "out_skims.export(join(fldr, assg_fldr, 'rt_skim'+'.omx'))\n",
    "demand.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate summary statistics\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DEMAND FILE Shape: (449, 449)    Tables: ['matrix']    Mappings: ['taz']\n",
      "SP BASE SKIM FILE Shape: (449, 449)    Tables: ['distance', 'free_flow_time']    Mappings: ['main_index']\n",
      "RT BASE SKIM FILE Shape: (449, 449)    Tables: ['distance_blended', 'time_final']    Mappings: ['main_index']\n"
     ]
    }
   ],
   "source": [
    "f = omx.open_file(join(fldr, dt_fldr, 'demand.omx'),'r')  #change\n",
    "print('DEMAND FILE Shape:',f.shape(),'   Tables:',f.list_matrices(),'   Mappings:',f.list_mappings())\n",
    "dem = f['matrix']\n",
    "\n",
    "spbf = omx.open_file(join(fldr, skm_fldr,'sp_skim.omx'),'r')    #change\n",
    "print('SP BASE SKIM FILE Shape:',spbf.shape(),'   Tables:',spbf.list_matrices(),'   Mappings:',spbf.list_mappings())\n",
    "spbt = spbf['free_flow_time']\n",
    "spbd = spbf['distance']\n",
    "\n",
    "rtbf = omx.open_file(join(fldr, assg_fldr, 'rt_skim.omx'),'r')\n",
    "print('RT BASE SKIM FILE Shape:',rtbf.shape(),'   Tables:',rtbf.list_matrices(),'   Mappings:',rtbf.list_mappings())\n",
    "\n",
    "rtbt = rtbf['time_final']\n",
    "rtbd = rtbf['distance_blended']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sum of demand trips 32041.0\n"
     ]
    }
   ],
   "source": [
    "#Summary information on the input trip tables\n",
    "print('sum of demand trips','{:.9}'.format(np.sum(dem)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Skims as .csv files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "32041.0 735395949.0 735382285.6666669 18659217.216666665 18677288.385667212\n"
     ]
    }
   ],
   "source": [
    "outfile = open(\"combined_skim.txt\",\"w\")  #change\n",
    "spb_cumtripcount = 0.0;\n",
    "spb_cumtime = 0.0;\n",
    "spb_cumdist = 0.0;\n",
    "rtb_cumtime = 0.0;\n",
    "rtb_cumdist = 0.0;\n",
    "largeval = 999999;\n",
    "\n",
    "#Shortest path base times and distances\n",
    "print(\"i j demand sp_dist rt_dist sp_time rt_time\",file=outfile)\n",
    "for i in range(matrix_size):\n",
    "    tripcount = 0.0;\n",
    "    sp_timecount = 0.0;\n",
    "    sp_distcount = 0.0;\n",
    "    rt_timecount = 0.0;\n",
    "    rt_distcount = 0.0;\n",
    "    for j in range(matrix_size):\n",
    "        if(dem[i][j]>0):\n",
    "            tripcount = tripcount + dem[i][j]\n",
    "            sp_timecount = sp_timecount + dem[i][j]*spbt[i][j]\n",
    "            sp_distcount = sp_distcount + dem[i][j]*spbd[i][j]\n",
    "            rt_timecount = rt_timecount + dem[i][j]*rtbt[i][j]\n",
    "            rt_distcount = rt_distcount + dem[i][j]*rtbd[i][j]\n",
    "            print(i,j,dem[i][j],spbd[i][j],rtbd[i][j],spbt[i][j],rtbt[i][j],file=outfile)\n",
    "    #print(\"SP Base Row\",i,'{:.6} {:.6} {:.6}'.format(tripcount,distcount,timecount),file=outfile)\n",
    "    spb_cumtripcount = spb_cumtripcount + tripcount;\n",
    "    spb_cumtime = spb_cumtime + sp_timecount;\n",
    "    spb_cumdist = spb_cumdist + sp_distcount;\n",
    "    rtb_cumtime = rtb_cumtime + rt_timecount;\n",
    "    rtb_cumdist = rtb_cumdist + rt_distcount;\n",
    "    #print(\"Row\",i,tripcount,timecount,distcount)\n",
    "#print(\"Shortest path base totals\",'{:.8} {:.8} {:.8}'.format(cumtripcount,cumdist,cumtime),file=outfile)\n",
    "#print(\"Shortest path base totals\",'{:.8} {:.8} {:.8}'.format(spb_cumtripcount,spb_cumdist,spb_cumtime))\n",
    "print(spb_cumtripcount,spb_cumdist,rtb_cumdist,spb_cumtime/60,rtb_cumtime/60)\n",
    "outfile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative calculations using numpy array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total pht 18659217.216666665   average per trip 582.3543964503813\n",
      "total pmt 735395949.0   average per trip 22951.71651945944\n"
     ]
    }
   ],
   "source": [
    "sp_pht = np.array(dem)*np.array(spbt)/60\n",
    "sp_pmt = np.array(dem)*np.array(spbd)\n",
    "print('total pht',np.sum(sp_pht),'  average per trip',np.sum(sp_pht)/np.sum(dem))\n",
    "print('total pmt',np.sum(sp_pmt),'  average per trip',np.sum(sp_pmt)/np.sum(dem))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total pht 18677288.38566721   average per trip 582.9183978548488\n",
      "total pmt 735382285.6666669   average per trip 22951.290086659807\n"
     ]
    }
   ],
   "source": [
    "rt_pht = np.array(dem)*np.array(rtbt)/60\n",
    "rt_pmt = np.array(dem)*np.array(rtbd)\n",
    "print('total pht',np.sum(rt_pht),'  average per trip',np.sum(rt_pht)/np.sum(dem))\n",
    "print('total pmt',np.sum(rt_pmt),'  average per trip',np.sum(rt_pmt)/np.sum(dem))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Close the files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.close()\n",
    "spbf.close()\n",
    "rtbf.close()\n",
    "outfile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}