code/X_paper_figures_multicity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extra notebook - Paper figures, multi-city analysis\n",
    "## Project: Growing Urban Bicycle Networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook creates the figures for the research paper and runs further multi-city analysis. Before running this notebook, the notebooks 03, 04, and 05 need to have been run for all 6 parameter sets, for all the cities in cities_all.csv. Results on car-constricted networks, if available from X_carconstricted_metrics, can be plotted by \n",
    "setting `plotconstricted = True`.\n",
    "\n",
    "Contact: Michael Szell (michael.szell@gmail.com)  \n",
    "Created: 2021-02-08  \n",
    "Last modified: 2024-09-24"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preliminaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "debug = False # If True, will produce plots and/or verbose output to double-check\n",
    "%run -i \"../parameters/parameters.py\"\n",
    "plotnotitle = True # If True, will not plot a text title for the plots on top left\n",
    "plotconstricted = False # If True, will add plots about constricted street network metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run -i path.py\n",
    "%run -i setup.py\n",
    "plt.style.use(PATH[\"parameters\"] + 'plotstyle.mplstyle')\n",
    "\n",
    "%load_ext watermark\n",
    "%watermark -n -v -m -g -iv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run -i functions.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run all parameter sets\n",
    "poi_source_list = [\"grid\", \"railwaystation\"]\n",
    "prune_measure_list = [\"betweenness\", \"closeness\", \"random\"]\n",
    "combs = list(itertools.product(poi_source_list, prune_measure_list))\n",
    "print(combs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "numcitiestotal = len(cities.keys())\n",
    "numcitiestotal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load all results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_result = {}\n",
    "for p in poi_source_list:\n",
    "    analysis_result[p] = {}\n",
    "    for m in prune_measure_list:\n",
    "        analysis_result[p][m] = {}\n",
    "        \n",
    "        for placeid, placeinfo in tqdm(cities.items(), desc=\"Cities\"):\n",
    "            filename = placeid + '_poi_' + p + \"_\" + m + \".csv\"\n",
    "            analysis_result[p][m][placeid] = np.genfromtxt(PATH[\"results\"] + placeid + \"/\" + filename, delimiter=',', names=True)\n",
    "            if len(analysis_result[p][m][placeid]) == 0:\n",
    "                analysis_result[p][m][placeid] = analysis_result[p][m][list(cities.keys())[0]]\n",
    "                for n in analysis_result[p][m][placeid].dtype.names:\n",
    "                    analysis_result[p][m][placeid][n] = [-1]*len(analysis_result[p][m][placeid][n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_existing = {}\n",
    "for placeid, placeinfo in tqdm(cities.items(), desc=\"Cities\"):\n",
    "    filename = placeid + '_existing.csv'\n",
    "    analysis_existing[placeid] = np.genfromtxt(PATH[\"results\"] + placeid + \"/\" + filename, delimiter=',', names=True, usecols = (1,2,3,4,5,6,7,8,9,10,11))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_constricted = {}\n",
    "if plotconstricted:\n",
    "    for p in poi_source_list:\n",
    "        analysis_constricted[p] = {}\n",
    "        for m in prune_measure_list:\n",
    "            analysis_constricted[p][m] = {}\n",
    "            \n",
    "            for placeid, placeinfo in tqdm(cities.items(), desc=\"Cities\"):\n",
    "                f = PATH[\"results_constricted\"] + \"results_\" + p + \"_\" + m + \"/metrics_\" + p + \"_\" + m + \"/\" + placeid + \"_carconstrictedbike_poi_\" + p + \"_\" + m + \".csv\"\n",
    "                if os.path.isfile(f):\n",
    "                    temp = np.loadtxt(f, delimiter=',', usecols = (2,3,4,5,6,7,8,9,10), skiprows=2)\n",
    "                    if np.shape(temp)[0] != 2: # we dont consider large cities for the avg where we only calculated 2 values\n",
    "                        analysis_constricted[p][m][placeid] = temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Directness and Efficiency for all cities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we want to check the hypothesis that there is a \"dip\" in the metrics of directness and efficiency, i.e. that there is a U-shaped form where the metric starts at a high value, then falls  (due to percolation / emergence of the giant component) and then grows back to an intermediate value. We do that for each combination of [('grid', 'betweenness'), ('grid', 'closeness'), ('railwaystation', 'betweenness'), ('railwaystation', 'closeness')]\n",
    "\n",
    "To do that, first we select all cities where the minimum values of the metric are to the right of the maximum value. For those cities we then only plot the (x,y) pairs of the min, max, and end values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Directness min/max/end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "directness_lcc = {}\n",
    "for p in poi_source_list:\n",
    "    directness_lcc[p] = {}\n",
    "    for m in prune_measure_list:\n",
    "        directness_lcc[p][m] = {}\n",
    "        directness_lcc[p][m][\"x\"] = {}\n",
    "        directness_lcc[p][m][\"y\"] = {}\n",
    "\n",
    "        directness_lcc[p][m][\"y\"][\"min\"] = [min(analysis_result[p][m][placeid][\"directness_lcc\"]) for placeid in cities.keys()]\n",
    "        directness_lcc[p][m][\"y\"][\"max\"] = [max(analysis_result[p][m][placeid][\"directness_lcc\"]) for placeid in cities.keys()]\n",
    "        directness_lcc[p][m][\"y\"][\"end\"] = [analysis_result[p][m][placeid][\"directness_lcc\"][-1] for placeid in cities.keys()]\n",
    "\n",
    "        directness_lcc[p][m][\"x\"][\"min\"] = [np.where(analysis_result[p][m][placeid][\"directness_lcc\"] == min(analysis_result[p][m][placeid][\"directness_lcc\"]))[0][-1] for placeid in cities.keys()]\n",
    "        directness_lcc[p][m][\"x\"][\"max\"] = [np.where(analysis_result[p][m][placeid][\"directness_lcc\"] == max(analysis_result[p][m][placeid][\"directness_lcc\"]))[0][-1] for placeid in cities.keys()]\n",
    "        directness_lcc[p][m][\"x\"][\"end\"] = [np.where(analysis_result[p][m][placeid][\"directness_lcc\"] == analysis_result[p][m][placeid][\"directness_lcc\"][-1])[0][-1] for placeid in cities.keys()]       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows = 2, ncols = 2, figsize = (10, 8), squeeze = True)\n",
    "axes = axes.flatten()\n",
    "\n",
    "# Plot data\n",
    "for i,ax in enumerate(axes):\n",
    "    ind = np.where(np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"x\"][\"min\"]) > np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"x\"][\"max\"]))[0]\n",
    "    print(str(len(ind)) + \"/\" + str(numcitiestotal) + \" cities found with x_min>x_max for \"+ str(combs[i]))\n",
    "    ax.plot([np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"x\"][\"max\"])[ind], np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"x\"][\"min\"])[ind]], [np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"y\"][\"max\"])[ind], np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"y\"][\"min\"])[ind]], ':', color=\"red\", alpha=0.3);\n",
    "    \n",
    "#     ind = np.where(np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"y\"][\"min\"]) >= 0)[0]\n",
    "    ax.plot([np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"x\"][\"min\"])[ind], np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"x\"][\"end\"])[ind]], [np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"y\"][\"min\"])[ind], np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"y\"][\"end\"])[ind]], ':', color=\"green\", alpha=0.3);\n",
    "    \n",
    "    ax.plot(np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"x\"][\"max\"])[ind], np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"y\"][\"max\"])[ind], '^r', label='max');\n",
    "    ax.plot(np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"x\"][\"min\"])[ind], np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"y\"][\"min\"])[ind], 'vg', label='min');\n",
    "    ax.plot(np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"x\"][\"end\"])[ind], np.asarray(directness_lcc[combs[i][0]][combs[i][1]][\"y\"][\"end\"])[ind], 'ok', label='end', markerfacecolor='none');\n",
    "    \n",
    "    ax.set_xlabel(combs[i][1] + ' quantile')\n",
    "    ax.set_ylabel('Directness of LCC')\n",
    "    ax.set_xlim([-1,40])\n",
    "    ax.set_ylim([0.4,1]) #0.35\n",
    "    if i == 0: ax.legend(loc='lower right');\n",
    "    if i == 1 or i == 3:\n",
    "        ax.set_ylabel('')\n",
    "        ax.set_yticklabels([])\n",
    "        ax.text(42, 0.8, combs[i][0], rotation = 90, horizontalalignment = \"center\", verticalalignment='center')\n",
    "    if i == 0 or i == 1:\n",
    "        ax.set_xlabel('')\n",
    "        ax.set_xticklabels([])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Efficiency min/max/end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "efficiency_global = {}\n",
    "for p in poi_source_list:\n",
    "    efficiency_global[p] = {}\n",
    "    for m in prune_measure_list:\n",
    "        efficiency_global[p][m] = {}\n",
    "        efficiency_global[p][m][\"x\"] = {}\n",
    "        efficiency_global[p][m][\"y\"] = {}\n",
    "\n",
    "        efficiency_global[p][m][\"y\"][\"min\"] = [min(analysis_result[p][m][placeid][\"efficiency_global\"]) for placeid in cities.keys()]\n",
    "        efficiency_global[p][m][\"y\"][\"max\"] = [max(analysis_result[p][m][placeid][\"efficiency_global\"]) for placeid in cities.keys()]\n",
    "        efficiency_global[p][m][\"y\"][\"end\"] = [analysis_result[p][m][placeid][\"efficiency_global\"][-1] for placeid in cities.keys()]\n",
    "\n",
    "        efficiency_global[p][m][\"x\"][\"min\"] = [np.where(analysis_result[p][m][placeid][\"efficiency_global\"] == min(analysis_result[p][m][placeid][\"efficiency_global\"]))[0][-1] for placeid in cities.keys()]\n",
    "        efficiency_global[p][m][\"x\"][\"max\"] = [np.where(analysis_result[p][m][placeid][\"efficiency_global\"] == max(analysis_result[p][m][placeid][\"efficiency_global\"]))[0][-1] for placeid in cities.keys()]\n",
    "        efficiency_global[p][m][\"x\"][\"end\"] = [np.where(analysis_result[p][m][placeid][\"efficiency_global\"] == analysis_result[p][m][placeid][\"efficiency_global\"][-1])[0][-1] for placeid in cities.keys()]       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows = 2, ncols = 2, figsize = (10, 8), squeeze = True)\n",
    "axes = axes.flatten()\n",
    "\n",
    "# Plot data\n",
    "for i,ax in enumerate(axes):\n",
    "   \n",
    "    ind = np.where(np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"min\"]) > np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"max\"]))[0]\n",
    "    print(str(len(ind)) + \"/\" + str(numcitiestotal) + \" cities found with x_min>x_max for \"+ str(combs[i]))\n",
    "    ax.plot([np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"max\"])[ind], np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"min\"])[ind]], [np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"y\"][\"max\"])[ind], np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"y\"][\"min\"])[ind]], ':', color=\"red\", alpha=0.3);\n",
    "    \n",
    "#     ind = np.where(np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"y\"][\"min\"]) >= 0)[0]\n",
    "    ax.plot([np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"min\"])[ind], np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"end\"])[ind]], [np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"y\"][\"min\"])[ind], np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"y\"][\"end\"])[ind]], ':', color=\"green\", alpha=0.3);\n",
    "    \n",
    "    ax.plot(np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"max\"])[ind], np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"y\"][\"max\"])[ind], '^r', label='max');\n",
    "    ax.plot(np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"min\"])[ind], np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"y\"][\"min\"])[ind], 'vg', label='min');\n",
    "    ax.plot(np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"end\"])[ind], np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"y\"][\"end\"])[ind], 'ok', label='end', markerfacecolor='none');\n",
    "    \n",
    "    ax.set_xlabel(combs[i][1] + ' quantile')\n",
    "    ax.set_ylabel('Global Efficiency')\n",
    "    ax.set_xlim([-1,40])\n",
    "    ax.set_ylim([0,1])\n",
    "    if i == 0: ax.legend(loc='lower center');\n",
    "    if i == 1 or i == 3:\n",
    "        ax.set_ylabel('')\n",
    "        ax.set_yticklabels([])\n",
    "        ax.text(42, 0.55, combs[i][0], rotation = 90, horizontalalignment = \"center\", verticalalignment='center')\n",
    "    if i == 0 or i == 1:\n",
    "        ax.set_xlabel('')\n",
    "        ax.set_xticklabels([])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting full curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows = 1, ncols = 4, figsize=(1200/plotparam[\"dpi\"], 300/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"], squeeze = True)\n",
    "axes = axes.flatten()\n",
    "\n",
    "for i, comb in enumerate(combs):\n",
    "    p = comb[0]\n",
    "    m = comb[1]\n",
    "    ind = np.where(np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"min\"]) > np.asarray(efficiency_global[combs[i][0]][combs[i][1]][\"x\"][\"max\"]))[0]\n",
    "    print(str(len(ind)) + \"/\" + str(numcitiestotal) + \" cities found with x_min>x_max for \"+ str(combs[i]))\n",
    "    print([list(cities.keys())[i] for i in ind])\n",
    "    \n",
    "    ax = axes[i]\n",
    "\n",
    "    for i in ind:\n",
    "        placeid = list(cities.keys())[i]\n",
    "        tmp, = ax.plot(prune_quantiles, analysis_result[p][m][placeid][\"efficiency_global_routed\"])\n",
    "        tmp.set_label('_hidden')\n",
    "    ax.set_title(comb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next: Plotting single cities with the biggest dips for only grid,betweenness and railwaystation, betweenness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[combs[0],combs[2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows = 2, ncols = 15, figsize=(1800/plotparam[\"dpi\"], 300/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"], squeeze = True)\n",
    "axes = axes.flatten()\n",
    "\n",
    "for i, comb in enumerate([combs[0],combs[2]]):\n",
    "    p = comb[0]\n",
    "    m = comb[1]\n",
    "    ind = np.where(np.asarray(efficiency_global[combs[i*2][0]][combs[i*2][1]][\"x\"][\"min\"]) > np.asarray(efficiency_global[combs[i*2][0]][combs[i*2][1]][\"x\"][\"max\"]))[0]\n",
    "    print(str(len(ind)) + \"/\" + str(numcitiestotal) + \" cities found with x_min>x_max for \"+ str(combs[i]))\n",
    "    print([list(cities.keys())[i] for i in ind])\n",
    "    \n",
    "    \n",
    "    for j in range(len(ind)):\n",
    "        ax = axes[j + i*len(ind)]\n",
    "        placeid = list(cities.keys())[ind[j]]\n",
    "        tmp, = ax.plot(prune_quantiles, analysis_result[p][m][placeid][\"efficiency_global_routed\"])\n",
    "        tmp.set_label('_hidden')\n",
    "        ax.set_title(placeid)\n",
    "        ax.set_ylim([0,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot only the biggest dips"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cities.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bigdips = [\"delft\", \"copenhagen\", \"boston\", \"paris\", \"barcelona\", \"sheffield\"]\n",
    "\n",
    "fig, axes = plt.subplots(nrows = 2, ncols = 3, figsize=(300/plotparam[\"dpi\"], 220/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"], squeeze = True)\n",
    "axes = axes.flatten()\n",
    "\n",
    "i = 0\n",
    "comb = combs[0]\n",
    "p = comb[0]\n",
    "m = comb[1]\n",
    "\n",
    "for j,placeid in enumerate(bigdips):  \n",
    "    ax = axes[j]\n",
    "    tmp, = ax.plot(prune_quantiles, analysis_result[p][m][placeid][\"efficiency_global_routed\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "    ax.text(0.96, 0.08, cities[placeid][\"name\"], fontsize=8, horizontalalignment='right')\n",
    "    ax.set_xlim(0, 1)\n",
    "    ax.set_ylim(0, 1)\n",
    "    ax.set_xticks([0, 0.5, 1])\n",
    "    if j==4:\n",
    "        ax.set_xlabel('Betweenness quantile')\n",
    "    if j <= 2:\n",
    "        ax.set_xlabel('')\n",
    "        ax.set_xticklabels([])\n",
    "    else:\n",
    "        ax.set_xticklabels([0,0.5,1])\n",
    "    ax.set_yticks([0,0.25,0.5,0.75,1])\n",
    "    if j % 3 == 0:\n",
    "        ax.set_yticklabels([0,0.25,0.5,0.75,1])\n",
    "    else:\n",
    "        ax.set_ylabel('')\n",
    "        ax.set_yticklabels([])\n",
    "    if j == 1:\n",
    "        ax.set_title('Dips in global efficiency')\n",
    "plt.subplots_adjust(top = 0.90, bottom = 0.16, left = 0.12, right = 0.97, wspace = 0.22, hspace = 0.22)\n",
    "fig.savefig(PATH[\"plots\"] + \"/\" + 'dipsglobalefficiency_poi_' + p + '.eps', facecolor = \"white\", edgecolor = 'none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now figure for paper: Directness and Global eff, for Boston (big left) and 4 cities (small right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.font_manager import findfont, FontProperties\n",
    "font = findfont(FontProperties(family=['sans-serif']))\n",
    "font"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smallmultiples = [\"montreal\", \"mumbai\", \"paris\", \"tokyo\"]\n",
    "checkpoints = [0.025, 0.1, 0.2]\n",
    "checkpoints_x = [0,3,7]\n",
    "xticks = [0, 0.1, 0.2,0.3,0.4]\n",
    "yticks_top = [0.6,0.7,0.8,0.9]\n",
    "yticks_bot = [0,0.2,0.4,0.6,0.8,1]\n",
    "ytextpos_top = [0.87, 0.835]\n",
    "ytextpos_bot = [0.04, 0.08]\n",
    "\n",
    "fig, axes = plt.subplots(nrows = 4, ncols = 4, figsize=(300/plotparam[\"dpi\"], 300/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"], squeeze = True)\n",
    "\n",
    "i = 0\n",
    "comb = combs[0]\n",
    "p = comb[0]\n",
    "m = comb[1]\n",
    "smallmultipleindex = 0\n",
    "\n",
    "for row in [1,2,3,4]:\n",
    "    if row < 3:\n",
    "        yticks = yticks_top\n",
    "        ytextpos = ytextpos_top\n",
    "    else:\n",
    "        yticks = yticks_bot\n",
    "        ytextpos = ytextpos_bot\n",
    "    for col in [1,2,3,4]:\n",
    "        i += 1\n",
    "        placeid = \"boston\"\n",
    "        if col >= 3: # Right half: small multiple\n",
    "            ax = plt.subplot(4,4,i)\n",
    "            placeid = smallmultiples[smallmultipleindex]\n",
    "            smallmultipleindex += 1\n",
    "            smallmultipleindex = smallmultipleindex % len(smallmultiples)\n",
    "            ax.text(0.96*xticks[-1], ytextpos[1], cities[placeid][\"name\"], fontsize=8, horizontalalignment='right')\n",
    "            if row <= 2:\n",
    "                ax.plot(prune_quantiles, analysis_result[p][m][placeid][\"directness_all_linkwise\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "            else:\n",
    "                ax.plot(prune_quantiles, analysis_result[p][m][placeid][\"efficiency_global_routed\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "            ax.set_yticklabels([])\n",
    "            ax.set_yticks(yticks)\n",
    "            ax.set_xticks(xticks)\n",
    "            if row == 4:\n",
    "                ax.set_xticklabels('')\n",
    "            else:\n",
    "                ax.set_xticklabels('')\n",
    "            \n",
    "        elif (row == 1 and col == 1): # Top left\n",
    "            ax = plt.subplot(4,4,(1,6))\n",
    "            ax.plot(prune_quantiles, analysis_result[p][m][placeid][\"directness_all_linkwise\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "            ax.text(0.98*xticks[-1], ytextpos[0], cities[placeid][\"name\"], fontsize=8, horizontalalignment='right')\n",
    "            ax.set_yticks(yticks)\n",
    "            ax.set_yticklabels(yticks)\n",
    "            ax.set_ylabel('Directness')\n",
    "            ax.set_xticks(xticks)\n",
    "            ax.set_xticklabels([])\n",
    "            for j,c in enumerate(checkpoints):\n",
    "                ax.plot([c, c], [0, 1], \":k\", linewidth = 1)\n",
    "                ax.plot([c], [analysis_result[p][m][placeid][\"directness_all_linkwise\"][checkpoints_x[j]]], \"ok\", fillstyle='none', markersize = 5)\n",
    "            \n",
    "        elif (row == 3 and col == 1): # Bottom left\n",
    "            ax = plt.subplot(4,4,(9,14))\n",
    "            ax.text(0.98*xticks[-1], ytextpos[0], cities[placeid][\"name\"], fontsize=8, horizontalalignment='right')\n",
    "            ax.plot(prune_quantiles, analysis_result[p][m][placeid][\"efficiency_global_routed\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "            ax.set_yticks(yticks)\n",
    "            ax.set_yticklabels(yticks)\n",
    "            ax.set_ylabel('Global efficiency')\n",
    "            ax.set_xlabel('Betweenness quantile $q_B$')\n",
    "            ax.set_xticks(xticks)\n",
    "            ax.set_xticklabels(xticks)\n",
    "            for j,c in enumerate(checkpoints):\n",
    "                ax.plot([c, c], [0, 1], \":k\", linewidth = 1)\n",
    "                ax.plot([c], [analysis_result[p][m][placeid][\"efficiency_global_routed\"][checkpoints_x[j]]], \"ok\", fillstyle='none', markersize = 5)\n",
    "            \n",
    "        ax.set_xlim([0, 0.4])\n",
    "        if row < 3:\n",
    "            ax.set_ylim(0.6, 0.9)\n",
    "        else:\n",
    "            ax.set_ylim(0, 1)\n",
    "        \n",
    "plt.subplots_adjust(top = 0.98, bottom = 0.13, left = 0.14, right = 0.97, wspace = 0.22, hspace = 0.22)\n",
    "            \n",
    "            \n",
    "\n",
    "fig.savefig(PATH[\"plots\"] + \"/\" + 'smallmultiples_poi_' + p + '.eps', facecolor = \"white\", edgecolor = 'none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing grown with existing metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we plot for some metrics how much higher the metric of the grown network is compared to the existing network of same length. This only works for cities that have a small enough existing length which is at some point reached by the grown network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comp_keys = [\"length_lcc\", \"efficiency_global_routed\", \"efficiency_local_routed\", \"coverage\"]\n",
    "comp_labels = [\"Length of LCC\", \"Global Efficiency\", \"Local Efficiency\", \"Coverage\"]\n",
    "\n",
    "cities_here = []\n",
    "\n",
    "for p in poi_source_list:\n",
    "    for m in prune_measure_list:\n",
    "        numcities = 0\n",
    "        values = []\n",
    "\n",
    "        if debug: fig = plt.figure(figsize=(400/plotparam[\"dpi\"], 400/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"])\n",
    "        if debug: plt.semilogy([min(x), max(x)], [1,1], \"k--\")\n",
    "        for placeid, placeinfo in tqdm(cities.items(), desc=\"Cities\"):\n",
    "            length_existing = analysis_existing[placeid][analysis_existing_rowkeys[\"biketrack\"]][\"length\"]\n",
    "            if length_existing and np.argmax(analysis_result[p][m][placeid][\"length\"] > length_existing):\n",
    "                numcities += 1\n",
    "                id_samelen = np.argmax(analysis_result[p][m][placeid][\"length\"] > length_existing) - 1\n",
    "\n",
    "                x = list(range(len(comp_keys)))\n",
    "                y = [[analysis_result[p][m][placeid][comp_keys[i]][id_samelen] / analysis_existing[placeid][analysis_existing_rowkeys[\"biketrack\"]][comp_keys[i]]] for i in x]\n",
    "                values.append(y)\n",
    "                cities_here.append(placeid)\n",
    "                if debug: plt.semilogy(x, y, \"o\")\n",
    "                print(id_samelen)\n",
    "\n",
    "        print(str(numcities) + \" cities found where L_grown=L for \" + p + \" | \" + m)\n",
    "        print(cities_here)\n",
    "        values = np.log10(np.array(values))\n",
    "        fig = plt.figure(figsize=(200/plotparam[\"dpi\"], 100/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"])\n",
    "        axes = fig.add_axes([0, 0, 1, 1])\n",
    "        axes.plot([0,0], [min(x), max(x)+2],\"k:\", linewidth=0.5)\n",
    "        # https://stackoverflow.com/questions/18500011/horizontal-box-plots-in-matplotlib-pandas#56088231\n",
    "        axes.boxplot([values[np.isfinite(values[:,i]).flatten(),i].flatten().tolist() for i in x], vert=False, showfliers=False); #with isfinite we exclude zeros, infs, or nans which sometimes happen for e.g. local efficiency\n",
    "        axes.set_xlabel('$M_{syn}/M_{real}$ | $L_{syn}=L_{real}$')\n",
    "        axes.set_title('Synthetic versus real networks')\n",
    "        axes.set_yticklabels([comp_labels[i] for i in x]);\n",
    "        axes.set_ylim([min(x)+0.5, max(x)+1.5])\n",
    "        # fig.autofmt_xdate(rotation=45)\n",
    "\n",
    "        axes.set_xlim([-0.7, 1.55])\n",
    "#         axes.set_xlim([-0.35, 1.55])\n",
    "        axes.set_xticks([-0.301,0,0.301, 0.699, 1, 1.301])\n",
    "        axes.set_xticklabels([0.5, 1, 2,5, 10, 20])\n",
    "        axes.xaxis.set_minor_locator(matplotlib.ticker.FixedLocator(np.log10(np.concatenate((np.linspace(0,1,10, endpoint=False),np.linspace(1, 10,9, endpoint=False),np.linspace(10, 100,9, endpoint=False))))))\n",
    "        axes.text(1.5, 4, p + \" | \" + m, fontsize=8, horizontalalignment='right');\n",
    "        fig.savefig(PATH[\"plots\"] + 'grownvsexisting_poi_' + p + \"_\" + m + '.eps', facecolor = \"white\", edgecolor = 'none', bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Existing cover plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do this for Milan."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Real Milan\n",
    "\n",
    "# %%capture\n",
    "    \n",
    "for placeid, placeinfo in tqdm(cities.items(), desc = \"Cities\"):\n",
    "    print(placeid + \": Plotting network covers\")\n",
    "\n",
    "    # Load networks\n",
    "    G_biketrack = csv_to_ig(PATH[\"data\"] + placeid + \"/\", placeid, 'biketrack')\n",
    "    G_carall = csv_to_ig(PATH[\"data\"] + placeid + \"/\", placeid, 'carall')\n",
    "    G_biketrackcarall = csv_to_ig(PATH[\"data\"] + placeid + \"/\", placeid, 'biketrackcarall')\n",
    "    G_bikeable = csv_to_ig(PATH[\"data\"] + placeid + \"/\", placeid, 'bikeable')\n",
    "    map_center = nxdraw(G_carall, \"carall\")\n",
    "    \n",
    "    \n",
    "    # Load results\n",
    "    filename = placeid + '_poi_' + poi_source + \"_\" + prune_measure + \".pickle\"\n",
    "    with open(PATH[\"results\"] + placeid + \"/\" + filename, 'rb') as f:\n",
    "        res = pickle.load(f)\n",
    "    \n",
    "    # Load covers\n",
    "    filename = placeid + \"_\"  + \"existing_covers\"\n",
    "    with open(PATH[\"results\"] + placeid + \"/\" + filename + \".pickle\",'rb') as f:\n",
    "        cov_car = pickle.load(f)['carall']\n",
    "    with open(PATH[\"results\"] + placeid + \"/\" + filename + \".pickle\",'rb') as f:\n",
    "        cov_biketrack_onstreet = pickle.load(f)['biketrack_onstreet']\n",
    "    with open(PATH[\"results\"] + placeid + \"/\" + filename + \".pickle\",'rb') as f:\n",
    "        cov_bikeable_offstreet = pickle.load(f)['bikeable_offstreet']\n",
    "    # Merge on and offstreet\n",
    "    cov_bike = ops.unary_union([cov_biketrack_onstreet, cov_bikeable_offstreet])\n",
    "    \n",
    "    # Construct and plot patches from covers\n",
    "    patchlist_car, patchlist_car_holes = cov_to_patchlist(cov_car, map_center)\n",
    "    fig = initplot()\n",
    "\n",
    "    # Covers\n",
    "    axes = fig.add_axes([0, 0, 1, 1]) # left, bottom, width, height (range 0 to 1)\n",
    "    patchlist_bike, patchlist_bike_holes = cov_to_patchlist(cov_bike, map_center)\n",
    "\n",
    "    # We have this contrived order due to alphas, holes, and matplotlib's inability to draw polygon patches with holes. This only works because the car network is a superset of the bike network.\n",
    "    # car orange, bike white, bike blue, bike holes white, bike holes orange, car holes white\n",
    "    patchlist_combined = patchlist_car + patchlist_bike + patchlist_bike + patchlist_bike_holes+ patchlist_bike_holes + patchlist_car_holes\n",
    "    pc = PatchCollection(patchlist_combined)\n",
    "    colors = np.array([[255/255,115/255,56/255,0.2] for _ in range(len(patchlist_car))]) # car orange\n",
    "    if len(patchlist_bike):\n",
    "        colors = np.append(colors, [[1,1,1,1] for _ in range(len(patchlist_bike))], axis = 0) # bike white\n",
    "        colors = np.append(colors, [[86/255,220/255,244/255,0.4] for _ in range(len(patchlist_bike))], axis = 0) # bike blue\n",
    "    if len(patchlist_bike_holes):\n",
    "        colors = np.append(colors, [[1,1,1,1] for _ in range(len(patchlist_bike_holes))], axis = 0) # bike holes white\n",
    "    if len(patchlist_bike_holes):\n",
    "        colors = np.append(colors, [[255/255,115/255,56/255,0.2] for _ in range(len(patchlist_bike_holes))], axis = 0) # bike holes orange\n",
    "    if len(patchlist_car_holes):\n",
    "        colors = np.append(colors, [[1,1,1,1] for _ in range(len(patchlist_car_holes))], axis = 0) # car holes white\n",
    "    pc.set_facecolors(colors)\n",
    "    pc.set_edgecolors(np.array([[0,0,0,0.4] for _ in range(len(patchlist_combined))])) # remove this line if the outline of the full coverage should remain\n",
    "    axes.add_collection(pc)\n",
    "    axes.set_aspect('equal')\n",
    "    axes.set_xmargin(0.01)\n",
    "    axes.set_ymargin(0.01)\n",
    "    axes.plot()\n",
    "\n",
    "    # Networks\n",
    "    nxdraw(G_carall, \"carall\", map_center)\n",
    "    nxdraw(G_biketrack, \"bikegrown\", map_center, nodesize = nodesize_grown)\n",
    "    nxdraw(G_biketrack.clusters().giant(), \"highlight_biketrack\", map_center, nodesize = nodesize_grown)\n",
    "    plt.savefig(PATH[\"plots_networks\"] + placeid + \"/\" + placeid + '_bikecarcover.png', bbox_inches=\"tight\", dpi=plotparam[\"dpi\"])\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Synthetic Milan\n",
    "\n",
    "# %%capture\n",
    "for placeid, placeinfo in tqdm(cities.items(), desc = \"Cities\"):\n",
    "    print(placeid + \": Plotting network covers\")\n",
    "\n",
    "    # Load networks\n",
    "    G_biketrack = csv_to_ig(PATH[\"data\"] + placeid + \"/\", placeid, 'biketrack')\n",
    "    G_carall = csv_to_ig(PATH[\"data\"] + placeid + \"/\", placeid, 'carall')\n",
    "    G_biketrackcarall = csv_to_ig(PATH[\"data\"] + placeid + \"/\", placeid, 'biketrackcarall')\n",
    "    G_bikeable = csv_to_ig(PATH[\"data\"] + placeid + \"/\", placeid, 'bikeable')\n",
    "    map_center = nxdraw(G_carall, \"carall\")\n",
    "    \n",
    "    # Load POIs\n",
    "    with open(PATH[\"data\"] + placeid + \"/\" + placeid + '_poi_' + poi_source + '_nnidscarall.csv') as f:\n",
    "        nnids = [int(line.rstrip()) for line in f]\n",
    "    nodesize_poi = nodesize_from_pois(nnids)\n",
    "    \n",
    "    # Load results\n",
    "    filename = placeid + '_poi_' + poi_source + \"_\" + prune_measure + \".pickle\"\n",
    "    with open(PATH[\"results\"] + placeid + \"/\" + filename, 'rb') as f:\n",
    "        res = pickle.load(f)\n",
    "    \n",
    "    # Load covers\n",
    "    filename = placeid + '_poi_' + poi_source + \"_\" + prune_measure + \"_covers\"\n",
    "    with open(PATH[\"results\"] + placeid + \"/\" + filename + \".pickle\",'rb') as f:\n",
    "        covs = pickle.load(f)\n",
    "    filename = placeid + \"_\"  + \"existing_covers\"\n",
    "    with open(PATH[\"results\"] + placeid + \"/\" + filename + \".pickle\",'rb') as f:\n",
    "        cov_car = pickle.load(f)['carall']\n",
    "    \n",
    "    # Construct and plot patches from covers\n",
    "    patchlist_car, patchlist_car_holes = cov_to_patchlist(cov_car, map_center)\n",
    "    GT = res[\"GTs\"][16]\n",
    "    prune_quantile = res[\"prune_quantiles\"][16]\n",
    "    cov = list(covs.values())[16]\n",
    "    fig = initplot()\n",
    "\n",
    "    # Covers\n",
    "    axes = fig.add_axes([0, 0, 1, 1]) # left, bottom, width, height (range 0 to 1)\n",
    "    patchlist_bike, patchlist_bike_holes = cov_to_patchlist(cov, map_center)\n",
    "\n",
    "    # We have this contrived order due to alphas, holes, and matplotlib's inability to draw polygon patches with holes. This only works because the car network is a superset of the bike network.\n",
    "    # car orange, bike white, bike blue, bike holes white, bike holes orange, car holes white\n",
    "    patchlist_combined = patchlist_car + patchlist_bike + patchlist_bike + patchlist_bike_holes+ patchlist_bike_holes + patchlist_car_holes\n",
    "    pc = PatchCollection(patchlist_combined)\n",
    "    colors = np.array([[255/255,115/255,56/255,0.2] for _ in range(len(patchlist_car))]) # car orange\n",
    "    if len(patchlist_bike):\n",
    "        colors = np.append(colors, [[1,1,1,1] for _ in range(len(patchlist_bike))], axis = 0) # bike white\n",
    "        colors = np.append(colors, [[86/255,220/255,244/255,0.4] for _ in range(len(patchlist_bike))], axis = 0) # bike blue\n",
    "    if len(patchlist_bike_holes):\n",
    "        colors = np.append(colors, [[1,1,1,1] for _ in range(len(patchlist_bike_holes))], axis = 0) # bike holes white\n",
    "    if len(patchlist_bike_holes):\n",
    "        colors = np.append(colors, [[255/255,115/255,56/255,0.2] for _ in range(len(patchlist_bike_holes))], axis = 0) # bike holes orange\n",
    "    if len(patchlist_car_holes):\n",
    "        colors = np.append(colors, [[1,1,1,1] for _ in range(len(patchlist_car_holes))], axis = 0) # car holes white\n",
    "    pc.set_facecolors(colors)\n",
    "    pc.set_edgecolors(np.array([[0,0,0,0.4] for _ in range(len(patchlist_combined))])) # remove this line if the outline of the full coverage should remain\n",
    "    axes.add_collection(pc)\n",
    "    axes.set_aspect('equal')\n",
    "    axes.set_xmargin(0.01)\n",
    "    axes.set_ymargin(0.01)\n",
    "    axes.plot()\n",
    "\n",
    "    # Networks\n",
    "    nxdraw(G_carall, \"carall\", map_center)\n",
    "    nxdraw(GT, \"highlight_biketrack\", map_center, nodesize = nodesize_grown)\n",
    "    plt.savefig(PATH[\"plots_networks\"] + placeid + \"/\" + placeid + '_bikesynthcarcover.png', bbox_inches=\"tight\", dpi=plotparam[\"dpi\"])\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Average analysis plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we plot one analysis figure as an average over all cities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Each parameter set as a single figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%capture\n",
    "\n",
    "for poi_source, prune_measure in combs:\n",
    "    print(poi_source, prune_measure)\n",
    "        \n",
    "    analysis_result_city = {}\n",
    "    numcities = 0\n",
    "    for placeid, placeinfo in tqdm(cities.items(), desc=\"Cities\"):\n",
    "\n",
    "        # PLOT Analysis\n",
    "        filename = placeid + '_poi_' + poi_source + \"_\" + prune_measure + \".csv\"\n",
    "        analysis_result_city_temp = np.genfromtxt(PATH[\"results\"] + placeid + \"/\" + filename, delimiter=',', names=True)\n",
    "        if len(analysis_result_city_temp) == 0: # Discard if no results (for example no railwaystations)\n",
    "            print(placeid + \": No analysis results available\")\n",
    "            continue\n",
    "        else:\n",
    "            numcities += 1\n",
    "            analysis_result_city[placeid] = analysis_result_city_temp\n",
    "            metric_keys = analysis_result_city[placeid].dtype.names\n",
    "    \n",
    "    # All cities are loaded, now create the average values\n",
    "    analysis_result_AVG = {}\n",
    "    for metric in metric_keys:\n",
    "        temp = np.zeros([40, numcities])\n",
    "        for i, placeid_this in enumerate(analysis_result_city.keys()):\n",
    "            temp[:, i] = analysis_result_city[placeid_this][metric]\n",
    "        analysis_result_AVG[metric] = np.mean(temp, axis = 1)\n",
    "            \n",
    "    \n",
    "                    \n",
    "    # Plot\n",
    "    nc = 5\n",
    "    fig, axes = plt.subplots(nrows = 2, ncols = nc, figsize = (16, 6))\n",
    "    keys_metrics = {\"length\": \"Length [km]\",\"coverage\": \"Coverage [km$^2$]\",\"overlap_biketrack\": \"Overlap Protected\",\"directness\": \"Directness\",\"efficiency_global\": \"Global Efficiency\",\n",
    "            \"length_lcc\": \"Length of LCC [km]\",\"poi_coverage\": \"Seed points covered\",\"overlap_bikeable\": \"Overlap Bikeable\",\"components\": \"Components\",\"efficiency_local\": \"Local Efficiency\"}\n",
    "\n",
    "    for i, ax in enumerate(axes[0]):\n",
    "        key = list(keys_metrics.keys())[i]\n",
    "        if key in [\"overlap_biketrack\", \"overlap_bikeable\"]:\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[key] / analysis_result_AVG[\"length\"], **plotparam_analysis[\"bikegrown\"])\n",
    "        elif key in [\"efficiency_global\", \"efficiency_local\"]:\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[key], **plotparam_analysis[\"bikegrown_abstract\"])\n",
    "            tmp, = ax.plot(prune_quantiles, analysis_result_AVG[key+\"_routed\"], **plotparam_analysis[\"bikegrown\"])\n",
    "            tmp.set_label('_hidden')\n",
    "        elif key in [\"length\", \"length_lcc\"]: # Convert m->km\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[key]/1000, **plotparam_analysis[\"bikegrown\"])\n",
    "        else:\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[key], **plotparam_analysis[\"bikegrown\"])\n",
    "\n",
    "        if i == 0:\n",
    "            ymax0 = ax.get_ylim()[1]\n",
    "            ax.set_ylim(0, ymax0)\n",
    "            ax.text(-0.15, ymax0*1.25, \"Average city\" + \" (\" + poi_source + \" | \" + prune_measure + \")\", fontsize=16, horizontalalignment='left')\n",
    "            ax.legend(loc = 'upper left')\n",
    "        if i == 4:\n",
    "            ax.legend(loc = 'best')\n",
    "\n",
    "\n",
    "        set_analysissubplot(key)\n",
    "        ax.set_title(list(keys_metrics.values())[i])\n",
    "        ax.set_xlabel('')\n",
    "        ax.set_xticklabels([])\n",
    "\n",
    "\n",
    "    for i, ax in enumerate(axes[1]):\n",
    "        key = list(keys_metrics.keys())[i+nc]\n",
    "        if key in [\"overlap_biketrack\", \"overlap_bikeable\"]:\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[key] / analysis_result_AVG[\"length\"], **plotparam_analysis[\"bikegrown\"])\n",
    "        elif key in [\"efficiency_global\", \"efficiency_local\"]:\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[key], **plotparam_analysis[\"bikegrown_abstract\"])\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[key+\"_routed\"], **plotparam_analysis[\"bikegrown\"])\n",
    "        elif key in [\"length\", \"length_lcc\"]: # Convert m->km\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[key]/1000, **plotparam_analysis[\"bikegrown\"])\n",
    "        else:\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[key], **plotparam_analysis[\"bikegrown\"])\n",
    "\n",
    "        if i == 0:\n",
    "            ax.set_ylim(0, ymax0)\n",
    "        set_analysissubplot(key)\n",
    "        ax.set_title(list(keys_metrics.values())[i+nc])\n",
    "        ax.set_xlabel(prune_measure + ' quantile')\n",
    "        if key in [\"poi_coverage\"]:\n",
    "            # https://stackoverflow.com/questions/30914462/matplotlib-how-to-force-integer-tick-labels\n",
    "            ax.yaxis.set_major_locator(MaxNLocator(integer=True)) \n",
    "\n",
    "    plt.subplots_adjust(top = 0.87, bottom = 0.09, left = 0.05, right = 0.97, wspace = 0.25, hspace = 0.25)\n",
    "    if plotconstricted:\n",
    "        fig.savefig(PATH[\"plots\"] + \"/\" + 'averagecity_analysis_poi_' + poi_source + \"_\" + prune_measure + '.png', facecolor = \"white\", edgecolor = 'none')\n",
    "    else:\n",
    "        fig.savefig(PATH[\"plots\"] + \"/\" + 'averagecity_analysis_poi_' + poi_source + \"_\" + prune_measure + '_noconstr.png', facecolor = \"white\", edgecolor = 'none')\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figures with bundled parameter sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_analysissubplot_special(key, p):\n",
    "    if key in [\"overlap_bikeable\"]:\n",
    "        ax.set_ylim(0.08,0.22)\n",
    "    if key in [\"overlap_biketrack\"]:\n",
    "        ax.set_ylim(0.015,0.065)\n",
    "    \n",
    "    if p == \"grid\":\n",
    "        if key in [\"components\"]:\n",
    "            ax.set_ylim(top = 25)\n",
    "        if key in [\"efficiency_local\"]:\n",
    "            ax.set_ylim(top = 0.25)\n",
    "#         if key in [\"overlap_bikeable\"]:\n",
    "#             ax.set_ylim(top = 0.25)\n",
    "#         if key in [\"overlap_biketrack\"]:\n",
    "#             ax.set_ylim(top = 0.06)\n",
    "        if key in [\"directness\"]:\n",
    "            ax.set_ylim(bottom = 0.43)\n",
    "            ax.set_ylim(top = 0.81)\n",
    "        if key in [\"directness_all_linkwise\"]:\n",
    "            ax.set_ylim(bottom = 0.47)\n",
    "            ax.set_ylim(top = 0.84)\n",
    "    elif p == \"railwaystation\":\n",
    "        if key in [\"components\"]:\n",
    "            ax.set_ylim(top = 25)\n",
    "        if key in [\"efficiency_local\"]:\n",
    "            ax.set_ylim(top = 0.25)\n",
    "#         if key in [\"overlap_bikeable\"]:\n",
    "#             ax.set_ylim(top = 0.25)\n",
    "#         if key in [\"overlap_biketrack\"]:\n",
    "#             ax.set_ylim(top = 0.06)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%capture\n",
    "plotparam_analysis[\"constricted\"] = {\"linewidth\": 1.5, \"color\": '#D22A0E', \"linestyle\": \"solid\", \"label\": \"Betw. cars\"}\n",
    "plotparam_analysis[\"mst\"] = {\"linewidth\": 1, \"color\": '#aaaaaa', \"linestyle\": \"dotted\", \"label\": \"MST\"}\n",
    "plotparam_analysis[\"bikegrown_betweenness\"] = {\"linewidth\": 2.5, \"color\": '#0eb6d2', \"linestyle\": \"solid\", \"label\": \"Betw.\"}\n",
    "plt.style.use(PATH[\"parameters\"] + 'plotstyle.mplstyle')\n",
    "# Run all parameter sets\n",
    "prune_measure_list = [\"betweenness\", \"closeness\", \"random\"]\n",
    "poi_source = \"grid\" # railwaystation grid\n",
    "\n",
    "analysis_result_city = {}\n",
    "analysis_result_AVG = {}\n",
    "analysis_constricted_AVG = {}\n",
    "analysis_mst_result_city = {}\n",
    "for m in prune_measure_list:\n",
    "    analysis_result_city[m] = {}\n",
    "    numcities = 0\n",
    "    for placeid, placeinfo in tqdm(cities.items(), desc=\"Cities\"):\n",
    "\n",
    "        # PLOT Analysis\n",
    "        filename = placeid + '_poi_' + poi_source + \"_mst.csv\"\n",
    "        analysis_mst_result_city_temp = np.genfromtxt(PATH[\"results\"] + placeid + \"/\" + filename, delimiter=',', names=True)\n",
    "        \n",
    "        filename = placeid + '_poi_' + poi_source + \"_\" + m + \".csv\"\n",
    "        analysis_result_city_temp = np.genfromtxt(PATH[\"results\"] + placeid + \"/\" + filename, delimiter=',', names=True)\n",
    "        if len(analysis_result_city_temp) == 0: # Discard if no results (for example no railwaystations)\n",
    "            print(placeid + \": No analysis results available\")\n",
    "            continue\n",
    "        else:\n",
    "            numcities += 1\n",
    "            analysis_result_city[m][placeid] = analysis_result_city_temp\n",
    "            metric_keys = analysis_result_city[m][placeid].dtype.names\n",
    "            analysis_mst_result_city[placeid] = analysis_mst_result_city_temp\n",
    "\n",
    "    # All cities are loaded, now create the average values\n",
    "    analysis_result_AVG[m] = {}\n",
    "    for metric in metric_keys:\n",
    "        temp = np.zeros([40, numcities])\n",
    "        for i, placeid_this in enumerate(analysis_result_city[m].keys()):\n",
    "            temp[:, i] = analysis_result_city[m][placeid_this][metric]\n",
    "        analysis_result_AVG[m][metric] = np.mean(temp, axis = 1)\n",
    "        \n",
    "    analysis_mst_result_AVG = {}\n",
    "    for metric in metric_keys:\n",
    "        temp = np.zeros([1, numcities])\n",
    "        for i, placeid_this in enumerate(analysis_mst_result_city.keys()):\n",
    "            temp[:, i] = analysis_mst_result_city[placeid_this][metric]\n",
    "        analysis_mst_result_AVG[metric] = np.mean(temp, axis = 1)\n",
    "    print(analysis_mst_result_AVG)\n",
    "\n",
    "    analysis_constricted_AVG[m] = {}\n",
    "    for metric in list(range(9)): # 9 constricted fields are loaded\n",
    "        temp = np.zeros([40, len(analysis_constricted[poi_source][m].keys())]) \n",
    "        for i, placeid_this in enumerate(analysis_constricted[poi_source][m].keys()):\n",
    "            temp[:, i] = analysis_constricted[poi_source][m][placeid_this][:, metric]\n",
    "        analysis_constricted_AVG[m][metric] = np.mean(temp, axis = 1)\n",
    "\n",
    "# Plot\n",
    "nc = 4\n",
    "fig, axes = plt.subplots(nrows = 2, ncols = nc, figsize=(640/plotparam[\"dpi\"], 320/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"])\n",
    "keys_metrics = {\"length\": \"Length [km]\",\"coverage\": \"Coverage [km$^2$]\", \"directness_all_linkwise\": \"Directness\", \"efficiency_global\": \"Global Efficiency\",\n",
    "        \"length_lcc\": \"Length of LCC [km]\",\"poi_coverage\": \"Seed points covered\",\"components\": \"Components\",\"efficiency_local\": \"Local Efficiency\"}\n",
    "\n",
    "for i, ax in enumerate(axes[0]):\n",
    "    key = list(keys_metrics.keys())[i]\n",
    "    for m in prune_measure_list:\n",
    "        if key in [\"overlap_biketrack\", \"overlap_bikeable\"]:\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[m][key] / analysis_result_AVG[m][\"length\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "        elif key in [\"efficiency_global\", \"efficiency_local\"]:\n",
    "            if m == \"betweenness\":\n",
    "                tmp, = ax.plot([0,1], [analysis_mst_result_AVG[key], analysis_mst_result_AVG[key]], **plotparam_analysis[\"mst\"])\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[m][key+\"_routed\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "        elif key in [\"length\", \"length_lcc\"]: # Convert m->km\n",
    "            if m == \"betweenness\":\n",
    "                ax.plot([0,1], [analysis_mst_result_AVG[key]/1000, analysis_mst_result_AVG[key]/1000], **plotparam_analysis[\"mst\"])\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[m][key]/1000, **plotparam_analysis[\"bikegrown_\" + m])\n",
    "        else:\n",
    "            if m == \"betweenness\":\n",
    "                ax.plot([0,1], [analysis_mst_result_AVG[key], analysis_mst_result_AVG[key]], **plotparam_analysis[\"mst\"])\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[m][key], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "            \n",
    "        if key == \"directness_all_linkwise\" and plotconstricted and m == \"betweenness\": # Only show betweenness - the others are similar\n",
    "            tmp, = ax.plot(prune_quantiles, analysis_constricted_AVG[m][8], **plotparam_analysis[\"constricted\"])\n",
    "            tmp.set_label('_hidden')\n",
    "        if key == \"efficiency_global\" and plotconstricted and m == \"betweenness\": # Only show betweenness - the others are similar\n",
    "            ax.plot(prune_quantiles, analysis_constricted_AVG[m][0], **plotparam_analysis[\"constricted\"])\n",
    "                \n",
    "        if key in [\"directness_all_linkwise\"]:\n",
    "            ax.arrow(0.33, 0.62, 0, -0.04, color = \"orange\", width = 0.01, head_width = 0.05, head_length = 0.02)\n",
    "\n",
    "        if i == 0:\n",
    "            ymax0 = ax.get_ylim()[1]\n",
    "            ax.set_ylim(0, ymax0)\n",
    "            if not plotnotitle:\n",
    "                ax.text(-0.15, ymax0*1.25, \"Average city\" + \" (\" + poi_source + \")\", fontsize=16, horizontalalignment='left')\n",
    "                \n",
    "        if i == 3:\n",
    "            leg = ax.legend(loc = 'best')\n",
    "            leg.get_frame().set_linewidth(0.5)\n",
    "            ax.set_zorder(100)\n",
    "            # https://stackoverflow.com/questions/23238041/move-and-resize-legends-box-in-matplotlib\n",
    "            bb = leg.get_bbox_to_anchor().inverse_transformed(ax.transAxes)\n",
    "            # Change to location of the legend. \n",
    "            xOffset = -1.26\n",
    "            bb.x0 += xOffset\n",
    "            bb.x1 += xOffset\n",
    "            yOffset = -1.26\n",
    "            bb.y0 += yOffset\n",
    "            bb.y1 += yOffset\n",
    "            leg.set_bbox_to_anchor(bb, transform = ax.transAxes)\n",
    "            \n",
    "            \n",
    "        \n",
    "\n",
    "    set_analysissubplot(key)\n",
    "    set_analysissubplot_special(key, poi_source)\n",
    "    ax.set_title(list(keys_metrics.values())[i])\n",
    "    ax.set_xlabel('')\n",
    "    ax.set_xticklabels([])\n",
    "\n",
    "\n",
    "for i, ax in enumerate(axes[1]):\n",
    "    key = list(keys_metrics.keys())[i+nc]\n",
    "    for m in prune_measure_list:\n",
    "        if key in [\"overlap_biketrack\", \"overlap_bikeable\"]:\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[m][key] / analysis_result_AVG[m][\"length\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "        elif key in [\"efficiency_global\", \"efficiency_local\"]:\n",
    "            if m == \"betweenness\":\n",
    "                tmp, = ax.plot([0,1], [analysis_mst_result_AVG[key], analysis_mst_result_AVG[key]], **plotparam_analysis[\"mst\"])\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[m][key+\"_routed\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "        elif key in [\"length\", \"length_lcc\"]: # Convert m->km\n",
    "            if m == \"betweenness\":\n",
    "                tmp, = ax.plot([0,1], [analysis_mst_result_AVG[key]/1000, analysis_mst_result_AVG[key]/1000], **plotparam_analysis[\"mst\"])\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[m][key]/1000, **plotparam_analysis[\"bikegrown_\" + m])\n",
    "        else:\n",
    "            if m == \"betweenness\":\n",
    "                tmp, = ax.plot([0,1], [analysis_mst_result_AVG[key], analysis_mst_result_AVG[key]], **plotparam_analysis[\"mst\"])\n",
    "            ax.plot(prune_quantiles, analysis_result_AVG[m][key], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "        \n",
    "        if key == \"efficiency_local\" and plotconstricted and m == \"betweenness\": # Only show betweenness - the others are similar\n",
    "            tmp, = ax.plot(prune_quantiles, analysis_constricted_AVG[m][1], **plotparam_analysis[\"constricted\"])\n",
    "            tmp.set_label('_hidden')\n",
    "                \n",
    "        if i == 0:\n",
    "            ax.set_ylim(0, ymax0)\n",
    "        set_analysissubplot(key)\n",
    "        set_analysissubplot_special(key, poi_source)\n",
    "        ax.set_title(list(keys_metrics.values())[i+nc])\n",
    "        ax.set_xlabel('Growth strategy quantile')\n",
    "        if key in [\"poi_coverage\"]:\n",
    "            # https://stackoverflow.com/questions/30914462/matplotlib-how-to-force-integer-tick-labels\n",
    "            ax.set_yticks([0,40,80,120,160])\n",
    "        \n",
    "\n",
    "if plotnotitle:\n",
    "    plt.subplots_adjust(top = 0.92, bottom = 0.14, left = 0.05, right = 0.97, wspace = 0.28, hspace = 0.28)\n",
    "else:\n",
    "    plt.subplots_adjust(top = 0.82, bottom = 0.14, left = 0.05, right = 0.97, wspace = 0.28, hspace = 0.28)\n",
    "if plotconstricted:\n",
    "    fig.savefig(PATH[\"plots\"] + \"/\" + 'averagecity_analysis_poi_' + poi_source + '.eps', facecolor = \"white\", edgecolor = 'none')\n",
    "else:\n",
    "    fig.savefig(PATH[\"plots\"] + \"/\" + 'averagecity_analysis_poi_' + poi_source + '_noconstr.eps', facecolor = \"white\", edgecolor = 'none')\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Directness comparisons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%capture\n",
    "plotparam_analysis[\"constricted\"] = {\"linewidth\": 1.5, \"color\": '#D22A0E', \"linestyle\": \"solid\", \"label\": \"Betw. cars\"}\n",
    "plotparam_analysis[\"mst\"] = {\"linewidth\": 1, \"color\": '#aaaaaa', \"linestyle\": \"dotted\", \"label\": \"MST\"}\n",
    "plotparam_analysis[\"bikegrown_betweenness\"] = {\"linewidth\": 2.5, \"color\": '#0eb6d2', \"linestyle\": \"solid\", \"label\": \"Betw.\"}\n",
    "plt.style.use(PATH[\"parameters\"] + 'plotstyle.mplstyle')\n",
    "# Run all parameter sets\n",
    "prune_measure_list = [\"betweenness\", \"closeness\", \"random\"]\n",
    "poi_source = \"grid\" # railwaystation grid\n",
    "\n",
    "analysis_result_city = {}\n",
    "analysis_result_AVG = {}\n",
    "analysis_constricted_AVG = {}\n",
    "analysis_mst_result_city = {}\n",
    "for m in prune_measure_list:\n",
    "    analysis_result_city[m] = {}\n",
    "    numcities = 0\n",
    "    for placeid, placeinfo in tqdm(cities.items(), desc=\"Cities\"):\n",
    "\n",
    "        # PLOT Analysis\n",
    "        filename = placeid + '_poi_' + poi_source + \"_mst.csv\"\n",
    "        analysis_mst_result_city_temp = np.genfromtxt(PATH[\"results\"] + placeid + \"/\" + filename, delimiter=',', names=True)\n",
    "        \n",
    "        filename = placeid + '_poi_' + poi_source + \"_\" + m + \".csv\"\n",
    "        analysis_result_city_temp = np.genfromtxt(PATH[\"results\"] + placeid + \"/\" + filename, delimiter=',', names=True)\n",
    "        if len(analysis_result_city_temp) == 0: # Discard if no results (for example no railwaystations)\n",
    "            print(placeid + \": No analysis results available\")\n",
    "            continue\n",
    "        else:\n",
    "            numcities += 1\n",
    "            analysis_result_city[m][placeid] = analysis_result_city_temp\n",
    "            metric_keys = analysis_result_city[m][placeid].dtype.names\n",
    "            analysis_mst_result_city[placeid] = analysis_mst_result_city_temp\n",
    "\n",
    "    # All cities are loaded, now create the average values\n",
    "    analysis_result_AVG[m] = {}\n",
    "    for metric in metric_keys:\n",
    "        temp = np.zeros([40, numcities])\n",
    "        for i, placeid_this in enumerate(analysis_result_city[m].keys()):\n",
    "            temp[:, i] = analysis_result_city[m][placeid_this][metric]\n",
    "        analysis_result_AVG[m][metric] = np.mean(temp, axis = 1)\n",
    "        \n",
    "    analysis_mst_result_AVG = {}\n",
    "    for metric in metric_keys:\n",
    "        temp = np.zeros([1, numcities])\n",
    "        for i, placeid_this in enumerate(analysis_mst_result_city.keys()):\n",
    "            temp[:, i] = analysis_mst_result_city[placeid_this][metric]\n",
    "        analysis_mst_result_AVG[metric] = np.mean(temp, axis = 1)\n",
    "    print(analysis_mst_result_AVG)\n",
    "\n",
    "    analysis_constricted_AVG[m] = {}\n",
    "    for metric in list(range(9)): # 9 constricted fields are loaded\n",
    "        temp = np.zeros([40, len(analysis_constricted[poi_source][m].keys())]) \n",
    "        for i, placeid_this in enumerate(analysis_constricted[poi_source][m].keys()):\n",
    "            temp[:, i] = analysis_constricted[poi_source][m][placeid_this][:, metric]\n",
    "        analysis_constricted_AVG[m][metric] = np.mean(temp, axis = 1)\n",
    "\n",
    "# Plot\n",
    "nc = 4\n",
    "fig, axes = plt.subplots(nrows = 1, ncols = nc, figsize=(640/plotparam[\"dpi\"], 180/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"])\n",
    "keys_metrics = {\"directness_all_linkwise\": \"All CC, avg. of ratios\", \"directness\": \"All CC, ratio of totals\", \"directness_lcc_linkwise\": \"LCC, avg. of ratios\", \"directness_lcc\": \"LCC, ratio of totals\", }\n",
    "\n",
    "for i, ax in enumerate(axes):\n",
    "    for m in prune_measure_list:\n",
    "        key = list(keys_metrics.keys())[i]\n",
    "        ax.plot(prune_quantiles, analysis_result_AVG[m][key], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "    tmp, = ax.plot([0,1], [analysis_mst_result_AVG[key], analysis_mst_result_AVG[key]], **plotparam_analysis[\"mst\"])\n",
    "    ax.set_ylim(bottom = 0.43)\n",
    "    ax.set_ylim(top = 0.84)\n",
    "    ax.set_title(list(keys_metrics.values())[i])\n",
    "    ax.set_xlabel('Growth strategy quantile')\n",
    "    ax.set_xticks([0, 0.2, 0.4, 0.6, 0.8, 1])\n",
    "    ax.set_xlim(0, 1)\n",
    "    if i == 0:\n",
    "        ax.set_ylabel('Directness')\n",
    "#     if i == 3:\n",
    "#         leg = ax.legend(loc = 'best')\n",
    "#         leg.get_frame().set_linewidth(0.5)\n",
    "#         ax.set_zorder(100)\n",
    "\n",
    "if plotnotitle:\n",
    "    plt.subplots_adjust(top = 0.86, bottom = 0.2, left = 0.063, right = 0.98, wspace = 0.28, hspace = 0.28)\n",
    "else:\n",
    "    plt.subplots_adjust(top = 0.86, bottom = 0.2, left = 0.063, right = 0.98, wspace = 0.28, hspace = 0.28)\n",
    "if plotconstricted:\n",
    "    fig.savefig(PATH[\"plots\"] + \"/\" + 'averagecity_directnesscomparisons_poi_' + poi_source + '.eps', facecolor = \"white\", edgecolor = 'none')\n",
    "else:\n",
    "    fig.savefig(PATH[\"plots\"] + \"/\" + 'averagecity_directnesscomparisons_poi_' + poi_source + '_noconstr.eps', facecolor = \"white\", edgecolor = 'none')\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overlaps with existing bike infra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.ticker import PercentFormatter\n",
    "\n",
    "# %%capture\n",
    "plt.style.use(PATH[\"parameters\"] + 'plotstyle.mplstyle')\n",
    "# Run all parameter sets\n",
    "prune_measure_list = [\"betweenness\", \"closeness\", \"random\"]\n",
    "poi_source = \"railwaystation\" # railwaystation grid\n",
    "\n",
    "analysis_result_city = {}\n",
    "analysis_result_AVG = {}\n",
    "for m in prune_measure_list:\n",
    "    analysis_result_city[m] = {}\n",
    "    numcities = 0\n",
    "    for placeid, placeinfo in tqdm(cities.items(), desc=\"Cities\"):\n",
    "\n",
    "        # PLOT Analysis\n",
    "        filename = placeid + '_poi_' + poi_source + \"_\" + m + \".csv\"\n",
    "        analysis_result_city_temp = np.genfromtxt(PATH[\"results\"] + placeid + \"/\" + filename, delimiter=',', names=True)\n",
    "        if len(analysis_result_city_temp) == 0: # Discard if no results (for example no railwaystations)\n",
    "            print(placeid + \": No analysis results available\")\n",
    "            continue\n",
    "        else:\n",
    "            numcities += 1\n",
    "            analysis_result_city[m][placeid] = analysis_result_city_temp\n",
    "            metric_keys = analysis_result_city[m][placeid].dtype.names\n",
    "\n",
    "    # All cities are loaded, now create the average values\n",
    "    analysis_result_AVG[m] = {}\n",
    "    for metric in metric_keys:\n",
    "        temp = np.zeros([40, numcities])\n",
    "        for i, placeid_this in enumerate(analysis_result_city[m].keys()):\n",
    "            temp[:, i] = analysis_result_city[m][placeid_this][metric]\n",
    "        analysis_result_AVG[m][metric] = np.mean(temp, axis = 1)\n",
    "\n",
    "# Plot\n",
    "nc = 1\n",
    "fig, axes = plt.subplots(nrows = 2, ncols = nc, figsize=(140/plotparam[\"dpi\"], 320/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"])\n",
    "keys_metrics = {\"overlap_biketrack\": \"Overlap Protected\", \"overlap_bikeable\": \"Overlap Bikeable\"}\n",
    "\n",
    "i = 0\n",
    "ax = axes[i]\n",
    "key = list(keys_metrics.keys())[i]\n",
    "for m in prune_measure_list:\n",
    "    ax.plot(prune_quantiles, analysis_result_AVG[m][key] / analysis_result_AVG[m][\"length\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "    if not plotnotitle:\n",
    "        ax.text(-0.15, ymax0*1.25, \"Average city\" + \" (\" + poi_source + \")\", fontsize=16, horizontalalignment='left')\n",
    "\n",
    "set_analysissubplot(key)\n",
    "# set_analysissubplot_special(key, poi_source)\n",
    "ax.set_title(list(keys_metrics.values())[i])\n",
    "ax.set_xlabel('')\n",
    "ax.set_xticklabels([])\n",
    "ax.yaxis.set_major_formatter(matplotlib.ticker.PercentFormatter(1, 0))\n",
    "\n",
    "# Inset Copenhagen\n",
    "ax = fig.add_axes([0.54, 0.71, 0.39, 0.20]) # inset axes\n",
    "cityexample = \"copenhagen\"\n",
    "for m in prune_measure_list:\n",
    "    ax.plot(prune_quantiles, analysis_result[p][m][cityexample][key] / analysis_result[p][m][cityexample][\"length\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "    ax.set_ylim(0.42,0.88)\n",
    "    ax.set_yticks([0.5,0.6,0.7,0.8])\n",
    "    ax.yaxis.set_major_formatter(matplotlib.ticker.PercentFormatter(1, 0))\n",
    "    ax.text(0.96, 0.8, \"Copenhagen\", horizontalalignment = \"right\", fontsize = 6)\n",
    "set_analysissubplot(key)\n",
    "ax.set_xticklabels([])\n",
    "\n",
    "\n",
    "\n",
    "i = 1\n",
    "ax = axes[i]\n",
    "key = list(keys_metrics.keys())[i]\n",
    "for m in prune_measure_list:\n",
    "    ax.plot(prune_quantiles, analysis_result_AVG[m][key] / analysis_result_AVG[m][\"length\"], **plotparam_analysis[\"bikegrown_\" + m])\n",
    "#     leg = ax.legend(loc = 'best')\n",
    "#     leg.get_frame().set_linewidth(0.5)\n",
    "    \n",
    "set_analysissubplot(key)\n",
    "# set_analysissubplot_special(key, poi_source)\n",
    "ax.set_title(list(keys_metrics.values())[i])\n",
    "ax.yaxis.set_major_formatter(matplotlib.ticker.PercentFormatter(1, 0))\n",
    "ax.set_xlabel('Growth strategy quantile')\n",
    "\n",
    "if plotnotitle:\n",
    "    plt.subplots_adjust(top = 0.92, bottom = 0.11, left = 0.21, right = 0.95, wspace = 0.28, hspace = 0.28)\n",
    "else:\n",
    "    plt.subplots_adjust(top = 0.82, bottom = 0.11, left = 0.21, right = 0.95, wspace = 0.28, hspace = 0.28)\n",
    "\n",
    "fig.savefig(PATH[\"plots\"] + \"/\" + 'averagecity_analysisoverlap_poi_' + poi_source + '.eps', facecolor = \"white\", edgecolor = 'none')\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constricted analysis plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_constrictedsubplot_special(key):\n",
    "    if key == 0:\n",
    "        ax.set_ylim(0.6,0.8)\n",
    "    elif key == 1:\n",
    "        ax.set_ylim(0.095,0.115)\n",
    "    elif key == 2:\n",
    "        ax.set_ylim(0.6,0.82)\n",
    "    elif key == 3:\n",
    "        ax.set_ylim(0.58,0.8)\n",
    "    elif key == 4:\n",
    "        ax.set_ylim(0.48,0.57)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%capture\n",
    "\n",
    "fig, axess = plt.subplots(nrows = 6, ncols = 5, figsize=(640/plotparam[\"dpi\"], 6*140/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"])\n",
    "\n",
    "\n",
    "analysis_constricted_AVG = {}\n",
    "axescount = 0\n",
    "for p in [\"grid\",\"railwaystation\"]:\n",
    "    analysis_constricted_AVG[p] = {}\n",
    "    for m in prune_measure_list:\n",
    "        analysis_constricted_AVG[p][m] = {}\n",
    "        for metric in list(range(9)): # 9 constricted fields are loaded\n",
    "            temp = np.zeros([40, len(analysis_constricted[p][m].keys())]) \n",
    "            for i, placeid_this in enumerate(analysis_constricted[p][m].keys()):\n",
    "                temp[:, i] = analysis_constricted[p][m][placeid_this][:, metric]\n",
    "            analysis_constricted_AVG[p][m][metric] = np.mean(temp, axis = 1)\n",
    "\n",
    "        axes = axess[axescount]\n",
    "        axescount += 1\n",
    "        for i, ax in enumerate(axes):\n",
    "            if axescount == 1:\n",
    "                if i != 2:\n",
    "                    ax.set_title(constricted_plotinfo[\"title\"][i])\n",
    "                else:\n",
    "                    ax.set_title(\"Directness\")\n",
    "            if axescount == 6:\n",
    "                ax.set_xlabel('Growth strategy quantile')\n",
    "            ax.set_xticks([0, 0.2, 0.4, 0.6, 0.8, 1])\n",
    "            if axescount != 6:\n",
    "                ax.set_xticklabels([])\n",
    "            if i == 0:\n",
    "                ax.set_ylabel(p + \"\\n\" + m)\n",
    "            ax.set_xlim(-0.05, 1.05)\n",
    "            ax.spines['right'].set_visible(False)\n",
    "            ax.spines['top'].set_visible(False)\n",
    "            set_constrictedsubplot_special(i)\n",
    "\n",
    "            # Efficiency global\n",
    "            if i == 0:\n",
    "                ax.plot(prune_quantiles, analysis_constricted_AVG[p][m][0], **plotparam_analysis[\"constricted_SI\"])\n",
    "                ymin, ymax = ax.get_ylim()\n",
    "\n",
    "            # Efficiency local\n",
    "            elif i == 1:\n",
    "                ax.plot(prune_quantiles, analysis_constricted_AVG[p][m][1], **plotparam_analysis[\"constricted_SI\"])\n",
    "\n",
    "            # Directness\n",
    "            elif i == 2:\n",
    "                ax.plot(prune_quantiles, analysis_constricted_AVG[p][m][8], **plotparam_analysis[\"constricted_SI\"])\n",
    "\n",
    "            # Clustering\n",
    "            elif i == 3:\n",
    "                ax.plot(prune_quantiles, analysis_constricted_AVG[p][m][2], **plotparam_analysis[\"constricted_10\"])\n",
    "                ax.plot(prune_quantiles, analysis_constricted_AVG[p][m][3], **plotparam_analysis[\"constricted_5\"])\n",
    "                ax.plot(prune_quantiles, analysis_constricted_AVG[p][m][4], **plotparam_analysis[\"constricted_3\"])\n",
    "#                 ax.legend(loc = 'best')\n",
    "                if axescount == 1:\n",
    "                    ax.text(0.45, 0.735, plotparam_analysis[\"constricted_10\"][\"label\"], fontsize = 7, color = plotparam_analysis[\"constricted_10\"][\"color\"])\n",
    "                    ax.text(0.45, 0.655, plotparam_analysis[\"constricted_5\"][\"label\"], fontsize = 7, color = plotparam_analysis[\"constricted_5\"][\"color\"])\n",
    "                    ax.text(0.45, 0.6, plotparam_analysis[\"constricted_3\"][\"label\"], fontsize = 7, color = plotparam_analysis[\"constricted_3\"][\"color\"])\n",
    "\n",
    "            # Anisotropy\n",
    "            elif i == 4:\n",
    "                ax.plot(prune_quantiles, analysis_constricted_AVG[p][m][5], **plotparam_analysis[\"constricted_10\"])\n",
    "                ax.plot(prune_quantiles, analysis_constricted_AVG[p][m][6], **plotparam_analysis[\"constricted_5\"])\n",
    "                ax.plot(prune_quantiles, analysis_constricted_AVG[p][m][7], **plotparam_analysis[\"constricted_3\"])\n",
    "\n",
    "plt.subplots_adjust(top = 0.95, bottom = 0.05, left = 0.09, right = 0.97, wspace = 0.35)\n",
    "fig.savefig(PATH[\"plots\"] + \"/\" + 'averagecity_constricted.eps', facecolor = \"white\", edgecolor = 'none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic stats on all cities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### World map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "world = world[(world.name!=\"Antarctica\")] #### exclude Antartica\n",
    "\n",
    "\n",
    "cities = pd.read_csv(PATH['parameters'] + 'cities_all_lat_lng.csv') #### read cities lat and lon\n",
    "\n",
    "cities_x = np.array(cities['lng'].tolist())\n",
    "cities_y = np.array(cities['lat'].tolist())\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(640/plotparam[\"dpi\"], 320/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"], squeeze = True)\n",
    "ax.set_facecolor('aliceblue')  #### color background with blue to represent oceans\n",
    "\n",
    "\n",
    "world.plot(ax=ax,color='lightgrey', edgecolor='whitesmoke', linewidth=0.5)\n",
    "plt.plot(cities_x,cities_y,'.',color='black', markersize=6) # 'o',markersize=5, markerfacecolor='black',markeredgewidth=0.8, markeredgecolor='grey'\n",
    "ax.set_ylim([-65, 90])\n",
    "ax.set_xlim([-180, 180])\n",
    "plt.xticks([])\n",
    "plt.yticks([])\n",
    "\n",
    "\n",
    "\n",
    "###################  Europe inset\n",
    "axins = ax.inset_axes([0.01, 0.01, 0.251, 0.48]) ### location and size of the inset [x_0,y_0, width, height]\n",
    "\n",
    "world.plot(color='lightgrey', edgecolor='whitesmoke',ax=axins, linewidth=0.6)\n",
    "axins.plot(cities_x,cities_y,'.',color='black', markersize=7)\n",
    "axins.set_facecolor('aliceblue')\n",
    "# sub region of the original image\n",
    "x1, x2, y1, y2 = -5, 22, 40.5, 61 #### coordinates to zoom in\n",
    "axins.set_xlim(x1, x2)\n",
    "axins.set_ylim(y1, y2)\n",
    "axins.set_xticks([])\n",
    "axins.set_yticks([])\n",
    "# ax.indicate_inset_zoom(axins, edgecolor=\"black\")  ### lines towards the inset\n",
    "plt.tight_layout()\n",
    "# plt.subplots_adjust(top = 1.12, bottom = 0, left = 0, right = 1)\n",
    "fig.savefig('world_insets.eps',bbox='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "world = world[(world.name!=\"Antarctica\")] #### exclude Antartica\n",
    "\n",
    "\n",
    "cities = pd.read_csv(PATH['parameters'] + 'cities_all_lat_lng.csv') #### read cities lat and lon\n",
    "\n",
    "cities_x = np.array(cities['lng'].tolist())\n",
    "cities_y = np.array(cities['lat'].tolist())\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(310/plotparam[\"dpi\"], 240/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"])\n",
    "\n",
    "ax = fig.add_subplot(3,2,(1,4))\n",
    "ax.set_facecolor('aliceblue')  #### color background with blue to represent oceans\n",
    "\n",
    "\n",
    "world.plot(ax=ax,color='lightgrey', edgecolor='whitesmoke', linewidth=0.4)\n",
    "plt.plot(cities_x,cities_y,'.',color='black', markersize=1.5) # 'o',markersize=5, markerfacecolor='black',markeredgewidth=0.8, markeredgecolor='grey'\n",
    "ax.set_ylim([-65, 90])\n",
    "ax.set_xlim([-180, 180])\n",
    "plt.xticks([])\n",
    "plt.yticks([])\n",
    "\n",
    "\n",
    "###################  Europe inset\n",
    "axins = ax.inset_axes([0.01, 0.01, 0.251, 0.48]) ### location and size of the inset [x_0,y_0, width, height]\n",
    "\n",
    "world.plot(color='lightgrey', edgecolor='whitesmoke',ax=axins, linewidth=0.5)\n",
    "axins.plot(cities_x,cities_y,'.',color='black', markersize=2)\n",
    "axins.set_facecolor('aliceblue')\n",
    "# sub region of the original image\n",
    "x1, x2, y1, y2 = -5, 22, 40.5, 61 #### coordinates to zoom in\n",
    "axins.set_xlim(x1, x2)\n",
    "axins.set_ylim(y1, y2)\n",
    "axins.set_xticks([])\n",
    "axins.set_yticks([])\n",
    "# ax.indicate_inset_zoom(axins, edgecolor=\"black\")  ### lines towards the inset\n",
    "plt.tight_layout()\n",
    "# plt.subplots_adjust(top = 1.12, bottom = 0, left = 0, right = 1)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Left subplot: Lengths\n",
    "axes = fig.add_subplot(325)\n",
    "axes.hist([i/1000 for i in datastats[\"length_biketrack\"]], bins = np.linspace(0,1400,15), edgecolor=\"w\")\n",
    "axes.set_ylim([0, 19])\n",
    "axes.set_yticks([0,5,10,15]);\n",
    "axes.set_xlabel(\"Bicycle network length [km]\")\n",
    "axes.set_ylabel(\"Cities\")\n",
    "axes.set_xticks([0,400,800,1200]);\n",
    "labels = axes.set_xticklabels([0,400,800,1200])\n",
    "# https://stackoverflow.com/questions/14852821/aligning-rotated-xticklabels-with-their-respective-xticks\n",
    "for i, label in enumerate(labels):\n",
    "    labelypos = label.get_position()[1]\n",
    "\n",
    "# Right subplot: Components\n",
    "axes = fig.add_subplot(326)\n",
    "axes.hist([i for i in datastats[\"components\"] if i != 0], bins=np.logspace(0, 4, 17), edgecolor=\"w\")\n",
    "\n",
    "axes.set_ylim([0, 19])\n",
    "axes.set_yticks([0,5,10,15]);\n",
    "axes.set_yticklabels([])\n",
    "axes.set_xlabel(\"Disconnected components\");\n",
    "# labels = axes.set_xticklabels([1,10,100,1000,10000])\n",
    "axes.set_xscale('log') \n",
    "axes.set_xticks([1,10,100,1000,10000]);\n",
    "# for i, label in enumerate(labels):\n",
    "#     label.set_y(labelypos)\n",
    "\n",
    "# # Right subplot: Street Lengths\n",
    "# axes = fig.add_subplot(326)\n",
    "# axes.hist([i/1000 for i in datastats[\"length_carall\"]], bins = np.linspace(0,14000,15), edgecolor=\"w\", color=\"red\")\n",
    "# axes.set_ylim([0, 19])\n",
    "# axes.set_yticks([0,5,10,15]);\n",
    "# axes.set_yticklabels([])\n",
    "# axes.set_xlabel(\"Street network length [km]\")\n",
    "# axes.set_xticks([0,4000,8000,12000]);\n",
    "# labels = axes.set_xticklabels([0,4000,8000,12000])\n",
    "\n",
    "    \n",
    "\n",
    "plt.subplots_adjust(top = 1, bottom = 0.16, left = 0.12, right = 0.99, wspace = 0.1, hspace = 0.1)\n",
    "fig.savefig(PATH[\"plots\"] + 'world_insets.eps',bbox='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stats / Histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datastats = {}\n",
    "datastats[\"length_biketrack\"] = []\n",
    "datastats[\"length_biketrack_onstreet\"] = []\n",
    "datastats[\"length_carall\"] = []\n",
    "datastats[\"components\"] = []\n",
    "for placeid, placeinfo in cities.items():\n",
    "    datastats[\"length_biketrack\"].append(analysis_existing[placeid][\"length\"][analysis_existing_rowkeys[\"biketrack\"]])\n",
    "    datastats[\"length_biketrack_onstreet\"].append(analysis_existing[placeid][\"length\"][analysis_existing_rowkeys[\"biketrack_onstreet\"]])\n",
    "    datastats[\"length_carall\"].append(analysis_existing[placeid][\"length\"][analysis_existing_rowkeys[\"carall\"]])\n",
    "    datastats[\"components\"].append(analysis_existing[placeid][\"components\"][analysis_existing_rowkeys[\"biketrack\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datastats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(320/plotparam[\"dpi\"], 120/plotparam[\"dpi\"]), dpi=plotparam[\"dpi\"], squeeze = True)\n",
    "\n",
    "# Left subplot: Lengths\n",
    "axes = ax[0]\n",
    "axes.hist([i/1000 for i in datastats[\"length_biketrack\"]], bins = np.linspace(0,1400,15), edgecolor=\"w\")\n",
    "axes.set_ylim([0, 19])\n",
    "axes.set_xlabel(\"Bicycle network length [km]\")\n",
    "axes.set_ylabel(\"Number of cities\")\n",
    "axes.set_xticks([0,400,800,1200]);\n",
    "\n",
    "# Right subplot: Components\n",
    "axes = ax[1]\n",
    "axes.hist([i for i in datastats[\"components\"] if i != 0], bins=np.logspace(0, 4, 17), edgecolor=\"w\")\n",
    "axes.set_xscale('log') \n",
    "axes.set_ylim([0, 19])\n",
    "axes.set_yticklabels([])\n",
    "axes.set_xlabel(\"Disconnected components\");\n",
    "axes.set_xticks([1,10,100,1000,10000]);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(datastats[\"length_biketrack\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.hist(np.log10([i for i in datastats[\"components\"] if i != 0]), bins = 20)\n",
    "\n",
    "fig = plt.figure(figsize=(4, 3))\n",
    "axes = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "axes.hist([i for i in datastats[\"components\"] if i != 0], bins=np.logspace(0, 4, 20), edgecolor=\"w\")\n",
    "\n",
    "axes.set_xscale('log') \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist([i/j for (i,j) in zip(datastats[\"length_biketrack\"], datastats[\"length_carall\"])], bins = 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[(i/j, p) for (i,j,p) in zip(datastats[\"length_biketrack\"], datastats[\"length_carall\"], list(cities.keys()))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist([i/j for (i,j) in zip(datastats[\"length_biketrack_onstreet\"], datastats[\"length_carall\"])])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "OSMNX",
   "language": "python",
   "name": "osmnx"
  },
  "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}