docs/walkthrough/sparse.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c871b4c1",
   "metadata": {},
   "source": [
    "# Using Sparse MAZ Skims"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8401e0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import sharrow as sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20bc86e7",
   "metadata": {},
   "source": [
    "This notebook walks through using sparse MAZ to MAZ skims with sharrow.\n",
    "The example data we'll use to demonstrate this feature starts with regular\n",
    "TAZ-based skims."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8982a01",
   "metadata": {},
   "outputs": [],
   "source": [
    "skims = sh.example_data.get_skims()\n",
    "skims"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35396a40",
   "metadata": {},
   "source": [
    "We'll also load a MAZ-to-TAZ mapping file, which defines the MAZ's and \n",
    "which TAZ is used for each MAZ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dae3b485",
   "metadata": {},
   "outputs": [],
   "source": [
    "maz_taz = sh.example_data.get_maz_to_taz()\n",
    "maz_taz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4639626f",
   "metadata": {},
   "source": [
    "Lastly, we'll load a sparse MAZ-to-MAZ skim table.  This table\n",
    "defines origin and destination MAZ's, and the walk distance\n",
    "between them.  The data is \"sparse\" in that only a limited number\n",
    "of OMAZ-DMAZ pairs are included.  Unlike traditional sparse arrays,\n",
    "the missing elements are not assumed to be zero, but instead we \n",
    "implicitly use the walk distance from the matching TAZ's in the\n",
    "TAZ-based skims for those zone pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e10b2e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "maz_to_maz_walk = sh.example_data.get_maz_to_maz_walk()\n",
    "maz_to_maz_walk"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8eb5269",
   "metadata": {},
   "source": [
    "To integrate these data sources, we will set a redirection on the skims.\n",
    "This will add the MAZ dimensions to the skims, MAZ id's as additional \n",
    "coordinates, and will set attribute flags to tell sharrow which dimensions\n",
    "have been redirected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8071c236",
   "metadata": {},
   "outputs": [],
   "source": [
    "skims.redirection.set(\n",
    "    maz_taz, \n",
    "    map_to='otaz', \n",
    "    name=\"omaz\",\n",
    "    map_also={'dtaz': \"dmaz\"}, \n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4323483",
   "metadata": {},
   "outputs": [],
   "source": [
    "skims"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1105ef9",
   "metadata": {},
   "source": [
    "Next, we can attach the sparse skims using `redirection.sparse_blender`.\n",
    "This formats the sparse skim table into compressed sparse row format,\n",
    "and attaches the resulting arrays to the Dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2269a8c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "skims.redirection.sparse_blender(\n",
    "    'DISTWALK', \n",
    "    maz_to_maz_walk.OMAZ, \n",
    "    maz_to_maz_walk.DMAZ, \n",
    "    maz_to_maz_walk.DISTWALK,\n",
    "    max_blend_distance=1.0,\n",
    "    index=maz_taz.index,\n",
    ")\n",
    "skims"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff55a65d",
   "metadata": {},
   "source": [
    "Now the skims are ready to use!\n",
    "\n",
    "For demonstration purposes, let's construct a trips dataframe with just a few\n",
    "origin-destination pairs. Note that we're using the zone id's from the more \n",
    "detailed MAZ system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4eaf4dbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "trips = pd.DataFrame({\n",
    "    'orig_maz': [100, 100, 100, 200, 200],\n",
    "    'dest_maz': [100, 101, 103, 201, 202],\n",
    "})\n",
    "trips"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "998d16be",
   "metadata": {},
   "source": [
    "We'll then put the trips together with the skims into a DataTree, as\n",
    "usual for sharrow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60a6848c",
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = sh.DataTree(\n",
    "    base=trips,\n",
    "    skims=skims,\n",
    "    relationships=(\n",
    "        \"base.orig_maz @ skims.omaz\",\n",
    "        \"base.dest_maz @ skims.dmaz\",\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3523e22",
   "metadata": {},
   "source": [
    "Now we can setup flows on this tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be2e86f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "flow = tree.setup_flow({\n",
    "    'plain_distance': 'DISTWALK',\n",
    "}, boundscheck=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a43c2b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "flow.load()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d7b3902",
   "metadata": {},
   "source": [
    "Where the sparse (maz) data is missing or exceeds the max blending distance,\n",
    "the dense (taz) data is returned.  Otherwise, the output is not strictly taken \n",
    "from the sparse or dense skims, but it is a blended mixture of the two. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa6dc462",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TEST\n",
    "from pytest import approx\n",
    "sparse_dat = np.array([0.01, 0.2, np.nan, 3.2, np.nan])\n",
    "dense_dat = np.array([0.12,0.12,0.12,0.17,0.17])\n",
    "def blend(s,d, max_s):\n",
    "    out = np.zeros_like(d)\n",
    "    ratio = s/max_s\n",
    "    out = d*ratio + s*(1-ratio)\n",
    "    out = np.where(s>max_s, d, out)\n",
    "    out = np.where(np.isnan(s), d, out)\n",
    "    return out\n",
    "assert blend(sparse_dat, dense_dat, 1.0) == approx(flow.load().ravel())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b4604e3",
   "metadata": {},
   "source": [
    "We can apply all the transformation we like, as usual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b9aa948",
   "metadata": {},
   "outputs": [],
   "source": [
    "flow2 = tree.setup_flow({\n",
    "    'plain_distance': 'DISTWALK',\n",
    "    'clip_distance': 'DISTWALK.clip(upper=0.15)',\n",
    "    'square_distance': 'DISTWALK**2',\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1526e25b",
   "metadata": {},
   "outputs": [],
   "source": [
    "flow2.load_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2fd8300",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TEST\n",
    "assert flow2.load_dataframe().values == approx(np.array([\n",
    "    [ 1.1100e-02,  1.1100e-02,  1.2321e-04],\n",
    "    [ 1.8400e-01,  1.5000e-01,  3.3856e-02],\n",
    "    [ 1.2000e-01,  1.2000e-01,  1.4400e-02],\n",
    "    [ 1.7000e-01,  1.5000e-01,  2.8900e-02],\n",
    "    [ 1.7000e-01,  1.5000e-01,  2.8900e-02]], dtype=np.float32)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0d40e91",
   "metadata": {},
   "source": [
    "## Using at and iat\n",
    "\n",
    "The `at` and `iat` accessors work even when sparse matrix tables are\n",
    "attached to a Dataset, with a few caveats.  First, only 2-dimension\n",
    "sparse tables are supported at this time.  Second, these accessors \n",
    "rely on the ability to reference the sparse data, which is lost if \n",
    "the dataset is naively filtered for variable names; filtering should\n",
    "instead be done in the `_names` argument, which filters the \n",
    "output of the accessor instead of the input, without needing to build\n",
    "the entire filtered dataset first.  For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a3971fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "skims.at(\n",
    "    omaz=trips.orig_maz,\n",
    "    dmaz=trips.dest_maz,\n",
    "    _names=['DIST', 'DISTWALK'],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48c42f7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TEST\n",
    "out = skims.at(\n",
    "    omaz=trips.orig_maz,\n",
    "    dmaz=trips.dest_maz,\n",
    "    _names=['DIST', 'DISTWALK'], _load=True,\n",
    ")\n",
    "np.testing.assert_array_almost_equal(\n",
    "    out['DIST'].to_numpy(), \n",
    "    np.array([0.12, 0.12, 0.12, 0.17, 0.17], dtype=np.float32)\n",
    ")\n",
    "np.testing.assert_array_almost_equal(\n",
    "    out['DISTWALK'].to_numpy(), \n",
    "    np.array([0.0111, 0.184, 0.12,  0.17, 0.17], dtype=np.float32)\n",
    ")\n",
    "\n",
    "from pytest import raises\n",
    "with raises(NotImplementedError):\n",
    "    skims.at(\n",
    "        omaz=trips.orig_maz,\n",
    "        dmaz=trips.dest_maz,\n",
    "        time_period=['AM', 'AM', 'AM', 'AM', 'AM'],\n",
    "        _names=['DIST', 'DISTWALK', 'SOV_TIME'], _load=True,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20f06ab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "skims.iat(\n",
    "    omaz=[  0,   0,   0, 100, 100],\n",
    "    dmaz=[  0,   1,   3, 101, 102],\n",
    "    _names=['DIST', 'DISTWALK'],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eaf576ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TEST\n",
    "out = skims.iat(\n",
    "    omaz=[  0,   0,   0, 100, 100],\n",
    "    dmaz=[  0,   1,   3, 101, 102],\n",
    "    _names=['DIST', 'DISTWALK'], _load=True,\n",
    ")\n",
    "np.testing.assert_array_almost_equal(\n",
    "    out['DIST'].to_numpy(), \n",
    "    np.array([0.12, 0.12, 0.12, 0.17, 0.17], dtype=np.float32)\n",
    ")\n",
    "np.testing.assert_array_almost_equal(\n",
    "    out['DISTWALK'].to_numpy(), \n",
    "    np.array([0.0111, 0.184, 0.12,  0.17, 0.17], dtype=np.float32)\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68dcc914",
   "metadata": {},
   "source": [
    "To circumvent the redirection, and sparse lookup and blending,\n",
    "simply point the accessor lookups to the dense dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9117a3d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "skims.at(\n",
    "    otaz=[1,1,1,16,16],\n",
    "    dtaz=[1,1,1,16,16],\n",
    "    _names=['DIST', 'DISTWALK'], _load=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34b334d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "skims.at(\n",
    "    otaz=[1,1,1,16,16],\n",
    "    dtaz=[1,1,1,16,16],\n",
    "    _name='DISTWALK',\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1078ea4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TEST\n",
    "import sys\n",
    "if sys.version_info > (3,8):\n",
    "    import secrets\n",
    "    token = \"skims-with-sparse\" + secrets.token_hex(5)\n",
    "    readback0 = skims.shm.to_shared_memory(token)\n",
    "    assert readback0.attrs == skims.attrs\n",
    "    readback = sh.Dataset.shm.from_shared_memory(token)\n",
    "    assert readback.attrs == skims.attrs\n",
    "    \n",
    "    out = readback.iat(\n",
    "        omaz=[  0,   0,   0, 100, 100],\n",
    "        dmaz=[  0,   1,   3, 101, 102],\n",
    "        _names=['DIST', 'DISTWALK'], _load=True,\n",
    "    )\n",
    "    np.testing.assert_array_almost_equal(\n",
    "        out['DIST'].to_numpy(), \n",
    "        np.array([0.12, 0.12, 0.12, 0.17, 0.17], dtype=np.float32)\n",
    "    )\n",
    "    np.testing.assert_array_almost_equal(\n",
    "        out['DISTWALK'].to_numpy(), \n",
    "        np.array([0.0111, 0.184, 0.12,  0.17, 0.17], dtype=np.float32)\n",
    "    )\n",
    "\n",
    "    out = readback.at(\n",
    "        omaz=trips.orig_maz,\n",
    "        dmaz=trips.dest_maz,\n",
    "        _names=['DIST', 'DISTWALK'], _load=True,\n",
    "    )\n",
    "    np.testing.assert_array_almost_equal(\n",
    "        out['DIST'].to_numpy(), \n",
    "        np.array([0.12, 0.12, 0.12, 0.17, 0.17], dtype=np.float32)\n",
    "    )\n",
    "    np.testing.assert_array_almost_equal(\n",
    "        out['DISTWALK'].to_numpy(), \n",
    "        np.array([0.0111, 0.184, 0.12,  0.17, 0.17], dtype=np.float32)\n",
    "    )\n",
    "    \n",
    "    assert readback.redirection.blenders == {'DISTWALK': {'max_blend_distance': 1.0, 'blend_distance_name': None}}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a25cafa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TEST\n",
    "assert skims.redirection.blenders == {'DISTWALK': {'max_blend_distance': 1.0, 'blend_distance_name': None}}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96a7c554",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TEST\n",
    "# reverse skims in sparse\n",
    "flow3 = tree.setup_flow({\n",
    "    'plain_distance': 'DISTWALK',\n",
    "    'reverse_distance': 'skims.reverse(\"DISTWALK\")',\n",
    "})\n",
    "\n",
    "assert flow3.load() == approx(np.array([[ 0.0111,  0.0111],\n",
    "       [ 0.184 ,  0.12  ],\n",
    "       [ 0.12  ,  0.12  ],\n",
    "       [ 0.17  ,  0.17  ],\n",
    "       [ 0.17  ,  0.17  ]], dtype=np.float32))\n",
    "\n",
    "z = skims.iat(\n",
    "    omaz=[  0,   1,   3, 101, 102],\n",
    "    dmaz=[  0,   0,   0, 100, 100],\n",
    "    _names=['DIST', 'DISTWALK'], _load=True,\n",
    ")\n",
    "assert z['DISTWALK'].data == approx(np.array([ 0.0111,  0.12  ,  0.12  ,  0.17  ,  0.17  ]))\n",
    "assert z['DIST'].data == approx(np.array([ 0.12,  0.12  ,  0.12  ,  0.17  ,  0.17  ]))"
   ]
  }
 ],
 "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.9.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "264.105px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}