{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Trip Distribution\n\nIn this example, we calibrate a Synthetic Gravity Model that same model plus IPF (Fratar/Furness).\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\nfrom aequilibrae.utils.create_example import create_example\nimport pandas as pd\nimport numpy as np"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We create the example project inside our temp folder\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fldr = join(gettempdir(), uuid4().hex)\n\nproject = create_example(fldr)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# We get the demand matrix directly from the project record\n# so let's inspect what we have in the project\nproj_matrices = project.matrices\nprint(proj_matrices.list())"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# We get the demand matrix\ndemand = proj_matrices.get_matrix(\"demand_omx\")\ndemand.computational_view([\"matrix\"])\n\n# And the impedance\nimpedance = proj_matrices.get_matrix(\"skims\")\nimpedance.computational_view([\"time_final\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's have a function to plot the Trip Length Frequency Distribution\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from math import log10, floor\nimport matplotlib.pyplot as plt\n\n\ndef plot_tlfd(demand, skim, name):\n    plt.clf()\n    b = floor(log10(skim.shape[0]) * 10)\n    n, bins, patches = plt.hist(\n        np.nan_to_num(skim.flatten(), 0),\n        bins=b,\n        weights=np.nan_to_num(demand.flatten()),\n        density=False,\n        facecolor=\"g\",\n        alpha=0.75,\n    )\n\n    plt.xlabel(\"Trip length\")\n    plt.ylabel(\"Probability\")\n    plt.title(\"Trip-length frequency distribution\")\n    plt.savefig(name, format=\"png\")\n    return plt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from aequilibrae.distribution import GravityCalibration"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for function in [\"power\", \"expo\"]:\n    gc = GravityCalibration(matrix=demand, impedance=impedance, function=function, nan_as_zero=True)\n    gc.calibrate()\n    model = gc.model\n    # We save the model\n    model.save(join(fldr, f\"{function}_model.mod\"))\n\n    # We can save an image for the resulting model\n    _ = plot_tlfd(gc.result_matrix.matrix_view, impedance.matrix_view, join(fldr, f\"{function}_tfld.png\"))\n\n    # We can save the result of applying the model as well\n    # We can also save the calibration report\n    with open(join(fldr, f\"{function}_convergence.log\"), \"w\") as otp:\n        for r in gc.report:\n            otp.write(r + \"\\n\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We save a trip length frequency distribution for the demand itself\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt = plot_tlfd(demand.matrix_view, impedance.matrix_view, join(fldr, \"demand_tfld.png\"))\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Forecast\nWe create a set of *'future'* vectors by applying some models\nand apply the model for both deterrence functions\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from aequilibrae.distribution import Ipf, GravityApplication, SyntheticGravityModel\nfrom aequilibrae.matrix import AequilibraeData\nimport numpy as np"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "zonal_data = pd.read_sql(\"Select zone_id, population, employment from zones order by zone_id\", project.conn)\n# We compute the vectors from our matrix\n\nargs = {\n    \"file_path\": join(fldr, \"synthetic_future_vector.aed\"),\n    \"entries\": demand.zones,\n    \"field_names\": [\"origins\", \"destinations\"],\n    \"data_types\": [np.float64, np.float64],\n    \"memory_mode\": True,\n}\n\nvectors = AequilibraeData()\nvectors.create_empty(**args)\n\nvectors.index[:] = zonal_data.zone_id[:]\n\n# We apply a trivial regression-based model and balance the vectors\nvectors.origins[:] = zonal_data.population[:] * 2.32\nvectors.destinations[:] = zonal_data.employment[:] * 1.87\nvectors.destinations *= vectors.origins.sum() / vectors.destinations.sum()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# We simply apply the models to the same impedance matrix now\nfor function in [\"power\", \"expo\"]:\n    model = SyntheticGravityModel()\n    model.load(join(fldr, f\"{function}_model.mod\"))\n\n    outmatrix = join(proj_matrices.fldr, f\"demand_{function}_model.aem\")\n    args = {\n        \"impedance\": impedance,\n        \"rows\": vectors,\n        \"row_field\": \"origins\",\n        \"model\": model,\n        \"columns\": vectors,\n        \"column_field\": \"destinations\",\n        \"nan_as_zero\": True,\n    }\n\n    gravity = GravityApplication(**args)\n    gravity.apply()\n\n    # We get the output matrix and save it to OMX too,\n    gravity.save_to_project(name=f\"demand_{function}_model_omx\", file_name=f\"demand_{function}_model.omx\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# We update the matrices table/records and verify that the new matrices are indeed there\nproj_matrices.update_database()\nprint(proj_matrices.list())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now run IPF for the future vectors\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "args = {\n    \"matrix\": demand,\n    \"rows\": vectors,\n    \"columns\": vectors,\n    \"column_field\": \"destinations\",\n    \"row_field\": \"origins\",\n    \"nan_as_zero\": True,\n}\n\nipf = Ipf(**args)\nipf.fit()\n\nipf.save_to_project(name=\"demand_ipf\", file_name=\"demand_ipf.aem\")\nipf.save_to_project(name=\"demand_ipf_omx\", file_name=\"demand_ipf.omx\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(proj_matrices.list())"
      ]
    },
    {
      "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.9.16"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}