Data_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bike Sharing and Urban Mobility in a Post-Pandemic World #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import pickle\n",
    "import pygsp\n",
    "import mplleaflet\n",
    "import scipy as sp\n",
    "\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "from utils.indexes import days_index, hour_rides, transform_source\n",
    "from utils.utils import compute_flows, compute_flows_between, weighted_coverage, find_partitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from networkx.algorithms.community.modularity_max import greedy_modularity_communities\n",
    "from networkx.algorithms.community.quality import performance, coverage\n",
    "from networkx.algorithms.community.centrality import girvan_newman\n",
    "from networkx.algorithms.community.asyn_fluid import asyn_fluidc\n",
    "\n",
    "import community \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Timeline</h3>\n",
    "\n",
    "March 1st: first reported case in New York \n",
    "<br>\n",
    "March 4th: Yeshiva University campus closed, high school in the Bronx \n",
    "<br>\n",
    "March 9th 16 confirmed cases \n",
    "<br>\n",
    "March 10th: mitigation measures expanded (online classes for universities) \n",
    "<br>\n",
    "March 11th: CUNY and SUNY closed\n",
    "<br>\n",
    "March 14th: first reported death in New York \n",
    "<br>\n",
    "March 16th: NYC public schools closed \n",
    "<br>\n",
    "March 20th: stay-at-home order for non-essential workers, total over 7000 cases\n",
    "<br>\n",
    "March 22nd: shortages of PPE for health workers, recommendation that health facilities stop testing non-hospitalized patients\n",
    "<br>\n",
    "March 26th: USNS Comfort heading to NYC to assist local hospitals \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "[New York Weather Data](https://www.timeanddate.com/weather/usa/new-york/historic?month=3&year=2020)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Collect and organize Data #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collect Bike Trips Data\n",
    "\n",
    "df = pd.read_csv('./data/202003-citibike-tripdata.csv')\n",
    "n_start_station = len(df['start station id'].unique())\n",
    "n_end_station = len(df['end station id'].unique())\n",
    "\n",
    "\n",
    "total_station_id = set(df['start station id']).union(set(df['end station id']))\n",
    "n_tot_station = len(total_station_id)\n",
    "\n",
    "id_index = dict(zip(sorted(total_station_id), np.arange(n_tot_station)))\n",
    "\n",
    "# Create Dictionary Station : Position\n",
    "\n",
    "locations = dict()\n",
    "for e in id_index.keys():\n",
    "    if df[df['start station id'] == e]['start station latitude'].shape[0]:\n",
    "        locations[id_index[e]] = (df[df['start station id'] == e]['start station longitude'].iloc[0],\n",
    "                                  df[df['start station id'] == e]['start station latitude'].iloc[0])\n",
    "    else:\n",
    "        locations[id_index[e]] = (df[df['end station id'] == e]['end station longitude'].iloc[0],\n",
    "                                  df[df['end station id'] == e]['end station latitude'].iloc[0])\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Dictionary with March daily trips and store in file\n",
    "# !!! USE IT JUST ONCE AND THEN LOAD PICKLE FILE !!!\n",
    "\n",
    "\"\"\"d = days_index(df)\n",
    "with open('./data/variables/march_index.pickle', 'wb') as file:\n",
    "    pickle.dump(d, file)\"\"\"\n",
    "\n",
    "\n",
    "# Load March daily trips\n",
    "with open('./data/python_variables/march_index.pickle', 'rb') as file:\n",
    "    d = pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monthly View \n",
    "\n",
    "d_a = [e for e in d.values()]\n",
    "d_a.append(df.shape[0])\n",
    "d_a = np.array(d_a)\n",
    "rides_day = d_a[1:] - d_a[:-1]\n",
    "\n",
    "# Rides Per Day\n",
    "plt.figure(figsize=(12,8))\n",
    "_ = plt.plot(np.arange(1, 32), rides_day, color='blue', marker='*')\n",
    "# _ = plt.xticks(np.arange(1, 8), [1, 5, 10, 15, 20, 25, 31], fontsize=18)\n",
    "_ = plt.xlabel(\"Day of March\", fontsize=24, color='black')\n",
    "_ = plt.ylabel(\"Rides per day\", fontsize=24, color='black')\n",
    "_ = plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)\n",
    "plt.grid(False)\n",
    "plt.savefig('./plots/days_aggregated.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Maximum number of rides on day {}\".format(np.argmax(rides_day)))\n",
    "print(\"With {} rides\".format(np.max(rides_day)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Threshold for edge creation \n",
    "ths = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Discard stations below a specific Longitude (3 stations)\n",
    "out_nodes = []\n",
    "for k, v in locations.items():\n",
    "    if v[0] < -74.02:\n",
    "        out_nodes.append(k)\n",
    "out_nodes = set(out_nodes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> First Week of March (2-6) </h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "days = [2, 3, 4, 5, 6]\n",
    "\n",
    "# Find adjacency matrix\n",
    "adj_2_6_tot = np.zeros((n_tot_station, n_tot_station))\n",
    "for day in days:\n",
    "    for index, row in df.iloc[d[day]:d[day+1]].iterrows():      \n",
    "        adj_2_6_tot[id_index[row['start station id']], id_index[row['end station id']]] += 1\n",
    "        adj_2_6_tot[id_index[row['end station id']], id_index[row['start station id']]] +=1\n",
    "    print('Day {} loaded...'.format(day))\n",
    "    \n",
    "\n",
    "print(\"Total number of rides : {}\".format(np.sum(adj_2_6_tot) / 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove edges with weight <= threshold and self loops\n",
    "\n",
    "adj_2_6 = adj_2_6_tot.copy()\n",
    "adj_2_6[adj_2_6 <= ths] = 0\n",
    "for i in range(n_tot_station):\n",
    "    adj_2_6[i, i] = 0\n",
    "print(\"Unweighted number of rides (edges) : {}\".format(np.sum(adj_2_6) / 2))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Undirected graph and plot degree distribution\n",
    "\n",
    "plt.figure(figsize=(12,8))\n",
    "g_2_6 = nx.from_numpy_matrix(adj_2_6)\n",
    "g_2_6_tot = nx.from_numpy_matrix(adj_2_6)\n",
    "plt.xlabel('Degree')\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.grid(False)\n",
    "degs_2_6 = np.array([x[1] for x in g_2_6.degree])\n",
    "_ = plt.hist(degs_2_6, color='blue')\n",
    "plt.savefig('./plots/degree_2_6.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Entropy of Distribution\n",
    "print(\"Entropy {}\".format(sp.stats.entropy(np.histogram(degs_2_6, bins=12)[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Clustering </h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove isolated compontents (Discard very small percentsge of nodes) and analyse spectrum\n",
    "\n",
    "isolated_nodes_2_6 = []\n",
    "for e in nx.connected_components(g_2_6):\n",
    "    if len(e) < 3:\n",
    "        isolated_nodes_2_6.append(e.pop())\n",
    "\n",
    "isolated_nodes_2_6.extend(out_nodes)\n",
    "isolated_nodes_2_6 = list(set(isolated_nodes_2_6))\n",
    "g_2_6.remove_nodes_from(isolated_nodes_2_6)\n",
    "\n",
    "for e in nx.connected_components(g_2_6):\n",
    "    if len(e) < 3:\n",
    "        print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G_2_6 = pygsp.graphs.Graph(nx.adjacency_matrix(g_2_6))\n",
    "\n",
    "G_2_6.compute_laplacian()\n",
    "G_2_6.compute_fourier_basis()\n",
    "\n",
    "print(\"ACC: {}\".format(nx.algorithms.cluster.average_clustering(g_2_6)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find Best partition over 20 Iterations of Louvan Algorithm\n",
    "\n",
    "max_mod_2_6 = 0\n",
    "part_2_6 = None\n",
    "\n",
    "for i in range(20):\n",
    "\n",
    "    partition = community.best_partition(g_2_6, weight='weight')\n",
    "    c_g = [set() for e in range(len(partition))]\n",
    "    for k, v in partition.items():\n",
    "        c_g[v].add(k)\n",
    "    mod = community.modularity(partition, g_2_6, weight='weight')\n",
    "    if max_mod_2_6 < mod:\n",
    "        max_mod_2_6 = mod\n",
    "        part_2_6 = partition\n",
    "\n",
    "\n",
    "print(\"weighted coverage : {}\".format(weighted_coverage(g_2_6, part_2_6)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show Clusters\n",
    "\n",
    "color_cluster = ['yellow', 'blue', 'red', 'gray', 'green']\n",
    "vec_col_2_6 = [color_cluster[e] for e in part_2_6.values()]\n",
    "\n",
    "plt.figure(figsize=(8,12))\n",
    "\n",
    "nx.draw_networkx_nodes(g_2_6, locations, node_size=60, node_color=vec_col_2_6)\n",
    "nx.draw_networkx_edges(g_2_6, locations, alpha=0.8, width=0.1, edge_color='black')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Second Week of March (9-13) </h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ov = 0\n",
    "ov2 = 0\n",
    "for k,v in locations.items():\n",
    "    if v[0] < locations[ov][0]:\n",
    "        ov2 = ov\n",
    "        ov = k     \n",
    "    elif v[0] < locations[ov2][0]:\n",
    "        ov2 = k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9 -> Monday\n",
    "\n",
    "days = [9, 10, 11, 12, 13]\n",
    "\n",
    "# Find adjacency matrix\n",
    "adj_9_13_tot = np.zeros((n_tot_station, n_tot_station))\n",
    "for day in days:\n",
    "    for index, row in df.iloc[d[day]:d[day+1]].iterrows():\n",
    "        adj_9_13_tot[id_index[row['start station id']], id_index[row['end station id']]] += 1\n",
    "        adj_9_13_tot[id_index[row['end station id']], id_index[row['start station id']]] += 1\n",
    "    print('Day {} loaded...'.format(day))\n",
    "    \n",
    "\n",
    "print(\"Total number of rides : {}\".format(np.sum(adj_9_13_tot) / 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove Self-loops and edges with weight <= threshold\n",
    "\n",
    "adj_9_13 = adj_9_13_tot.copy()\n",
    "adj_9_13[adj_9_13 <= ths] = 0\n",
    "\n",
    "for i in range(n_tot_station):\n",
    "    adj_9_13[i, i] = 0\n",
    "\n",
    "print(\"Unweighted number of rides : {}\".format(np.sum(adj_9_13) / 2))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Undirected graph  and plot degree distribution \n",
    "\n",
    "plt.figure(figsize=(12,8))\n",
    "g_9_13 = nx.from_numpy_matrix(adj_9_13)\n",
    "g_9_13_tot = nx.from_numpy_matrix(adj_9_13)\n",
    "plt.xlabel('Degree')\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.grid(False)\n",
    "degs_9_13 = np.array([x[1] for x in g_9_13.degree])\n",
    "_ = plt.hist(degs_9_13, bins=12, color='blue')\n",
    "plt.savefig('./plots/degree_9_13.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot Entropy of Degree distribution\n",
    "\n",
    "print(\"Entropy {}\".format(sp.stats.entropy(np.histogram(degs_9_13, bins=12)[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Clustering </h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove isolated components (small percentage of nodes)\n",
    "\n",
    "isolated_nodes_9_13 = []\n",
    "for e in nx.connected_components(g_9_13):\n",
    "    if len(e) < 3:\n",
    "        isolated_nodes_9_13.append(e.pop())\n",
    "\n",
    "isolated_nodes_9_13.extend(out_nodes)\n",
    "isolated_nodes_9_13 = list(set(isolated_nodes_9_13))\n",
    "g_9_13.remove_nodes_from(isolated_nodes_9_13)\n",
    "\n",
    "for e in nx.connected_components(g_9_13):\n",
    "    if len(e) < 3:\n",
    "        print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Study spectrum\n",
    "G_9_13 = pygsp.graphs.Graph(nx.adjacency_matrix(g_9_13))\n",
    "\n",
    "G_9_13.compute_laplacian()\n",
    "G_9_13.compute_fourier_basis()\n",
    "\n",
    "print(\"ACC: {}\".format(nx.algorithms.cluster.average_clustering(g_9_13)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find best partition\n",
    "\n",
    "max_mod_9_13 = 0\n",
    "part_9_13 = None\n",
    "\n",
    "for i in range(20):\n",
    "\n",
    "    partition = community.best_partition(g_9_13, weight='weight')\n",
    "    c_g = [set() for e in range(len(partition))]\n",
    "    for k, v in partition.items():\n",
    "        c_g[v].add(k)\n",
    "    mod = community.modularity(partition, g_9_13, weight='weight')\n",
    "    if max_mod_9_13 < mod:\n",
    "        max_mod_9_13 = mod\n",
    "        part_9_13 = partition\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"weighted coverage : {}\".format(weighted_coverage(g_9_13, part_9_13)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show Clusters\n",
    "\n",
    "color_cluster = ['yellow', 'blue', 'red', 'gray', 'green']\n",
    "vec_col_9_13 = [color_cluster[e] for e in part_9_13.values()]\n",
    "\n",
    "plt.figure(figsize=(8, 12))\n",
    "\n",
    "nx.draw_networkx_nodes(g_9_13, locations, node_size=40, node_color=vec_col_9_13)\n",
    "nx.draw_networkx_edges(g_9_13, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.axis('off')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Third Week of March (16-20) </h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 16 -> Monday\n",
    "\n",
    "days = [16, 17, 18, 19, 20]\n",
    "\n",
    "# Find adjacency matrix\n",
    "adj_16_20_tot = np.zeros((n_tot_station, n_tot_station))\n",
    "for day in days:\n",
    "    for index, row in df.iloc[d[day]:d[day+1]].iterrows():\n",
    "        adj_16_20_tot[id_index[row['start station id']], id_index[row['end station id']]] += 1\n",
    "        adj_16_20_tot[id_index[row['end station id']], id_index[row['start station id']]] += 1\n",
    "    print('Day {} loaded...'.format(day))\n",
    "    \n",
    "\n",
    "print(\"Total number of rides : {}\".format(np.sum(adj_16_20_tot) / 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove edges with weight <= threshold and Self-loops\n",
    "adj_16_20 = adj_16_20_tot.copy()\n",
    "adj_16_20[adj_16_20 <= ths] = 0\n",
    "\"\"\"adj_16_20[adj_16_20 >= ths] = 1\"\"\"\n",
    "\n",
    "print(\"Unweighted number of rides : {}\".format(np.sum(adj_16_20) / 2))\n",
    "\n",
    "# Remove Self-loops\n",
    "for i in range(n_tot_station):\n",
    "    adj_16_20[i, i] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Undirected graph and plot degree distribution\n",
    "plt.figure(figsize=(12,8))\n",
    "g_16_20 = nx.from_numpy_matrix(adj_16_20)\n",
    "g_16_20_tot = nx.from_numpy_matrix(adj_16_20)\n",
    "plt.xlabel('Degree')\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.grid(False)\n",
    "degs_16_20 = np.array([x[1] for x in g_16_20.degree])\n",
    "_ = plt.hist(degs_16_20, color='blue')\n",
    "plt.savefig('./plots/degree_16_20.eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Clustering </h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove isolated components (small percentage of nodes)\n",
    "isolated_nodes_16_20 = []\n",
    "for e in nx.connected_components(g_16_20):\n",
    "    if len(e) < 3:\n",
    "        isolated_nodes_16_20.append(e.pop())\n",
    "isolated_nodes_16_20.extend(out_nodes)\n",
    "isolated_nodes_16_20 = list(set(isolated_nodes_16_20))\n",
    "g_16_20.remove_nodes_from(isolated_nodes_16_20)\n",
    "\n",
    "for e in nx.connected_components(g_16_20):\n",
    "    if len(e) < 3:\n",
    "        print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolated_nodes_16_20_ = []\n",
    "for e in nx.connected_components(g_16_20):\n",
    "    if len(e) < 3:\n",
    "        isolated_nodes_16_20_.append(e.pop())\n",
    "isolated_nodes_16_20 = list(set(isolated_nodes_16_20))\n",
    "g_16_20.remove_nodes_from(isolated_nodes_16_20_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in nx.connected_components(g_16_20):\n",
    "    if len(e) < 3:\n",
    "        print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolated_nodes_16_20.extend(isolated_nodes_16_20_)\n",
    "isolated_nodes_16_20 = list(set(isolated_nodes_16_20))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G_16_20 = pygsp.graphs.Graph(nx.adjacency_matrix(g_16_20))\n",
    "\n",
    "G_16_20.compute_laplacian()\n",
    "G_16_20.compute_fourier_basis()\n",
    "\n",
    "\n",
    "print(\"ACC: {}\".format(nx.algorithms.cluster.average_clustering(g_16_20)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find best partition\n",
    "max_mod_16_20 = 0\n",
    "part_16_20 = None\n",
    "\n",
    "for i in range(20):\n",
    "\n",
    "    partition = community.best_partition(g_16_20, weight='weight')\n",
    "    c_g = [set() for e in range(len(partition))]\n",
    "    for k, v in partition.items():\n",
    "        c_g[v].add(k)\n",
    "    mod = community.modularity(partition, g_16_20, weight='weight')\n",
    "    if max_mod_16_20 < mod:\n",
    "        max_mod_16_20 = mod\n",
    "        part_16_20 = partition\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"weighted coverage : {}\".format(weighted_coverage(g_16_20, part_16_20)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show Clusters\n",
    "\n",
    "color_cluster = ['gray', 'yellow', 'red', 'blue', 'green']\n",
    "vec_col_16_20 = [color_cluster[e] for e in part_16_20.values()]\n",
    "\n",
    "plt.figure(figsize=(8, 12))\n",
    "\n",
    "nx.draw_networkx_nodes(g_16_20, locations, node_size=40, node_color=vec_col_16_20)\n",
    "nx.draw_networkx_edges(g_16_20, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Fourth Week of March (23-27) </h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 23 -> Monday\n",
    "\n",
    "days = [23, 24, 25, 26, 27]\n",
    "\n",
    "# Find adjacency matrix\n",
    "adj_23_27_tot = np.zeros((n_tot_station, n_tot_station))\n",
    "for day in days:\n",
    "    for index, row in df.iloc[d[day]:d[day+1]].iterrows():\n",
    "        adj_23_27_tot[id_index[row['start station id']], id_index[row['end station id']]] += 1\n",
    "        adj_23_27_tot[id_index[row['end station id']], id_index[row['start station id']]] += 1\n",
    "    print('Day {} loaded...'.format(day))\n",
    "    \n",
    "\n",
    "print(\"Total number of rides : {}\".format(np.sum(adj_23_27_tot) / 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove edges with weights <= threshold and Self-loops\n",
    "\n",
    "adj_23_27 = adj_23_27_tot.copy()\n",
    "adj_23_27[adj_23_27 <= ths] = 0\n",
    "\n",
    "for i in range(n_tot_station):\n",
    "    adj_23_27[i, i] = 0\n",
    "\n",
    "print(\"Unweighted number of rides : {}\".format(np.sum(adj_23_27) / 2))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Undirected graph and plot degree distribution\n",
    "plt.figure(figsize=(12,8))\n",
    "g_23_27 = nx.from_numpy_matrix(adj_23_27)\n",
    "g_23_27_tot = nx.from_numpy_matrix(adj_23_27)\n",
    "plt.xlabel('Degree')\n",
    "plt.ylabel(\"Frequency\")\n",
    "plt.grid(False)\n",
    "degs_23_27 = np.array([x[1] for x in g_23_27.degree])\n",
    "_ = plt.hist(degs_23_27, bins=12, color='blue')\n",
    "plt.savefig('./plots/degree_23_29.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Entropy of degree distribution\n",
    "print(\"Entropy {}\".format(sp.stats.entropy(np.histogram(degs_23_27, bins=12)[0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Clustering </h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolated_nodes_23_27 = []\n",
    "for e in nx.connected_components(g_23_27):\n",
    "    if len(e) < 3:\n",
    "        isolated_nodes_23_27.append(e.pop())\n",
    "\n",
    "isolated_nodes_23_27.extend(out_nodes)\n",
    "isolated_nodes_23_27 = list(set(isolated_nodes_23_27))\n",
    "g_23_27.remove_nodes_from(isolated_nodes_23_27)\n",
    "\n",
    "for e in nx.connected_components(g_23_27):\n",
    "    if len(e) < 3:\n",
    "        print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G_23_27 = pygsp.graphs.Graph(nx.adjacency_matrix(g_23_27))\n",
    "\n",
    "G_23_27.compute_laplacian()\n",
    "G_23_27.compute_fourier_basis()\n",
    "\n",
    "print(\"ACC: {}\".format(nx.algorithms.cluster.average_clustering(g_23_26)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find best partition\n",
    "\n",
    "max_mod_23_27 = 0\n",
    "part_23_27 = None\n",
    "\n",
    "for i in range(20):\n",
    "\n",
    "    partition = community.best_partition(g_23_27, weight='weight')\n",
    "    c_g = [set() for e in range(len(partition))]\n",
    "    for k, v in partition.items():\n",
    "        c_g[v].add(k)\n",
    "    mod = community.modularity(partition, g_23_27, weight='weight')\n",
    "    if max_mod_23_27 < mod:\n",
    "        max_mod_23_27 = mod\n",
    "        part_23_27 = partition\n",
    "\n",
    "\n",
    "print(\"weighted coverage : {}\".format(weighted_coverage(g_23_27, part_23_27)))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show Clusters\n",
    "color_cluster = ['yellow', 'blue', 'red', 'gray', 'green']\n",
    "vec_col_23_27 = [color_cluster[e] for e in part_23_27.values()]\n",
    "\n",
    "plt.figure(figsize=(8, 12))\n",
    "\n",
    "nx.draw_networkx_nodes(g_23_27, locations, node_size=60, node_color=vec_col_23_27)\n",
    "nx.draw_networkx_edges(g_23_27, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_man_idx = c_g[2]\n",
    "queens_idx = c_g[4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot all degree distributions\n",
    "\n",
    "plt.figure(figsize=(14,12))\n",
    "\n",
    "plt.subplot(2,2,1)\n",
    "plt.xlabel('Degree', fontsize=24)\n",
    "plt.ylabel(\"Frequency\", fontsize=24)\n",
    "plt.grid(False)\n",
    "plt.xticks(fontsize=18)\n",
    "plt.yticks(fontsize=18)\n",
    "plt.ylim(0, 250)\n",
    "_ = plt.hist(degs_2_6, bins=12, color='blue')\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "plt.xlabel('Degree', fontsize=24)\n",
    "plt.ylabel(\"Frequency\", fontsize=24)\n",
    "plt.grid(False)\n",
    "plt.xticks(fontsize=18)\n",
    "plt.yticks(fontsize=18)\n",
    "plt.ylim(0, 250)\n",
    "_ = plt.hist(degs_9_13, bins=12, color='blue')\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "plt.xlabel('Degree', fontsize=24)\n",
    "plt.ylabel(\"Frequency\", fontsize=24)\n",
    "plt.grid(False)\n",
    "plt.xticks(fontsize=18)\n",
    "plt.yticks(fontsize=18)\n",
    "plt.ylim(0, 250)\n",
    "_ = plt.hist(degs_16_20, bins=12, color='blue')\n",
    "\n",
    "plt.subplot(2,2,4)\n",
    "plt.xlabel('Degree', fontsize=24)\n",
    "plt.ylabel(\"Frequency\", fontsize=24)\n",
    "plt.grid(False)\n",
    "plt.xticks(fontsize=18)\n",
    "plt.yticks(fontsize=18)\n",
    "plt.ylim(0, 250)\n",
    "_ = plt.hist(degs_23_27, bins=12, color='blue')\n",
    "\n",
    "plt.savefig('plots/degs_hist.eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Degree Distributions sorted by First Week ascending order </h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sort_deg = np.argsort(degs_2_6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 20))\n",
    "\n",
    "plt.subplot(2,2,1)\n",
    "plt.scatter(np.arange(len(degs_2_6)), degs_2_6[sort_deg], color='blue')\n",
    "plt.xlabel('Sorted Nodes', fontsize=22)\n",
    "plt.ylabel('Degree', fontsize=22)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.ylim(0,380)\n",
    "plt.title(\"First Week\", fontsize=24)\n",
    "plt.grid(False)\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "plt.scatter(np.arange(n_tot_station), degs_9_13[sort_deg], color='blue')\n",
    "plt.xlabel('Sorted Nodes', fontsize=22)\n",
    "plt.ylabel('Degree', fontsize=22)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.ylim(0,380)\n",
    "plt.title(\"Second Week\", fontsize=24)\n",
    "plt.grid(False)\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "plt.scatter(np.arange(n_tot_station), degs_16_20[sort_deg], color='blue')\n",
    "plt.xlabel('Sorted Nodes', fontsize=22)\n",
    "plt.ylabel('Degree', fontsize=22)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.title(\"Third Week\", fontsize=24)\n",
    "plt.ylim(0,380)\n",
    "plt.grid(False)\n",
    "\n",
    "plt.subplot(2,2,4)\n",
    "plt.scatter(np.arange(n_tot_station), degs_23_27[sort_deg], color='blue')\n",
    "plt.xlabel('Sorted Nodes', fontsize=22)\n",
    "plt.ylabel('Degree', fontsize=22)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.title(\"Fourth Week\", fontsize=24)\n",
    "plt.ylim(0,380)\n",
    "plt.grid(False)\n",
    "\n",
    "plt.savefig('./plots/degree_ordered.eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute Flows Between Clusters and Regions #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "queens = []\n",
    "for k,v in part_2_6.items():\n",
    "    if part_2_6[k]==4 and part_23_27.get(k) == 2:\n",
    "        queens.append(k)\n",
    "\n",
    "brooklyn = []\n",
    "for k,v in part_2_6.items():\n",
    "    if part_2_6[k]==2 and part_23_27.get(k) == 2:\n",
    "        brooklyn.append(k)\n",
    "\n",
    "\n",
    "astoria = []\n",
    "for k,v in part_2_6.items():\n",
    "    if part_2_6[k]==4 and part_23_27.get(k) == 4:\n",
    "        astoria.append(k)\n",
    "\n",
    "manhattan = []\n",
    "for k,v in part_2_6.items():\n",
    "    man = {0,1,3}\n",
    "    if part_2_6[k] in man and part_23_27.get(k) in man:\n",
    "        manhattan.append(k)\n",
    "        \n",
    "blue_last = []\n",
    "for k,v in part_23_27.items():\n",
    "    if part_23_27.get(k) == 1:\n",
    "        blue_last.append(k)\n",
    "\n",
    "yellow_last = []\n",
    "for k,v in part_23_27.items():\n",
    "    if part_23_27.get(k) == 0:\n",
    "        yellow_last.append(k)\n",
    "\n",
    "blue_first = []\n",
    "for k,v in part_2_6.items():\n",
    "    if part_2_6.get(k) == 1:\n",
    "        blue_first.append(k)\n",
    "\n",
    "gray_first = []\n",
    "for k,v in part_2_6.items():\n",
    "    if part_2_6.get(k) == 0:\n",
    "        gray_first.append(k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inn_1, out_1  = compute_flows(g_2_6, queens)\n",
    "print(\"{:20s}   {:7d}     {:7d}\".format(\"First Flow Queens\", int(inn_1), int(out_1) ))\n",
    "\n",
    "inn_4, out_4  = compute_flows(g_23_27, queens)\n",
    "print(\"{:20s}   {:7d}     {:7d}\".format(\"Last Flow Queens\", int(inn_4), int(out_4) ))\n",
    "\n",
    "print(\"{:20s}   {:2.5f}     {:2.5f}\".format(\"Ratio\", inn_4/inn_1, out_4/out_1))\n",
    "print(\"\\n\")\n",
    "\n",
    "inn_1, out_1  = compute_flows(g_2_6, brooklyn)\n",
    "print(\"{:20s}   {:7d}     {:7d}\".format(\"First Flow brooklyn\", int(inn_1), int(out_1) ))\n",
    "\n",
    "inn_4, out_4  = compute_flows(g_23_27, brooklyn)\n",
    "print(\"{:20s}   {:7d}     {:7d}\".format(\"Last Flow brooklyn\", int(inn_4), int(out_4) ))\n",
    "\n",
    "print(\"{:20s}   {:2.5f}     {:2.5f}\".format(\"Ratio\", inn_4/inn_1, out_4/out_1))\n",
    "print(\"\\n\")\n",
    "\n",
    "inn_1, out_1  = compute_flows(g_2_6, manhattan)\n",
    "print(\"{:20s}   {:7d}     {:7d}\".format(\"First Flow manhattan\", int(inn_1), int(out_1) ))\n",
    "\n",
    "inn_4, out_4  = compute_flows(g_23_27, manhattan)\n",
    "print(\"{:20s}   {:7d}     {:7d}\".format(\"Last Flow manhattan\", int(inn_4), int(out_4) ))\n",
    "\n",
    "print(\"{:20s}   {:2.5f}     {:2.5f}\".format(\"Ratio\", inn_4/inn_1, out_4/out_1))\n",
    "print(\"\\n\")\n",
    "\n",
    "inn_1, out_1  = compute_flows(g_2_6, astoria)\n",
    "print(\"{:20s}   {:7d}     {:7d}\".format(\"First Flow astoria\", int(inn_1), int(out_1) ))\n",
    "\n",
    "inn_4, out_4  = compute_flows(g_23_27, astoria)\n",
    "print(\"{:20s}   {:7d}     {:7d}\".format(\"Last Flow astoria\", int(inn_4), int(out_4) ))\n",
    "\n",
    "print(\"{:20s}   {:2.5f}     {:2.5f}\".format(\"Ratio\", inn_4/inn_1, out_4/out_1))\n",
    "print(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flow_1 = compute_flows_between(g_2_6, queens, astoria)\n",
    "flow_4 = compute_flows_between(g_23_27, queens, astoria)\n",
    "print(\"{:25s}   {:7d}     {:7d}\".format(\"Flows Queens/Astoria\", int(flow_1), int(flow_4) ))\n",
    "\n",
    "print(\"{:25s}   {:2.5f}\".format(\"Ratio\", flow_4/flow_1))\n",
    "print(\"\\n\")\n",
    "\n",
    "flow_1 = compute_flows_between(g_2_6, queens, brooklyn)\n",
    "flow_4 = compute_flows_between(g_23_27, queens, brooklyn)\n",
    "print(\"{:25s}   {:7d}     {:7d}\".format(\"Flows Queens/Brooklyn\", int(flow_1), int(flow_4) ))\n",
    "\n",
    "print(\"{:25s}   {:2.5f}\".format(\"Ratio\", flow_4/flow_1))\n",
    "print(\"\\n\")\n",
    "\n",
    "flow_1 = compute_flows_between(g_2_6, astoria, manhattan)\n",
    "flow_4 = compute_flows_between(g_23_27, astoria, manhattan)\n",
    "print(\"{:25s}   {:7d}     {:7d}\".format(\"Flows Astoria/Manhattan\", int(flow_1), int(flow_4) ))\n",
    "\n",
    "print(\"{:25s}   {:2.5f}\".format(\"Ratio\", flow_4/flow_1))\n",
    "print(\"\\n\")\n",
    "\n",
    "flow_1 = compute_flows_between(g_2_6, queens, manhattan)\n",
    "flow_4 = compute_flows_between(g_23_27, queens, manhattan)\n",
    "print(\"{:25s}   {:7d}     {:7d}\".format(\"Flows Queens/Manhattan\", int(flow_1), int(flow_4) ))\n",
    "\n",
    "print(\"{:25s}   {:2.5f}\".format(\"Ratio\", flow_4/flow_1))\n",
    "print(\"\\n\")\n",
    "\n",
    "flow_1 = compute_flows_between(g_2_6, brooklyn, manhattan)\n",
    "flow_4 = compute_flows_between(g_23_27, brooklyn, manhattan)\n",
    "print(\"{:25s}   {:7d}     {:7d}\".format(\"Flows Brooklyn/Manhattan\", int(flow_1), int(flow_4) ))\n",
    "\n",
    "print(\"{:25s}   {:2.5f}\".format(\"Ratio\", flow_4/flow_1))\n",
    "print(\"\\n\")\n",
    "\n",
    "flow_1 = compute_flows_between(g_2_6, blue_last, yellow_last)\n",
    "flow_4 = compute_flows_between(g_23_27, blue_last, yellow_last)\n",
    "print(\"{:25s}   {:7d}     {:7d}\".format(\"Flows Blue/Yellow\", int(flow_1), int(flow_4) ))\n",
    "\n",
    "print(\"{:25s}   {:2.5f}\".format(\"Ratio\", flow_4/flow_1))\n",
    "print(\"\\n\")\n",
    "\n",
    "flow_1 = compute_flows_between(g_2_6, blue_first, gray_first)\n",
    "flow_4 = compute_flows_between(g_23_27, blue_first, gray_first)\n",
    "print(\"{:25s}   {:7d}     {:7d}\".format(\"Flows Blue/Gray\", int(flow_1), int(flow_4) ))\n",
    "\n",
    "print(\"{:25s}   {:2.5f}\".format(\"Ratio\", flow_4/flow_1))\n",
    "print(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Compare flows during first and last week </h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ORDER: UM LM BR MM AS WB\n",
    "cluster_sets = [set() for e in range(1 + len(set(part_2_6.values())))]\n",
    "out = []\n",
    "for k,v in part_2_6.items():\n",
    "    if k in set(list(g_23_27.nodes())):\n",
    "        cluster_sets[v].add(k)\n",
    "    else:\n",
    "        out.append(k)\n",
    "\n",
    "for e in cluster_sets[4]:\n",
    "    if part_2_6[e] == 4 and part_23_27.get(e) == 2:\n",
    "        cluster_sets[5].add(e)\n",
    "\n",
    "for e in cluster_sets[5]:\n",
    "    cluster_sets[4].remove(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First Week\n",
    "\n",
    "flows_vec_first = []\n",
    "for e in cluster_sets:\n",
    "    inn, out = compute_flows(g_2_6, e)\n",
    "    flows_vec_first.append(inn+out)\n",
    "    \n",
    "flows_vec_first = np.array(flows_vec_first)\n",
    "\n",
    "flow_matrix_first = np.zeros((len(flows_vec_first), len(flows_vec_first)))\n",
    "for i in range(flow_matrix_first.shape[0]):\n",
    "    for j in range(flow_matrix_first.shape[0]):\n",
    "        flow_matrix_first[i,j] = compute_flows_between(g_2_6, cluster_sets[i], cluster_sets[j])\n",
    "\n",
    "flow_matrix_first_ratio = flow_matrix_first / flows_vec_first[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Last Week\n",
    "\n",
    "flows_vec_second = []\n",
    "for e in cluster_sets:\n",
    "    inn, out = compute_flows(g_23_27, e)\n",
    "    flows_vec_second.append(inn+out)\n",
    "    \n",
    "flows_vec_second = np.array(flows_vec_second)\n",
    "\n",
    "flow_matrix_second = np.zeros((len(flows_vec_second), len(flows_vec_second)))\n",
    "for i in range(flow_matrix_second.shape[0]):\n",
    "    for j in range(flow_matrix_second.shape[0]):\n",
    "        flow_matrix_second[i,j] = compute_flows_between(g_23_27, cluster_sets[i], cluster_sets[j])\n",
    "\n",
    "flow_matrix_second_ratio = flow_matrix_second / flows_vec_second[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.colors as colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = ['UM', 'LM', 'MM', 'BR', 'WB', 'AS']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Relative Flows among regions in first and second week \n",
    "\n",
    "plt.figure(figsize=(18, 12))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(flow_matrix_first_ratio, cmap=cmap, interpolation='nearest', norm=colors.LogNorm(vmin=0.01, vmax=1))\n",
    "plt.grid(False)\n",
    "plt.colorbar()\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n",
    "\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(flow_matrix_second_ratio, cmap=cmap, interpolation='nearest', norm=colors.LogNorm(vmin=0.01, vmax=1))\n",
    "plt.grid(False)\n",
    "plt.colorbar()\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat = flow_matrix_second / flow_matrix_first\n",
    "\n",
    "mat[[2, 3]] = mat[[3, 2]]\n",
    "mat[:,[2, 3]] = mat[:,[3, 2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ratio between \n",
    "\n",
    "plt.figure(figsize=(14,12))\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(mat, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n",
    "_ = plt.colorbar().ax.tick_params(labelsize=22)\n",
    "\n",
    "plt.savefig('plots/heatmap_flows.eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3> Analyze Matrices </h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in mat:\n",
    "    s = ''\n",
    "    for i in e:\n",
    "        s += str(np.round(i, 4))\n",
    "        s += '  '\n",
    "    print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in flow_matrix_first_ratio:\n",
    "    s = ''\n",
    "    for i in e:\n",
    "        s += str(np.round(i, 4))\n",
    "        s += '  '\n",
    "    print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in flow_matrix_second_ratio:\n",
    "    s = ''\n",
    "    for i in e:\n",
    "        s += str(np.round(i, 4))\n",
    "        s += '  '\n",
    "    print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Comparison First week/weekend with last week\n",
    "flows_vec_weekend = []\n",
    "for e in cluster_sets:\n",
    "    inn, out = compute_flows(g_78, e)\n",
    "    flows_vec_weekend.append(inn+out)\n",
    "    \n",
    "flows_vec_weekend = np.array(flows_vec_weekend)\n",
    "\n",
    "flow_matrix_weekend = np.zeros((len(flows_vec_weekend), len(flows_vec_weekend)))\n",
    "for i in range(flow_matrix_weekend.shape[0]):\n",
    "    for j in range(flow_matrix_weekend.shape[0]):\n",
    "        flow_matrix_weekend[i,j] = compute_flows_between(g_78, cluster_sets[i], cluster_sets[j])\n",
    "\n",
    "flow_matrix_weekend_ratio = flow_matrix_weekend / flows_vec_weekend[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(18, 12))\n",
    "\n",
    "plt.subplot(2,2,1)\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(flow_matrix_weekend_ratio, cmap=cmap, interpolation='nearest', norm=colors.LogNorm(vmin=0.01, vmax=1))\n",
    "plt.grid(False)\n",
    "plt.colorbar()\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n",
    "\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(flow_matrix_first_ratio, cmap=cmap, interpolation='nearest', norm=colors.LogNorm(vmin=0.01, vmax=1))\n",
    "plt.grid(False)\n",
    "plt.colorbar()\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(flow_matrix_weekend_ratio, cmap=cmap, interpolation='nearest', norm=colors.LogNorm(vmin=0.01, vmax=1))\n",
    "plt.grid(False)\n",
    "plt.colorbar()\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n",
    "\n",
    "\n",
    "plt.subplot(2,2,4)\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(flow_matrix_second_ratio, cmap=cmap, interpolation='nearest', norm=colors.LogNorm(vmin=0.01, vmax=1))\n",
    "plt.grid(False)\n",
    "plt.colorbar()\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat_first_weekend = flow_matrix_first / flow_matrix_weekend\n",
    "\n",
    "mat_first_weekend[[2, 3]] = mat_first_weekend[[3, 2]]\n",
    "mat_first_weekend[:,[2, 3]] = mat_first_weekend[:,[3, 2]]\n",
    "\n",
    "mat_last_weekend = flow_matrix_second / flow_matrix_weekend\n",
    "\n",
    "mat_last_weekend[[2, 3]] = mat_last_weekend[[3, 2]]\n",
    "mat_last_weekend[:,[2, 3]] = mat_last_weekend[:,[3, 2]]\n",
    "\n",
    "mat_weekend = flow_matrix_weekend / flow_matrix_weekend\n",
    "\n",
    "mat_weekend[[2, 3]] = mat_weekend[[3, 2]]\n",
    "mat_weekend[:,[2, 3]] = mat_weekend[:,[3, 2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,12))\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(mat_first_weekend, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n",
    "_ = plt.colorbar().ax.tick_params(labelsize=22)\n",
    "\n",
    "# plt.savefig('plots/heatmap_flows.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,12))\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(mat_last_weekend, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n",
    "_ = plt.colorbar().ax.tick_params(labelsize=22)\n",
    "\n",
    "# plt.savefig('plots/heatmap_flows.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,12))\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(mat_weekend, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n",
    "_ = plt.colorbar().ax.tick_params(labelsize=22)\n",
    "\n",
    "# plt.savefig('plots/heatmap_flows.eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions on Graphs #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_2_6 = np.sum(adj_2_6_tot, axis=0)\n",
    "f_9_13 = np.sum(adj_9_13_tot, axis=0)\n",
    "f_16_20 = np.sum(adj_16_20_tot, axis=0)\n",
    "f_23_27 = np.sum(adj_23_27_tot, axis=0)\n",
    "\n",
    "rif = f_2_6 + 1\n",
    "\n",
    "r_2_6 = rif / rif\n",
    "r_9_13 = f_9_13 / rif\n",
    "r_16_20 = f_16_20 / rif\n",
    "r_23_27 = f_23_27 / rif\n",
    "\n",
    "rif = degs_2_6 + 1\n",
    "\n",
    "d_2_6 = rif / rif\n",
    "d_9_13 = degs_9_13 / rif\n",
    "d_16_20 = degs_16_20 / rif\n",
    "d_23_37 = degs_23_27 / rif\n",
    "\n",
    "f_2_6 = np.delete(f_2_6, isolated_nodes_2_6)\n",
    "f_9_13 = np.delete(f_9_13, isolated_nodes_9_13)\n",
    "f_16_20 = np.delete(f_16_20, isolated_nodes_16_20)\n",
    "f_23_27 = np.delete(f_23_27, isolated_nodes_23_27)\n",
    "\n",
    "r_2_6 = np.delete(r_2_6, isolated_nodes_2_6)\n",
    "r_9_13 = np.delete(r_9_13, isolated_nodes_9_13)\n",
    "r_16_20 = np.delete(r_16_20, isolated_nodes_16_20)\n",
    "r_23_27 = np.delete(r_23_27, isolated_nodes_23_27)\n",
    "\n",
    "d_2_6 = np.delete(d_2_6, isolated_nodes_2_6)\n",
    "d_9_13 = np.delete(d_9_13, isolated_nodes_9_13)\n",
    "d_16_20 = np.delete(d_16_20, isolated_nodes_16_20)\n",
    "d_23_27 = np.delete(d_23_37, isolated_nodes_23_27)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clusters #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,18))\n",
    "plt.subplot(2,2,1)\n",
    "\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('First Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_2_6, locations, node_size=40, node_color=vec_col_2_6)\n",
    "nx.draw_networkx_edges(g_2_6, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Second Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_9_13, locations, node_size=40, node_color=vec_col_9_13)\n",
    "nx.draw_networkx_edges(g_9_13, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Third Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_16_20, locations, node_size=40, node_color=vec_col_16_20)\n",
    "nx.draw_networkx_edges(g_16_20, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,4)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Fourth Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_23_27, locations, node_size=40, node_color=vec_col_23_27)\n",
    "nx.draw_networkx_edges(g_23_27, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.savefig('plots/clusters.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weighted_coverage(g_2_6, part_2_6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weighted_coverage(g_9_13, part_9_13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weighted_coverage(g_16_20, part_16_20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weighted_coverage(g_23_27, part_23_27)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Graph over means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_2_6 = f_2_6 / np.mean(f_2_6)\n",
    "m_9_13 = f_9_13 / np.mean(f_9_13)\n",
    "m_16_20 = f_16_20 / np.mean(f_16_20)\n",
    "m_23_27 = f_23_27 / np.mean(f_23_27)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap= matplotlib.cm.get_cmap('jet')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.mean(f_2_6))\n",
    "print(np.mean(f_9_13))\n",
    "print(np.mean(f_16_20))\n",
    "print(np.mean(f_23_27))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.max(m_2_6))\n",
    "print(np.max(m_9_13))\n",
    "print(np.max(m_16_20))\n",
    "print(np.max(m_23_27))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = plt.hist(m_2_6)\n",
    "_ = plt.hist(m_9_13)\n",
    "_ = plt.hist(m_16_20)\n",
    "_ = plt.hist(m_23_27)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,15))\n",
    "mat = np.arange(0, 4, 0.4)\n",
    "label = np.array([round(x, 2) for x in mat])\n",
    "mat = np.expand_dims(mat, axis=0)\n",
    "plt.imshow(mat, origin=\"lower\", cmap=cmap, interpolation='nearest')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9], label[[0, 3, 6, 9]], fontsize=22)\n",
    "plt.set_cmap(cmap=cmap)\n",
    "\n",
    "plt.savefig('plots/means_cmap.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap= matplotlib.cm.get_cmap('jet')\n",
    "plt.figure(figsize=(10,18))\n",
    "plt.subplot(2,2,1)\n",
    "\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "nx.draw_networkx_nodes(g_2_6, locations, node_size=80, node_color=m_2_6, cmap=cmap, vmin=0, vmax=7)\n",
    "# nx.draw_networkx_edges(g_2_6, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "mplleaflet.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,18))\n",
    "plt.subplot(2,2,1)\n",
    "\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('First Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_2_6, locations, node_size=40, node_color=m_2_6, cmap=cmap, vmin=0, vmax=4)\n",
    "nx.draw_networkx_edges(g_2_6, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Second Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_9_13, locations, node_size=40, node_color=m_9_13, cmap=cmap, vmin=0, vmax=4)\n",
    "nx.draw_networkx_edges(g_9_13, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Third Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_16_20, locations, node_size=40, node_color=m_16_20, cmap=cmap, vmin=0, vmax=4)\n",
    "nx.draw_networkx_edges(g_16_20, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,4)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Fourth Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_23_27, locations, node_size=40, node_color=m_23_27, cmap=cmap, vmin=0, vmax=4)\n",
    "nx.draw_networkx_edges(g_23_27, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.savefig('plots/over_means.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(m_2_6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rev = dict(zip(list(id_index.values()), list(id_index.keys())))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[df['start station id'] == rev[117]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap= matplotlib.cm.get_cmap('jet')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,15))\n",
    "mat = np.arange(0, 1.5, 0.1)\n",
    "label = np.array([round(x, 2) for x in mat])\n",
    "mat = np.expand_dims(mat, axis=0)\n",
    "plt.imshow(mat, origin=\"lower\", cmap=cmap, interpolation='nearest')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9, 12], label[[0, 3, 6, 9, 12]], fontsize=22)\n",
    "plt.set_cmap(cmap=cmap)\n",
    "\n",
    "plt.savefig('plots/relative_cmap.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,18))\n",
    "plt.subplot(2,2,1)\n",
    "\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('First Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_2_6, locations, node_size=40, node_color=r_2_6, cmap=cmap, vmin=0, vmax=1.5)\n",
    "nx.draw_networkx_edges(g_2_6, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Second Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_9_13, locations, node_size=40, node_color=r_9_13, cmap=cmap, vmin=0, vmax=1.5)\n",
    "nx.draw_networkx_edges(g_9_13, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Third Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_16_20, locations, node_size=40, node_color=r_16_20, cmap=cmap, vmin=0, vmax=1.5)\n",
    "nx.draw_networkx_edges(g_16_20, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,4)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Fourth Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_23_27, locations, node_size=40, node_color=r_23_27, cmap=cmap, vmin=0, vmax=1.5)\n",
    "nx.draw_networkx_edges(g_23_27, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.savefig('plots/ratio_first_week.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,18))\n",
    "plt.subplot(2,2,1)\n",
    "\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('First Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_2_6, locations, node_size=40, node_color=d_2_6, cmap=cmap, vmin=0, vmax=1.5)\n",
    "nx.draw_networkx_edges(g_2_6, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,2)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Second Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_9_13, locations, node_size=40, node_color=d_9_13, cmap=cmap, vmin=0, vmax=1.5)\n",
    "nx.draw_networkx_edges(g_9_13, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,3)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Third Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_16_20, locations, node_size=40, node_color=d_16_20, cmap=cmap, vmin=0, vmax=1.5)\n",
    "nx.draw_networkx_edges(g_16_20, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(2,2,4)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('Fourth Week', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_23_27, locations, node_size=40, node_color=d_23_27, cmap=cmap, vmin=0, vmax=1.5)\n",
    "nx.draw_networkx_edges(g_23_27, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.savefig('plots/ratio_deg_first_week.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(899):\n",
    "    adj_23_27_tot[i,i] = 0\n",
    "    adj_78_tot[i,i] = 0\n",
    "    adj_2_6_tot[i,i] = 0\n",
    "adj_23_27_tot /= 2\n",
    "adj_2_6_tot /= 2\n",
    "adj_78_tot /= 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flow_78 = np.sum(adj_78_tot, axis=0) \n",
    "flow_78 /= np.mean(flow_78)\n",
    "plt.plot(flow_78)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flow_26 = np.sum(adj_2_6_tot, axis=0) \n",
    "flow_26 /= np.mean(flow_26)\n",
    "plt.plot(flow_26)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flow_2327 = np.sum(adj_23_27_tot, axis=0) \n",
    "flow_2327 /= np.mean(flow_2327)\n",
    "plt.plot(flow_2327)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.norm(flow_2327 - flow_78, ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.norm(flow_26 - flow_78, ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.norm(flow_2327 - flow_26, ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.norm(flow_78, ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.norm(flow_2327, ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linalg.norm(flow_26, ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_flow_78 = flow_78+np.ones(899)*0.001 / np.sum(flow_78)\n",
    "n_flow_26 = flow_26 / np.sum(flow_26)\n",
    "n_flow_2327 = flow_2327 / np.sum(flow_2327)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp.stats.entropy(n_flow_26, n_flow_2327)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First Weekend #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9 -> Monday\n",
    "# Choose day \n",
    "days = [7, 8]\n",
    "\n",
    "# Find all connections with weights\n",
    "adj_78_tot = np.zeros((n_tot_station, n_tot_station))\n",
    "for day in days:\n",
    "    for index, row in df.iloc[d[day]:d[day+1]].iterrows():\n",
    "        adj_78_tot[id_index[row['start station id']], id_index[row['end station id']]] += 1\n",
    "        adj_78_tot[id_index[row['end station id']], id_index[row['start station id']]] += 1\n",
    "    print('Day {} loaded...'.format(day))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ths = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adj_78 = adj_78_tot.copy()\n",
    "adj_78[adj_78 <= ths] = 0\n",
    "\"\"\"adj_23_27[adj_23_27 >= ths] = 1\"\"\"\n",
    "print(\"Unweighted number of rides : {}\".format(np.sum(adj_78) / 2))\n",
    "\n",
    "# Remove Self-loops\n",
    "for i in range(n_tot_station):\n",
    "    adj_78[i, i] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g_78 = nx.from_numpy_matrix(adj_78)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolated_nodes_78 = []\n",
    "for e in nx.connected_components(g_78):\n",
    "    if len(e) < 3:\n",
    "        isolated_nodes_78.append(e.pop())\n",
    "\n",
    "isolated_nodes_78.extend(out_nodes)\n",
    "isolated_nodes_78 = list(set(isolated_nodes_78))\n",
    "g_78.remove_nodes_from(isolated_nodes_78)\n",
    "\n",
    "for e in nx.connected_components(g_78):\n",
    "    if len(e) < 3:\n",
    "        print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_mod_78 = 0\n",
    "part_78 = None\n",
    "\n",
    "for i in range(20):\n",
    "\n",
    "    partition = community.best_partition(g_78, weight='weight')\n",
    "    c_g = [set() for e in range(len(partition))]\n",
    "    for k, v in partition.items():\n",
    "        c_g[v].add(k)\n",
    "    mod = community.modularity(partition, g_78, weight='weight')\n",
    "    if max_mod_78 < mod:\n",
    "        max_mod_78 = mod\n",
    "        part_78 = partition\n",
    "\n",
    "\n",
    "print(\"weighted coverage : {}\".format(weighted_coverage(g_78, part_78)))\n",
    "\n",
    "color_cluster = ['yellow', 'blue', 'red', 'gray', 'green']\n",
    "vec_col_78 = [color_cluster[e] for e in part_78.values()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat = flow_matrix_second / flow_matrix_first\n",
    "\n",
    "mat[[2, 3]] = mat[[3, 2]]\n",
    "mat[:,[2, 3]] = mat[:,[3, 2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 12))\n",
    "\n",
    "nx.draw_networkx_nodes(g_78, locations, node_size=40, node_color=vec_col_78)\n",
    "nx.draw_networkx_edges(g_78, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.axis('off')\n",
    "\n",
    "# plt.savefig('plots/weighted_cluster_16_20.jpeg')\n",
    "\n",
    "#mplleaflet.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_78 = np.sum(adj_78_tot, axis=0)\n",
    "f_78 = np.delete(f_78, isolated_nodes_78)\n",
    "f_78 = f_78 / np.mean(f_78)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,15))\n",
    "mat = np.arange(0, 4, 0.4)\n",
    "label = np.array([round(x, 2) for x in mat])\n",
    "mat = np.expand_dims(mat, axis=0)\n",
    "plt.imshow(mat, origin=\"lower\", cmap=cmap, interpolation='nearest')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9], label[[0, 3, 6, 9]], fontsize=22)\n",
    "plt.set_cmap(cmap=cmap)\n",
    "\n",
    "plt.savefig('plots/means_cmap.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 12))\n",
    "\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('First Weekend', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_78, locations, node_size=40, node_color=f_78, cmap=cmap, vmin=0, vmax=4)\n",
    "nx.draw_networkx_edges(g_78, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.savefig('plots/first_weekend_flow.jpg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Second Weekdend #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 9 -> Monday\n",
    "# Choose day \n",
    "days = [14, 15]\n",
    "\n",
    "# Find all connections with weights\n",
    "adj_1415_tot = np.zeros((n_tot_station, n_tot_station))\n",
    "for day in days:\n",
    "    for index, row in df.iloc[d[day]:d[day+1]].iterrows():\n",
    "        adj_1415_tot[id_index[row['start station id']], id_index[row['end station id']]] += 1\n",
    "        adj_1415_tot[id_index[row['end station id']], id_index[row['start station id']]] += 1\n",
    "    print('Day {} loaded...'.format(day))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adj_1415 = adj_1415_tot.copy()\n",
    "adj_1415[adj_1415 <= ths] = 0\n",
    "\"\"\"adj_23_27[adj_23_27 >= ths] = 1\"\"\"\n",
    "print(\"Unweighted number of rides : {}\".format(np.sum(adj_78) / 2))\n",
    "\n",
    "# Remove Self-loops\n",
    "for i in range(n_tot_station):\n",
    "    adj_1415[i, i] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g_1415 = nx.from_numpy_matrix(adj_1415)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolated_nodes_1415 = []\n",
    "for e in nx.connected_components(g_1415):\n",
    "    if len(e) < 3:\n",
    "        isolated_nodes_1415.append(e.pop())\n",
    "\n",
    "isolated_nodes_1415.extend(out_nodes)\n",
    "isolated_nodes_1415 = list(set(isolated_nodes_1415))\n",
    "g_1415.remove_nodes_from(isolated_nodes_1415)\n",
    "\n",
    "for e in nx.connected_components(g_1415):\n",
    "    if len(e) < 3:\n",
    "        print(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_mod_1415 = 0\n",
    "part_1415 = None\n",
    "\n",
    "for i in range(20):\n",
    "\n",
    "    partition = community.best_partition(g_1415, weight='weight')\n",
    "    c_g = [set() for e in range(len(partition))]\n",
    "    for k, v in partition.items():\n",
    "        c_g[v].add(k)\n",
    "    mod = community.modularity(partition, g_1415, weight='weight')\n",
    "    if max_mod_1415 < mod:\n",
    "        max_mod_1415 = mod\n",
    "        part_1415 = partition\n",
    "\n",
    "\n",
    "print(\"weighted coverage : {}\".format(weighted_coverage(g_1415, part_1415)))\n",
    "\n",
    "color_cluster = ['yellow', 'blue', 'red', 'gray', 'green']\n",
    "vec_col_1415 = [color_cluster[e] for e in part_1415.values()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 12))\n",
    "\n",
    "nx.draw_networkx_nodes(g_1415, locations, node_size=40, node_color=vec_col_1415)\n",
    "nx.draw_networkx_edges(g_1415, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.axis('off')\n",
    "\n",
    "# plt.savefig('plots/weighted_cluster_16_20.jpeg')\n",
    "\n",
    "#mplleaflet.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_1415 = np.sum(adj_1415_tot, axis=0)\n",
    "f_1415 = np.delete(f_1415, isolated_nodes_1415)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 12))\n",
    "\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "plt.title('First Weekend', fontsize=22)\n",
    "nx.draw_networkx_nodes(g_1415, locations, node_size=40, node_color=f_1415, cmap=cmap)\n",
    "nx.draw_networkx_edges(g_1415, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.savefig('plots/second_weekend_flow.jpg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DiGraph Analysis #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIRST WEEK\n",
    "# Choose day \n",
    "days = [2, 3, 4, 5, 6]\n",
    "\n",
    "# Find all connections with weights\n",
    "dadj_2_6_tot = np.zeros((n_tot_station, n_tot_station))\n",
    "for day in days:\n",
    "    for index, row in df.iloc[d[day]:d[day+1]].iterrows():\n",
    "        dadj_2_6_tot[id_index[row['start station id']], id_index[row['end station id']]] += 1\n",
    "    print('Day {} loaded...'.format(day))\n",
    "    \n",
    "\n",
    "dadj_2_6 = dadj_2_6_tot.copy()\n",
    "dadj_2_6[dadj_2_6 <= ths] = 0\n",
    "\"\"\"adj_23_27[adj_23_27 >= ths] = 1\"\"\"\n",
    "\n",
    "\n",
    "# Remove Self-loops\n",
    "for i in range(n_tot_station):\n",
    "    dadj_2_6[i, i] = 0\n",
    "    \n",
    "print(\"Weighted number of rides : {}\".format(np.sum(dadj_2_6)))\n",
    "\n",
    "dg_2_6 = nx.from_numpy_matrix(dadj_2_6, create_using=nx.DiGraph)\n",
    "\n",
    "dg_2_6.remove_nodes_from(out_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# LAST WEEK\n",
    "# Choose day \n",
    "days = [23, 24, 25, 26, 27]\n",
    "\n",
    "# Find all connections with weights\n",
    "dadj_23_27_tot = np.zeros((n_tot_station, n_tot_station))\n",
    "for day in days:\n",
    "    for index, row in df.iloc[d[day]:d[day+1]].iterrows():\n",
    "        dadj_23_27_tot[id_index[row['start station id']], id_index[row['end station id']]] += 1\n",
    "    print('Day {} loaded...'.format(day))\n",
    "    \n",
    "\n",
    "dadj_23_27 = dadj_23_27_tot.copy()\n",
    "dadj_23_27[dadj_23_27 <= ths] = 0\n",
    "\"\"\"adj_23_27[adj_23_27 >= ths] = 1\"\"\"\n",
    "\n",
    "\n",
    "# Remove Self-loops\n",
    "for i in range(n_tot_station):\n",
    "    dadj_23_27[i, i] = 0\n",
    "    \n",
    "print(\"Weighted number of rides : {}\".format(np.sum(dadj_23_27)))\n",
    "\n",
    "dg_23_27 = nx.from_numpy_matrix(dadj_23_27, create_using=nx.DiGraph)\n",
    "\n",
    "dg_23_27.remove_nodes_from(out_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First Week\n",
    "\n",
    "d_flows_vec_first = []\n",
    "for e in cluster_sets:\n",
    "    inn, out = compute_flows(dg_2_6, e)\n",
    "    d_flows_vec_first.append(inn+out)\n",
    "    \n",
    "d_flows_vec_first = np.array(d_flows_vec_first)\n",
    "\n",
    "d_flow_matrix_first = np.zeros((len(d_flows_vec_first), len(d_flows_vec_first)))\n",
    "for i in range(d_flow_matrix_first.shape[0]):\n",
    "    for j in range(d_flow_matrix_first.shape[0]):\n",
    "        d_flow_matrix_first[i,j] = compute_flows_between(dg_2_6, cluster_sets[i], cluster_sets[j])\n",
    "\n",
    "d_flow_matrix_first_ratio = d_flow_matrix_first / d_flows_vec_first[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Last Week\n",
    "\n",
    "d_flows_vec_last = []\n",
    "for e in cluster_sets:\n",
    "    inn, out = compute_flows(dg_23_27, e)\n",
    "    d_flows_vec_last.append(inn+out)\n",
    "    \n",
    "d_flows_vec_last = np.array(d_flows_vec_last)\n",
    "\n",
    "d_flow_matrix_last = np.zeros((len(d_flows_vec_last), len(d_flows_vec_last)))\n",
    "for i in range(d_flow_matrix_last.shape[0]):\n",
    "    for j in range(d_flow_matrix_last.shape[0]):\n",
    "        d_flow_matrix_last[i,j] = compute_flows_between(dg_23_27, cluster_sets[i], cluster_sets[j])\n",
    "\n",
    "d_flow_matrix_last_ratio = d_flow_matrix_last / d_flows_vec_last[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat = d_flow_matrix_last / d_flow_matrix_first\n",
    "\n",
    "mat[[2, 3]] = mat[[3, 2]]\n",
    "mat[:,[2, 3]] = mat[:,[3, 2]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,12))\n",
    "cmap = plt.get_cmap('hot')\n",
    "plt.imshow(mat, cmap=cmap, interpolation='nearest', vmin=0, vmax=1)\n",
    "plt.grid(False)\n",
    "\n",
    "_= plt.xticks(np.arange(6), names, fontsize=22)\n",
    "_= plt.yticks(np.arange(6), names, fontsize=22)\n",
    "_ = plt.colorbar().ax.tick_params(labelsize=22)\n",
    "\n",
    "plt.savefig('plots/d_heatmap_flows.eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1) Temporal Analysis # "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weekly days #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First Week \n",
    "hour_full_2_8 = []\n",
    "j = 1\n",
    "# plt.figure(figsize=(12,8))\n",
    "\n",
    "for i in [2, 3, 4, 5, 6, 7, 8]:\n",
    "    # plt.subplot(2,3,j)\n",
    "    hours = hour_rides(d[i], d[i+1], df)\n",
    "\n",
    "    h_rides = [len(v) for k,v in hours.items()]\n",
    "    hour_full_2_8.extend(h_rides)\n",
    "    # plt.plot(h_rides)\n",
    "    j += 1\n",
    "    1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Second Week\n",
    "hour_full_9_15 = []\n",
    "j = 1\n",
    "# plt.figure(figsize=(12,8))\n",
    "\n",
    "for i in [9, 10, 11, 12, 13, 14, 15]:\n",
    "    # plt.subplot(2,3,j)\n",
    "    hours = hour_rides(d[i], d[i+1], df)\n",
    "\n",
    "    h_rides = [len(v) for k,v in hours.items()]\n",
    "    hour_full_9_15.extend(h_rides)\n",
    "    # plt.plot(h_rides)\n",
    "    j += 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Third Week\n",
    "hour_full_16_22 = []\n",
    "j = 1\n",
    "# plt.figure(figsize=(12,8))\n",
    "\n",
    "for i in [16, 17, 18, 19, 20, 21, 22]:\n",
    "    # plt.subplot(2,3,j)\n",
    "    hours = hour_rides(d[i], d[i+1], df)\n",
    "\n",
    "    h_rides = [len(v) for k,v in hours.items()]\n",
    "    hour_full_16_22.extend(h_rides)\n",
    "    # plt.plot(h_rides)\n",
    "    j += 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fourth Week\n",
    "hour_full_23_29 = []\n",
    "j = 1\n",
    "# plt.figure(figsize=(12,8))\n",
    "\n",
    "for i in [23, 24, 25, 26, 27, 28, 29]:\n",
    "    # plt.subplot(2,3,j)\n",
    "    hours = hour_rides(d[i], d[i+1], df)\n",
    "\n",
    "    h_rides = [len(v) for k,v in hours.items()]\n",
    "    hour_full_23_29.extend(h_rides)\n",
    "    # plt.plot(h_rides)\n",
    "    j += 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_month = []\n",
    "full_month.extend(hour_full_2_8)\n",
    "full_month.extend(hour_full_9_15)\n",
    "full_month.extend(hour_full_16_22)\n",
    "full_month.extend(hour_full_23_29)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rain = [[24,30],[36,42], [108,114],[114,120],[270,276],[366,372],[408,414],[414,420],[510,516],[516,520],[564,570], [630, 636]]\n",
    "weekend = [[120, 167], [287,335], [455, 503], [623, 671]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('./data/python_variables/weekedns.pickle', 'wb') as file:\n",
    "    pickle.dump(weekend, file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(24,12))\n",
    "_ = plt.plot(full_month, color='black')\n",
    "_ = plt.plot(np.arange(120, 167), hour_full_2_8[120:], color='green')\n",
    "_ = plt.plot(np.arange(287, 335), hour_full_9_15[120:], color='green')\n",
    "_ = plt.plot(np.arange(455, 503), hour_full_16_22[120:], color='green')\n",
    "_ = plt.plot(np.arange(623, 671), hour_full_23_29[120:], color='green')\n",
    "\n",
    "\n",
    "_ = plt.plot(np.arange(24, 30), full_month[24:30], color='blue')\n",
    "_ = plt.plot(np.arange(36, 42), full_month[36:42], color='blue')\n",
    "_ = plt.plot(np.arange(108, 114), full_month[108:114], color='blue')\n",
    "_ = plt.plot(np.arange(114, 120), full_month[114:120], color='blue')\n",
    "_ = plt.plot(np.arange(270, 276), full_month[270:276], color='blue')\n",
    "_ = plt.plot(np.arange(366, 372), full_month[366:372], color='blue')\n",
    "_ = plt.plot(np.arange(408, 414), full_month[408:414], color='blue')\n",
    "_ = plt.plot(np.arange(414, 420), full_month[414:420], color='blue')\n",
    "_ = plt.plot(np.arange(510, 516), full_month[510:516], color='blue')\n",
    "_ = plt.plot(np.arange(516, 520), full_month[516:520], color='blue')\n",
    "_ = plt.plot(np.arange(564, 570), full_month[564:570], color='blue')\n",
    "_ = plt.plot(np.arange(630, 636), full_month[630:636], color='blue')\n",
    "\n",
    "\n",
    "_ = plt.xticks([0, 167, 335, 503], np.arange(1, 5), fontsize=18)\n",
    "\n",
    "_ = plt.xlabel('Hours', fontsize=24)\n",
    "_ = plt.ylabel('# of Rides', fontsize=24)\n",
    "_ = plt.yticks(fontsize=18)\n",
    "plt.grid(False)\n",
    "plt.savefig('./plots/hours_rides.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rain + Weekends, handcrafted\n",
    "\n",
    "plt.figure(figsize=(24,12))\n",
    "_ = plt.plot(full_month, color='black')\n",
    "_ = plt.plot(np.arange(120, 167), hour_full_2_8[120:], color='green')\n",
    "_ = plt.plot(np.arange(287, 335), hour_full_9_15[120:], color='green')\n",
    "_ = plt.plot(np.arange(455, 503), hour_full_16_22[120:], color='green')\n",
    "_ = plt.plot(np.arange(623, 671), hour_full_23_29[120:], color='green')\n",
    "\n",
    "\n",
    "_ = plt.plot(np.arange(24, 30), full_month[24:30], color='blue')\n",
    "_ = plt.plot(np.arange(36, 42), full_month[36:42], color='blue')\n",
    "_ = plt.plot(np.arange(108, 114), full_month[108:114], color='blue')\n",
    "_ = plt.plot(np.arange(114, 120), full_month[114:120], color='blue')\n",
    "_ = plt.plot(np.arange(270, 276), full_month[270:276], color='blue')\n",
    "_ = plt.plot(np.arange(366, 372), full_month[366:372], color='blue')\n",
    "_ = plt.plot(np.arange(408, 414), full_month[408:414], color='blue')\n",
    "_ = plt.plot(np.arange(414, 420), full_month[414:420], color='blue')\n",
    "_ = plt.plot(np.arange(510, 516), full_month[510:516], color='blue')\n",
    "_ = plt.plot(np.arange(516, 520), full_month[516:520], color='blue')\n",
    "_ = plt.plot(np.arange(564, 570), full_month[564:570], color='blue')\n",
    "_ = plt.plot(np.arange(630, 636), full_month[630:636], color='blue')\n",
    "\n",
    "\n",
    "_ = plt.xticks([0, 167, 335, 503], np.arange(1, 5), fontsize=18)\n",
    "\n",
    "_ = plt.xlabel('Hours', fontsize=24)\n",
    "_ = plt.ylabel('# of Rides', fontsize=24)\n",
    "_ = plt.yticks(fontsize=18)\n",
    "plt.grid(False)\n",
    "plt.savefig('./plots/hours_rides.eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trip Duration #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set last day\n",
    "d[32] = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "td = []\n",
    "days = np.arange(1, 32)\n",
    "out = []\n",
    "for day in days:\n",
    "    ti = []\n",
    "    for index, row in df.iloc[d[day]:d[day+1]].iterrows():\n",
    "        if row['tripduration'] < 3600:\n",
    "            ti.append(row['tripduration'] )\n",
    "        else:\n",
    "            out.append(row['tripduration'])\n",
    "    td.append(ti)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = [np.mean(e) for e in td]\n",
    "\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.plot(np.arange(1, 32), m, color='blue', marker='*')\n",
    "# _ = plt.xticks(np.arange(1, 32), np.arange(1, 32))\n",
    "_ = plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)\n",
    "plt.xlabel('Day', fontsize=24)\n",
    "plt.ylabel('Avg Trip Duration', fontsize=24)\n",
    "plt.grid(False)\n",
    "plt.savefig('./plots/time_duration_means.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first = []\n",
    "for e in [2,3,4,5,6]:\n",
    "    first.extend(td[e])\n",
    "first = np.array(first)\n",
    "m_f = np.mean(first)\n",
    "s_f = np.std(first)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "second = []\n",
    "for e in [9, 10, 11, 12, 13]:\n",
    "    second.extend(td[e])\n",
    "second = np.array(second)\n",
    "m_s = np.mean(second)\n",
    "s_s = np.std(second)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "third = []\n",
    "for e in [16, 17, 18, 19, 20]:\n",
    "    third.extend(td[e])\n",
    "third = np.array(third)\n",
    "m_t = np.mean(third)\n",
    "s_t = np.std(third)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = []\n",
    "for e in [23, 24, 25, 26, 27] :\n",
    "    q.extend(td[e])\n",
    "q = np.array(q)\n",
    "m_q = np.mean(q)\n",
    "s_q = np.std(q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "means = [m_f, m_s, m_t, m_q]\n",
    "ss = [s_f, s_s, s_t, s_q]\n",
    "x = np.arange(4)\n",
    "\n",
    "plt.errorbar(x, means, ss, linestyle='None', marker='^')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat_time = np.zeros((36, 31))\n",
    "i = 0\n",
    "for e in td:\n",
    "    for a in e:\n",
    "        mat_time[int(a/100),i] += 1\n",
    "    i += 1\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat_time /= np.sum(mat_time, axis=0)\n",
    "\n",
    "mat_time = np.cumsum(mat_time, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = matplotlib.cm.get_cmap('hot')\n",
    "plt.figure(figsize=(14,15))\n",
    "mat = np.arange(0, 1, 0.1)\n",
    "label = np.array([round(x, 2) for x in mat])\n",
    "mat = np.expand_dims(mat, axis=0)\n",
    "plt.imshow(mat, origin=\"lower\", cmap=cmap, interpolation='nearest')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9], label[[0, 3, 6, 9]], fontsize=22)\n",
    "plt.set_cmap(cmap=cmap)\n",
    "\n",
    "plt.savefig('plots/trip_duration_cmap.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat_time = np.flip(mat_time, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_mat = mat_time[12:, :]\n",
    "plot_mat = np.flip(plot_mat, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "med = [np.median(e) for e in td]\n",
    "plt.figure(figsize=(14, 12))\n",
    "plt.imshow(plot_mat, cmap='hot', origin='lower')\n",
    "plt.grid(False)\n",
    "_ = plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)\n",
    "plt.xlabel('Day', fontsize=24)\n",
    "plt.ylabel('Empirical CDF (1 hour)', fontsize=24)\n",
    "_ = plt.colorbar().ax.tick_params(labelsize=22)\n",
    "\n",
    "\n",
    "plt.plot(np.arange(0, 31), np.array(med)/100, color='blue', marker='*')\n",
    "_ = plt.xticks([0, 4, 9, 14, 19, 24, 30], [1, 5, 10, 15, 20, 25, 31], fontsize=18)\n",
    "_ = plt.yticks([5, 11, 17, 23], ['10 min', '20 min', '30 min', '40 min'], fontsize=18)\n",
    "plt.xlabel('Day', fontsize=24)\n",
    "plt.ylabel('Avg Trip Duration', fontsize=24)\n",
    "plt.grid(False)\n",
    "plt.savefig('./plots/cdf_timeduration.eps')\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(nx.algorithms.connectivity.edge_connectivity(dg_2_6, s=600, t=40))\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_23_27, s=600, t=40))\n",
    "print(\"\\n\")\n",
    "\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_2_6, s=600, t=776))\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_23_27, s=600, t=776))\n",
    "print(\"\\n\")\n",
    "\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_2_6, s=40, t=776))\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_23_27, s=40, t=776))\n",
    "print(\"\\n\")\n",
    "\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_2_6, s=40, t=524))\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_23_27, s=40, t=524))\n",
    "print(\"\\n\")\n",
    "\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_2_6, s=600, t=524))\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_23_27, s=600, t=524))\n",
    "print(\"\\n\")\n",
    "\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_2_6, s=776, t=524))\n",
    "print(nx.algorithms.connectivity.edge_connectivity(dg_23_27, s=776, t=524))\n",
    "print(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Different Ratio Graphs # "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_4 = f_4 / f_4\n",
    "r_8 = f_8 / f_4\n",
    "r_11 = f_11 / f_4\n",
    "r_15 = f_15 / f_4\n",
    "r_18 = f_18 / f_4\n",
    "r_22 = f_22 / f_4\n",
    "r_24 = f_24 / f_4\n",
    "r_29 = f_29 / f_4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.max(r_4))\n",
    "print(np.max(r_8))\n",
    "print(np.max(r_11))\n",
    "print(np.max(r_15))\n",
    "print(np.max(r_18))\n",
    "print(np.max(r_22))\n",
    "print(np.max(r_24))\n",
    "print(np.max(r_29))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.max(f_4))\n",
    "print(np.max(f_11))\n",
    "print(np.max(f_18))\n",
    "print(np.max(f_24))\n",
    "\n",
    "print(np.min(f_4))\n",
    "print(np.min(f_11))\n",
    "print(np.min(f_18))\n",
    "print(np.min(f_24))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = plt.cm.get_cmap('jet')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 10))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.grid(False)\n",
    "plt.axis('off')\n",
    "nx.draw_networkx_nodes(g_4, locations, node_size=40, node_color=f_4, cmap=cmap, vmin=0, vmax=1197)\n",
    "nx.draw_networkx_edges(g_4, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "\"\"\"plt.grid(False)\n",
    "plt.axis('off')\n",
    "nx.draw_networkx_nodes(g_8, locations, node_size=40, node_color=r_8, cmap=cmap, vmin=0, vmax=1.5)\n",
    "nx.draw_networkx_edges(g_8, locations, alpha=0.8, width=0.2, edge_color='black')\"\"\"\n",
    "\n",
    "mat = np.arange(0,1197, 100)\n",
    "mat = np.expand_dims(mat, axis=1)\n",
    "plt.imshow(mat, origin=\"lower\", cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.xticks([0], [''])\n",
    "_ = plt.yticks([0, 3, 6, 9, 11], mat[[0,3, 6, 9, 11],0])\n",
    "plt.set_cmap(cmap=cmap)\n",
    "\n",
    "_ = plt.savefig('./plots/day_4_abs.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 10))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.grid(False)\n",
    "plt.axis('off')\n",
    "nx.draw_networkx_nodes(g_11, locations, node_size=40, node_color=f_11, cmap=cmap, vmin=0, vmax=1197)\n",
    "nx.draw_networkx_edges(g_11, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "\"\"\"plt.grid(False)\n",
    "plt.axis('off')\n",
    "nx.draw_networkx_nodes(g_15, locations, node_size=40, node_color=r_15, cmap=cmap, vmin=0, vmax=1)\n",
    "nx.draw_networkx_edges(g_15, locations, alpha=0.8, width=0.2, edge_color='black')\"\"\"\n",
    "\n",
    "mat = np.arange(0,1197, 100)\n",
    "mat = np.expand_dims(mat, axis=1)\n",
    "plt.imshow(mat, origin=\"lower\", cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.xticks([0], [''])\n",
    "_ = plt.yticks([0, 3, 6, 9, 11], mat[[0,3, 6, 9, 11],0])\n",
    "plt.set_cmap(cmap=cmap)\n",
    "_ = plt.savefig('./plots/day_11_abs.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 10))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.grid(False)\n",
    "plt.axis('off')\n",
    "nx.draw_networkx_nodes(g_18, locations, node_size=40, node_color=f_18, cmap=cmap, vmin=0, vmax=1197)\n",
    "nx.draw_networkx_edges(g_18, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "\"\"\"plt.grid(False)\n",
    "plt.axis('off')\n",
    "nx.draw_networkx_nodes(g_22, locations, node_size=40, node_color=r_22, cmap=cmap, vmin=0, vmax=1)\n",
    "nx.draw_networkx_edges(g_22, locations, alpha=0.8, width=0.2, edge_color='black')\"\"\"\n",
    "\n",
    "mat = np.arange(0,1.5, 0.1)\n",
    "mat = np.expand_dims(mat, axis=1)\n",
    "plt.imshow(mat, origin=\"lower\", cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.xticks([0], [''])\n",
    "_ = plt.yticks([0, 3, 6, 9, 12], [0,0.3,0.6,0.9,1.2])\n",
    "plt.set_cmap(cmap=cmap)\n",
    "\n",
    "\n",
    "_ = plt.savefig('./plots/day_18_abs.jpeg')\n",
    "\n",
    "#mplleaflet.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 10))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.grid(False)\n",
    "plt.axis('off')\n",
    "nx.draw_networkx_nodes(g_24, locations, node_size=40, node_color=f_24, cmap=cmap, vmin=0, vmax=1197)\n",
    "nx.draw_networkx_edges(g_24, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "\"\"\"plt.grid(False)\n",
    "plt.axis('off')\n",
    "nx.draw_networkx_nodes(g_29, locations, node_size=40, node_color=r_29, cmap=cmap, vmin=0, vmax=1)\n",
    "nx.draw_networkx_edges(g_29, locations, alpha=0.8, width=0.2, edge_color='black')\"\"\"\n",
    "\n",
    "mat = np.arange(0,1.5, 0.1)\n",
    "mat = np.expand_dims(mat, axis=1)\n",
    "plt.imshow(mat, origin=\"lower\", cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.xticks([0], [''])\n",
    "_ = plt.yticks([0, 3, 6, 9, 12], [0,0.3,0.6,0.9,1.2])\n",
    "plt.set_cmap(cmap=cmap)\n",
    "\n",
    "_ = plt.savefig('./plots/day_24_abs.jpeg')\n",
    "\n",
    "# mplleaflet.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = [len(e) for k,e in hour_rides(d[3], d[4], df).items()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Important nodes first week\n",
    "s = set(np.argsort(f_2_6)[-200:])\n",
    "g_top_2_6 = g_2_6.copy()\n",
    "rem = [e for e in g_top_2_6.nodes() if e not in s ]\n",
    "\n",
    "for e in rem:\n",
    "    g_top_2_6.remove_node(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,8))\n",
    "plt.grid(False)\n",
    "nx.draw_networkx_nodes(g_top_2_6, locations, node_size=40, node_color='black')\n",
    "nx.draw_networkx_edges(g_top_2_6, locations, alpha=0.8, width=0.2, edge_color='gray')\n",
    "plt.axis(False)\n",
    "\n",
    "mplleaflet.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Important nodes last week\n",
    "s = set(np.argsort(f_23_27)[-200:])\n",
    "g_top_23_27 = g_23_27.copy()\n",
    "rem = [e for e in g_top_23_27.nodes() if e not in s ]\n",
    "\n",
    "for e in rem:\n",
    "    g_top_23_27.remove_node(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,8))\n",
    "plt.grid(False)\n",
    "nx.draw_networkx_nodes(g_top_23_27, locations, node_size=40, node_color='black')\n",
    "nx.draw_networkx_edges(g_top_23_27, locations, alpha=0.8, width=0.2, edge_color='gray')\n",
    "plt.axis(False)\n",
    "\n",
    "mplleaflet.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heat Diffusion #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pygsp.filters import Heat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# good nodes , 40, 600, 320"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_nodes = [190]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t[time]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0, 3, 10)\n",
    "\n",
    "h = Heat(G_2_6, tau=t)\n",
    "\n",
    "source = np.zeros(G_2_6.U.shape[0])\n",
    "source[transform_source(source_nodes, isolated_nodes_2_6)] = 75\n",
    "\n",
    "y_2_6 = h.filter(source)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in sorted(isolated_nodes_2_6):\n",
    "    y_2_6 = np.insert(y_2_6, e, 0, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "h = Heat(G_23_27, tau=t)\n",
    "\n",
    "source = np.zeros(G_23_27.U.shape[0])\n",
    "source[transform_source(source_nodes, isolated_nodes_23_27)] = 75\n",
    "\n",
    "y_23_27 = h.filter(source)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in sorted(isolated_nodes_23_27):\n",
    "    y_23_27 = np.insert(y_23_27, e, 0, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hotmap = plt.cm.get_cmap('viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "node_size_2 = np.ones(g_2_6_tot.number_of_nodes()) * 40\n",
    "node_size_2[source_nodes[0]] = 60\n",
    "node_size_23 = np.ones(g_23_27_tot.number_of_nodes()) * 40\n",
    "node_size_23[source_nodes[0]] = 60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "node_size_2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = 9\n",
    "node = source_nodes[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vec_2_6 = y_2_6[:,time].copy()\n",
    "vec_2_6[vec_2_6 > 1] = 1\n",
    "colors_2_6 = [matplotlib.colors.to_hex(e, keep_alpha=False) for e in hotmap(vec_2_6)]\n",
    "\n",
    "vec_2_6_6 = y_2_6[:,9].copy()\n",
    "vec_2_6_6[vec_2_6_6 > 1] = 1\n",
    "colors_2_6_6 = [matplotlib.colors.to_hex(e, keep_alpha=False) for e in hotmap(vec_2_6_6)]\n",
    "\n",
    "vec_23_27 = y_23_27[:, time].copy()\n",
    "vec_23_27[vec_23_27 > 1] = 1\n",
    "colors_23_27 = [matplotlib.colors.to_hex(e, keep_alpha=False) for e in hotmap(vec_23_27)]\n",
    "\n",
    "vec_23_27_6 = y_23_27[:,9].copy()\n",
    "vec_23_27_6[vec_23_27_6 > 1] = 1\n",
    "colors_23_27_6 = [matplotlib.colors.to_hex(e, keep_alpha=False) for e in hotmap(vec_23_27_6)]\n",
    "\n",
    "colors_2_6[source_nodes[0]] = matplotlib.colors.to_hex('red', keep_alpha=False)\n",
    "colors_23_27[source_nodes[0]] = matplotlib.colors.to_hex('red', keep_alpha=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_nodes_2 = [x for x in g_2_6_tot.nodes() if x not in out_nodes]\n",
    "draw_edges_2 = [x for x in g_2_6_tot.edges() if x[0] not in out_nodes and x[1] not in out_nodes]\n",
    "\n",
    "draw_nodes_23 = [x for x in g_23_27_tot.nodes() if x not in out_nodes]\n",
    "draw_edges_23 = [x for x in g_23_27_tot.edges() if x[0] not in out_nodes and x[1] not in out_nodes]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in reversed(list(out_nodes)):\n",
    "    node_size_2 = np.delete(node_size_2, x)\n",
    "    colors_2_6 = np.delete(colors_2_6, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in reversed(list(out_nodes)):\n",
    "    node_size_23 = np.delete(node_size_23, x)\n",
    "    colors_23_27 = np.delete(colors_23_27, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,15))\n",
    "mat = np.arange(0, 1, 0.1)\n",
    "label = np.array([round(x, 2) for x in mat])\n",
    "mat = np.expand_dims(mat, axis=0)\n",
    "plt.imshow(mat, origin=\"lower\", cmap=cmap, interpolation='nearest')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9], label[[0, 3, 6, 9]], fontsize=22)\n",
    "plt.set_cmap(cmap=hotmap)\n",
    "\n",
    "plt.savefig('plots/heat_cmap.jpg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,8))\n",
    "plt.subplot(1,2,1)\n",
    "\n",
    "plt.title('First Week', fontsize=22)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "nx.draw_networkx_nodes(g_2_6_tot, locations, nodelist=draw_nodes_2,  node_size=node_size_2, node_color=colors_2_6, cmap=cmap, vmin=0, vmax=1)\n",
    "nx.draw_networkx_edges(g_2_6_tot, locations, edgelist=draw_edges_2, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.title('Fourth Week', fontsize=22)\n",
    "plt.grid(False)\n",
    "plt.axis(False)\n",
    "nx.draw_networkx_nodes(g_23_27_tot, locations, nodelist=draw_nodes_23, node_size=node_size_23, node_color=colors_23_27, cmap=cmap, vmin=0, vmax=1)\n",
    "nx.draw_networkx_edges(g_23_27_tot, locations, edgelist=draw_edges_23, alpha=0.8, width=0.2, edge_color='black')\n",
    "\n",
    "plt.savefig('plots/heat_{}.jpeg'.format(source_nodes[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(np.arange(len(degs_23_27)), degs_23_27[sort_deg])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "degs_2_6[sort_deg][700:800]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "degs_23_27[sort_deg][700:800]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "degs_23_27[sort_deg][702]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sort_deg[702]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "id_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "25/69"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "degs_2_6[895]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "degs_23_27[895]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[df['start station id'] == 3096]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx.algorithms.cluster.clustering(g_2_6, 320)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx.algorithms.cluster.clustering(g_23_27, 320)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logi_id['5335.03']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 35, 24, 874\n",
    "# 320 - 3096 Mccarren park\n",
    "# 449 - 3318 Marx Brothers Playground\n",
    "# 479 - 3349 Brooklyn Botanic Garden \n",
    "# 352 - 3137 Central Park Manhattan\n",
    "# 408 - 3256 Pier 40 - Hudson River Park\n",
    "# 895 - 24 Beooklyn Bridge\n",
    "# 190 - 115 WB Bridge\n",
    "source_nodes[0] = 190"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_2_6 = np.mean(f_2_6)\n",
    "avg_23_27 = np.mean(f_23_27)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_2_6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_23_27"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.max(cols))\n",
    "print(np.min(cols))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vmin = 0\n",
    "vmax = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.sum(nx.adj_matrix(g_2_6_tot)[source_nodes[0],:].todense()))\n",
    "\n",
    "print(np.sum(nx.adj_matrix(g_23_27_tot)[source_nodes[0],:].todense()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_2_6 = dict(zip(list(g_2_6.nodes), np.arange(len(g_2_6))))\n",
    "map_23_27 = dict(zip(list(g_23_27.nodes), np.arange(len(g_23_27))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "two_order_ns = set()\n",
    "n = source_nodes[0]\n",
    "two_order_ns.add(n)\n",
    "for e in nx.neighbors(g_2_6_tot, n):\n",
    "    two_order_ns.add(e)\n",
    "\"\"\"new_nodes = set()\n",
    "for e in two_order_ns:\n",
    "    for n in nx.neighbors(g_2_6_tot, e):\n",
    "        new_nodes.add(n)\n",
    "        \n",
    "two_order_ns = two_order_ns.union(new_nodes)\"\"\"\n",
    "    \n",
    "g = nx.subgraph(g_2_6, list(two_order_ns))\n",
    "\n",
    "cols=[]\n",
    "for e in g.nodes:\n",
    "    if e==source_nodes[0]:\n",
    "        #cols.append('red')\n",
    "        cols.append(0)\n",
    "    else:\n",
    "        # cols.append('black')\n",
    "        cols.append(f_2_6[map_2_6[e]] / avg_2_6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "two_order_ns = set()\n",
    "n = source_nodes[0]\n",
    "two_order_ns.add(n)\n",
    "for e in nx.neighbors(g_23_27_tot, n):\n",
    "    two_order_ns.add(e)\n",
    "new_nodes = set()\n",
    "\"\"\"for e in two_order_ns:\n",
    "    for n in nx.neighbors(g_23_27_tot, e):\n",
    "        new_nodes.add(n)\"\"\"\n",
    "\"\"\"        \n",
    "two_order_ns = two_order_ns.union(new_nodes)\"\"\"\n",
    "    \n",
    "g = nx.subgraph(g_23_27, list(two_order_ns))\n",
    "\n",
    "cols=[]\n",
    "for e in g.nodes:\n",
    "    if e==source_nodes[0]:\n",
    "        #cols.append('red')\n",
    "        cols.append(0)\n",
    "    else:\n",
    "        # cols.append('black')\n",
    "        cols.append(f_23_27[map_23_27[e]] / avg_23_27)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newcmap = matplotlib.cm.get_cmap('viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(1,2,1)\n",
    "plt.grid(False)\n",
    "nx.draw_networkx_nodes(g, locations, node_size=85, node_color=cols, cmap=newcmap, vmin=vmin, vmax=vmax)\n",
    "nx.draw_networkx_edges(g, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "nx.draw_networkx_nodes(g, locations, nodelist=[source_nodes[0]], node_size=85, node_color='red', cmap=newcmap)\n",
    "\n",
    "plt.axis(False)\n",
    "\n",
    "mplleaflet.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[df['start station id'] == 3398]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,10))\n",
    "mat = np.arange(0, 1, 0.075)\n",
    "mat = np.expand_dims(mat, 0)\n",
    "plt.imshow(mat, cmap=hotmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 4, 8, 12], [0,0.3,0.6,0.9], fontsize=18)\n",
    "plt.savefig('plots/heat_cmap.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index_23 = dict(zip(g_23_27.nodes(), np.arange(g_23_27.number_of_nodes())))\n",
    "cols_bridge = np.zeros(g_23_27.number_of_nodes())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_man_idx = list(upper_man_idx)\n",
    "queens_idx = list(queens_idx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Illegal Bridges\n",
    "for e in upper_man_idx[:100]:\n",
    "    cols_bridge[index_23[e]] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 12))\n",
    "\n",
    "plt.grid(False)\n",
    "nx.draw_networkx_nodes(g_23_27, locations, node_size=40, node_color=cols_bridge)\n",
    "nx.draw_networkx_edges(g_23_27, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "plt.axis(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index_2_6 = dict(zip(g_2_6.nodes(), np.arange(g_2_6.number_of_nodes())))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zoom = []\n",
    "for k, v in locations.items():\n",
    "    if v[0] > -73.96 and v[0] < -73.90 and v[1] > 40.77 and v[1] < 40.80:\n",
    "        zoom.append(k)\n",
    "zoom = np.array(zoom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    " v[0] > -73.96 and v[0] < -73.90 and v[1] > 40.76 and v[1] < 40.80\n",
    "manhattan = list(np.arange(55))\n",
    "manhattan.extend([60, 69, 96, 99])\n",
    "queens = list(set(np.arange(102)).difference(set(manhattan)))\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "manhattan = list(np.arange(20))\n",
    "# manhattan.extend([60, 69, 96, 99])\n",
    "queens = list(set(np.arange(40)).difference(set(manhattan)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols_bridge = ['blue' for e in np.arange(g_2_6.number_of_nodes())]\n",
    "\n",
    "for e in zoom[manhattan]:\n",
    "    if index_2_6.get(e) is not None:\n",
    "        cols_bridge[index_2_6[e]] = 'yellow'\n",
    "\n",
    "for e in zoom[queens]:\n",
    "    if index_2_6.get(e) is not None:\n",
    "        cols_bridge[index_2_6[e]] = 'red'\n",
    "        \n",
    "plt.figure(figsize=(8, 12))\n",
    "\n",
    "plt.grid(False)\n",
    "nx.draw_networkx_nodes(g_2_6, locations, node_size=40, node_color=cols_bridge, cmap=colors_bridge)\n",
    "nx.draw_networkx_edges(g_2_6, locations, alpha=0.8, width=0.2, edge_color='black')\n",
    "plt.axis(False) \n",
    "\n",
    "# mplleaflet.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first_week\n",
    "\n",
    "# start station id --> 3 \n",
    "# end station id --> 7\n",
    "set_manhattan = set(zoom[manhattan])\n",
    "set_queens = set(zoom[queens])\n",
    "\n",
    "start_row = d[2] \n",
    "end_row = d[7]\n",
    "\n",
    "first_cross = []\n",
    "\n",
    "for row in df.iloc[start_row:end_row].iterrows():\n",
    "    if id_index[row[1][3]] in set_manhattan and id_index[row[1][7]] in set_queens:\n",
    "        first_cross.append(row[0])\n",
    "    elif id_index[row[1][7]] in set_manhattan and id_index[row[1][3]] in set_queens:\n",
    "        first_cross.append(row[0])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# last\n",
    "\n",
    "# start station id --> 3 \n",
    "# end station id --> 7\n",
    "set_manhattan = set(zoom[manhattan])\n",
    "set_queens = set(zoom[queens])\n",
    "\n",
    "start_row = d[23] \n",
    "end_row = d[28]\n",
    "\n",
    "last_cross = []\n",
    "\n",
    "for row in df.iloc[start_row:end_row].iterrows():\n",
    "    if id_index[row[1][3]] in set_manhattan and id_index[row[1][7]] in set_queens:\n",
    "        last_cross.append(row[0])\n",
    "    elif id_index[row[1][7]] in set_manhattan and id_index[row[1][3]] in set_queens:\n",
    "        last_cross.append(row[0])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.iloc[first_cross]['tripduration'].describe()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = np.array(df.iloc[first_cross]['tripduration'])\n",
    "f = f[f<2400]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = np.array(df.iloc[last_cross]['tripduration'])\n",
    "l = l[l<2400]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = plt.hist(f, color='blue')\n",
    "plt.grid(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = plt.hist(l, color='blue')\n",
    "plt.grid(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.iloc[last_cross]['start station id']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Metro Data #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import shapefile as shp\n",
    "import seaborn as sns\n",
    "from pyproj import Proj, transform\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "\n",
    "from shapely.geometry import shape, Point\n",
    "\n",
    "import shapefile\n",
    "\n",
    "from utils.shapefile import read_shapefile, plot_shape, plot_map, plot_map_fill"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# UTILS\n",
    "\n",
    "def find_day_metro(day, month, df):\n",
    "    day = str_int(day)\n",
    "    month = str_int(month)\n",
    "    ids = []\n",
    "    for index, row in df.iterrows():\n",
    "        if row['date'].split('-')[-1] == day and row['date'].split('-')[-2] == month:\n",
    "            ids.append(index)\n",
    "    return ids\n",
    "        \n",
    "def str_int(i):\n",
    "    if i < 10:\n",
    "        return '0' + str(i)\n",
    "    else:\n",
    "        return str(i)\n",
    "\n",
    "def calc_color(data, num_bins=None, ref_bins=None, name='YlOrBr'):\n",
    "    \n",
    "    cmap=plt.get_cmap(name)\n",
    "    \n",
    "    if ref_bins is None:\n",
    "        new_data, bins = pd.qcut(data, num_bins, retbins=True, labels=list(range(num_bins)))\n",
    "    else:\n",
    "        bins = ref_bins\n",
    "        num_bins = len(ref_bins)-1\n",
    "        new_data = pd.cut(data, ref_bins, include_lowest=True, labels=list(range(num_bins)))\n",
    "    m = num_bins - 1\n",
    "    color_ton = []\n",
    "    for val in new_data:\n",
    "        color_ton.append(matplotlib.colors.to_hex(cmap(val/m), keep_alpha=False)) \n",
    "    \n",
    "    #plot color map\n",
    "    \n",
    "    \"\"\"cat = np.arange(num_bins)\n",
    "    ex = cat / (num_bins-1)\n",
    "    ex = np.expand_dims(ex, 1)\n",
    "    \n",
    "    plt.imshow(ex, cmap=cmap)\n",
    "    plt.grid(False)\n",
    "    plt.yticks([0], [''])\n",
    "    _ = plt.axis('on')\"\"\"\n",
    "\n",
    "    return color_ton, bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(style=\"whitegrid\", palette=\"pastel\", color_codes=True)\n",
    "sns.mpl.rc(\"figure\", figsize=(12,8))\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "sf = shapefile.Reader(\"./data/shapefiles_/ZIP_CODE_040114.dbf\")\n",
    "\n",
    "df_shape = read_shapefile(sf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metro_data = pd.read_csv('./data/metro/metro_whole_station.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metro_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find Geo-Coordinates for stations\n",
    "\n",
    "stat_geo = {}\n",
    "for index, row in metro_data.iterrows():\n",
    "    coor = np.array([row['gtfs_longitude'], row['gtfs_latitude']])\n",
    "    stat_geo[row['stop_name']] = coor\n",
    "    \n",
    "# Find xlim and ylim \n",
    "x_lim = [100, -100]\n",
    "y_lim = [100, 0]\n",
    "for k,v in locations.items():\n",
    "    \n",
    "    if v[0] < x_lim[0]:\n",
    "        x_lim[0] = v[0]\n",
    "    elif v[0] > x_lim[1]:\n",
    "        x_lim[1] = v[0]\n",
    "        \n",
    "    if v[1] < y_lim[0]:\n",
    "        y_lim[0] = v[1]\n",
    "    elif v[1] > y_lim[1]:\n",
    "        y_lim[1] = v[1]\n",
    "        \n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find Zip Code for Stations\n",
    "\n",
    "points = [(k, Point(v)) for k,v in stat_geo.items()]\n",
    "multis = [(e.record[0], shape(e.shape.__geo_interface__)) for e in sf.shapeRecords()]\n",
    "\n",
    "stat_zip = {}\n",
    "\n",
    "for stat, pt in points:\n",
    "    point = shape(pt)\n",
    "    for zip_code, poly in multis:\n",
    "        if point.within(poly): \n",
    "            stat_zip[stat] = zip_code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find Zip Code for BS - Stations\n",
    "\n",
    "points_bs = [(k, Point(v)) for k,v in locations.items()]\n",
    "multis_bs = [(e.record[0], shape(e.shape.__geo_interface__)) for e in sf.shapeRecords()]\n",
    "\n",
    "bs_zip = {}\n",
    "not_found = []\n",
    "for stat, pt in points_bs:\n",
    "    point = shape(pt)\n",
    "    for zip_code, poly in multis_bs:\n",
    "        if point.within(poly): \n",
    "            bs_zip[stat] = zip_code\n",
    "    if bs_zip.get(stat) is None:\n",
    "        not_found.append(stat)\n",
    "        \n",
    "# Find f_2_6_tot with all nodes \n",
    "\n",
    "f_24_tot = list(np.sum(adj_24_tot, axis=0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"for n in reversed(not_found):\n",
    "    del f_2_6_tot[n]\"\"\"\n",
    "\n",
    "f_24_tot = np.array(f_24_tot)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zip_index = dict(zip(df_shape['ZIPCODE'], np.arange(len(df_shape))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create indexes for March \n",
    "\n",
    "\"\"\"march_ids = []\n",
    "for d in range(1, 32):\n",
    "    march_ids.append(find_day_metro(d, 3, metro_data))\n",
    "\n",
    "with open('./data/python_variables/march_ids.pickle', 'wb') as file:\n",
    "    pickle.dump(march_ids, file)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('./data/python_variables/march_ids.pickle', 'rb') as file:\n",
    "    march_ids = pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find Daily runs\n",
    "\n",
    "days = [4, 11, 18, 24]\n",
    "\n",
    "days_metro = dict()\n",
    "\n",
    "for day in days:\n",
    "    days_metro[day] = {}\n",
    "    for index, row in metro_data.iloc[march_ids[day-1]].iterrows():\n",
    "        zip_code = stat_zip.get(row['stop_name'])\n",
    "        if zip_code is not None:\n",
    "            if days_metro[day].get(zip_code) is None:\n",
    "                days_metro[day][zip_code] = row['entries'] \n",
    "            else:\n",
    "                days_metro[day][zip_code] += row['entries'] \n",
    "    print(\"Aggregated day {}\".format(day))\n",
    "\n",
    "# Collect By indexes\n",
    "\n",
    "idxs = {}\n",
    "for day in days:\n",
    "    idxs[day] = []\n",
    "    idxs[day].append(np.array([zip_index[str(e)] for e in days_metro[day].keys() ]))\n",
    "    idxs[day].append(np.array(list(days_metro[day].values())))\n",
    "    idxs[day].append(idxs[day][1] / idxs[4][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, v in idxs.items():\n",
    "    v[1] = np.where(v[1] > 150000, 150000, v[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref, bins = pd.qcut(idxs[4][1], 12, retbins=True, labels=list(range(12)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_lab = np.array([int(e) for e in bins])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_lab = np.insert(bin_lab, 0, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For relative\n",
    "\n",
    "bins = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]\n",
    "cmap = plt.get_cmap('jet')\n",
    "plt.figure(figsize=(12,10))\n",
    "mat = bins.copy()\n",
    "mat = np.expand_dims(mat, 0)\n",
    "plt.imshow(mat, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 2, 5, 8, 10], (mat[0,[0, 2, 5, 8, 10]]), fontsize=18)\n",
    "# plt.savefig('plots/metro_relative_cmap.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = plt.get_cmap('jet')\n",
    "plt.figure(figsize=(12,10))\n",
    "mat = bin_lab.copy()\n",
    "mat = np.expand_dims(mat, 0)\n",
    "plt.imshow(mat, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9, 13 ], (mat[0,[0, 3, 6, 9, 13 ]]), fontsize=18)\n",
    "plt.savefig('plots/metro_relative_cmap.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, v in idxs.items():\n",
    "\n",
    "    color_ton, _ = calc_color(v[2], ref_bins=bins, name='jet')\n",
    "\n",
    "    ax = plot_map_fill(v[0], sf, color_ton, x_lim=x_lim, y_lim=y_lim, figsize=(8,6))\n",
    "    \n",
    "    _ = plt.xlabel(\"Longitude\", fontsize=22)\n",
    "    _ = plt.ylabel(\"Latidute\", fontsize=22)\n",
    "    \n",
    "    _ = plt.savefig('plots/metro{}.eps'.format(k))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(f_24_tot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(bs_zip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find BS runs\n",
    "bs_agg = dict()\n",
    "\n",
    "for k, v in enumerate(f_24):\n",
    "    zip_code = bs_zip.get(k)\n",
    "    if zip_code is not None:\n",
    "        if bs_agg.get(zip_code) is None:\n",
    "            bs_agg[zip_code] = v\n",
    "        else:\n",
    "            bs_agg[zip_code] += v \n",
    "print(\"Aggregated week one\")\n",
    "\n",
    "# Collect By indexes\n",
    "\n",
    "bs_tot = []\n",
    "\n",
    "bs_tot.append(np.array([zip_index[str(e)] for e in bs_agg.keys() ]))\n",
    "bs_tot.append(np.array(list(bs_agg.values())))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot BS Aggregated per ZIP CODE\n",
    "\n",
    "# bs_tot[1] = np.where(bs_tot[1] > 6000, 6000, bs_tot[1] )\n",
    "\n",
    "ref, bins = pd.qcut(bs_tot[1], 12, retbins=True, labels=list(range(12)))\n",
    "bin_lab = np.array([int(e) for e in bins])\n",
    "\n",
    "cmap = plt.get_cmap('jet')\n",
    "plt.figure(figsize=(12,10))\n",
    "mat = bins.copy()\n",
    "mat = np.expand_dims(mat, 0)\n",
    "plt.imshow(mat, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9, 12], (bin_lab[[0, 3, 6, 9, 12]]), fontsize=18)\n",
    "\n",
    "plt.savefig('plots/zip_bs_cmap_4.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# bs_tot[1] = np.where(bs_tot[1] > 2000, 2000, bs_tot[1])\n",
    "\n",
    "color_ton, _ = calc_color(bs_tot[1], ref_bins=bins, name='jet')\n",
    "\n",
    "ax = plot_map_fill(bs_tot[0], sf, color_ton, x_lim=x_lim, y_lim=y_lim, figsize=(8,6))\n",
    "\n",
    "_ = plt.xlabel(\"Longitude\", fontsize=22)\n",
    "_ = plt.ylabel(\"Latidute\", fontsize=22)\n",
    "\n",
    "_ = plt.savefig('plots/bs_zip_24.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Left\n",
    "\n",
    "zip_left = pd.read_csv('./data/zip_left.csv', names=['zipcode', 'perc'], header=None)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zip_left.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_left = [zip_index[str(e)] for e in zip_left['zipcode']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = [0, 1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_ton, _ = calc_color(list(zip_left['perc']), ref_bins=bins)\n",
    "\n",
    "ax = plot_map_fill(idx_left, sf, color_ton, x_lim=x_lim, y_lim=y_lim, figsize=(8,6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"import geopandas as gpd\n",
    "\n",
    "tmp = gpd.GeoDataFrame.from_file('./data/shapefiles/ZIP_CODE_040114.shp')\n",
    "\n",
    "tmpWGS84 = tmp.to_crs({'proj':'longlat', 'ellps':'WGS84', 'datum':'WGS84'})\n",
    "\n",
    "tmpWGS84.to_file('./data/shapefiles_/ZIP_CODE_040114.shp')\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Covid cases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_covid = pd.read_csv('data/coronavirus-data/data-by-modzcta.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_covid.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "covid_zip = {}\n",
    "\n",
    "set_zip = set(zip_index.keys())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, row in df_covid.iterrows():\n",
    "    if str(row['MODIFIED_ZCTA']) in set_zip:\n",
    "        covid_zip[zip_index[str(row['MODIFIED_ZCTA'])]] = row['COVID_CASE_COUNT']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_covid['COVID_CASE_COUNT'].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(list(covid_zip.values()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_covid, bins_covid = calc_color(list(covid_zip.values()), num_bins=14, name='jet')\n",
    "ax = plot_map_fill(list(covid_zip.keys()), sf, color_covid, x_lim=[-74.05, -73.90], y_lim=[40.6, 40.80], figsize=(8,6))\n",
    "\n",
    "# plt.savefig('plots/covid_cases.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = plt.get_cmap('jet')\n",
    "plt.figure(figsize=(12,10))\n",
    "bins_covid_int = [int(x) for x in bins_covid]\n",
    "mat = bins_covid_int.copy()\n",
    "mat = np.expand_dims(mat, 0)\n",
    "plt.imshow(mat, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9, 12], (mat[0, [0, 3, 6, 9, 12]]), fontsize=18)\n",
    "\n",
    "plt.savefig('plots/covid_cmap.eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analsys of new network "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio_adj = adj_23_27_tot / (adj_2_6_tot+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio_adj[ratio_adj < 1] = 0\n",
    "\n",
    "for i in range(ratio_adj.shape[0]):\n",
    "    ratio_adj[i,i] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(ratio_adj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find connections that has grown more than 5 times with at least 20 runs \n",
    "idx = np.argwhere(np.logical_and(ratio_adj >= 2, adj_23_27_tot >=  10))\n",
    "idx_remove = np.argwhere(np.logical_or(ratio_adj < 2, adj_23_27_tot < 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in idx_remove:\n",
    "    ratio_adj[e[0], e[1]] = 0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "g_ratio = nx.from_numpy_matrix(ratio_adj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolated_nodes_ratio = []\n",
    "for e in nx.connected_components(g_ratio):\n",
    "    if len(e) == 1:\n",
    "        isolated_nodes_ratio.append(e.pop())\n",
    "\n",
    "g_ratio.remove_nodes_from(isolated_nodes_ratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 12))\n",
    "\n",
    "nx.draw_networkx_nodes(g_ratio, locations, node_size=40, node_color='black')\n",
    "nx.draw_networkx_edges(g_ratio, locations, alpha=0.8, width=0.9, edge_color='black')\n",
    "plt.axis(False)\n",
    "plt.grid(False)\n",
    "\n",
    "mplleaflet.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decreasing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio_adj_dec = adj_23_27_tot / (adj_2_6_tot+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio_adj_dec[ratio_adj_dec > 1] = 0\n",
    "\n",
    "for i in range(ratio_adj_dec.shape[0]):\n",
    "    ratio_adj_dec[i,i] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find connections that has grown more than 5 times with at least 20 runs \n",
    "idx = np.argwhere(np.logical_and(ratio_adj_dec <= 0.5, adj_23_27_tot >=  10))\n",
    "idx_remove = np.argwhere(np.logical_or(ratio_adj_dec > 0.5, adj_23_27_tot < 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for e in idx_remove:\n",
    "    ratio_adj_dec[e[0], e[1]] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g_dec = nx.from_numpy_matrix(ratio_adj_dec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "isolated_nodes_ratio_dec = []\n",
    "for e in nx.connected_components(g_dec):\n",
    "    if len(e) == 1:\n",
    "        isolated_nodes_ratio_dec.append(e.pop())\n",
    "\n",
    "g_dec.remove_nodes_from(isolated_nodes_ratio_dec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 12))\n",
    "\n",
    "nx.draw_networkx_nodes(g_ratio, locations, node_size=20, node_color='black')\n",
    "nx.draw_networkx_edges(g_ratio, locations, alpha=0.8, width=0.9, edge_color='black')\n",
    "plt.axis(False)\n",
    "plt.grid(False)\n",
    "\n",
    "\n",
    "nx.draw_networkx_nodes(g_dec, locations, node_size=20, node_color='black')\n",
    "nx.draw_networkx_edges(g_dec, locations, alpha=0.8, width=0.9, edge_color='red')\n",
    "plt.axis(False)\n",
    "plt.grid(False)\n",
    "\n",
    "mplleaflet.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G_2_6.e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 10))\n",
    "plt.plot(G_2_6.e, color='blue', linewidth=2)\n",
    "_ = plt.plot(G_9_13.e, color='gray', linewidth=2)\n",
    "_ = plt.plot(G_16_20.e, color='red', linewidth=2)\n",
    "_ = plt.plot(G_23_27.e, color='green', linewidth=2)\n",
    "\n",
    "\n",
    "plt.grid(True)\n",
    "_ = plt.xticks(fontsize=22)\n",
    "_ = plt.yticks(fontsize=22)\n",
    "_ = plt.xlabel(\"Eigenvalue Index\", fontsize=26)\n",
    "\n",
    "plt.legend(['First', 'Second', 'Third', 'Last'], fontsize=22)\n",
    "\n",
    "plt.savefig('plots/eigens_tot.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,10))\n",
    "plt.plot(G_2_6.e[:5], color='blue', linewidth=2)\n",
    "_ = plt.plot(G_9_13.e[:5], color='gray', linewidth=2)\n",
    "_ = plt.plot(G_16_20.e[:5], color='red', linewidth=2)\n",
    "_ = plt.plot(G_23_27.e[:5], color='green', linewidth=2)\n",
    "\n",
    "plt.grid(True)\n",
    "_ = plt.xticks(fontsize=22)\n",
    "_ = plt.yticks(fontsize=22)\n",
    "_ = plt.xlabel(\"Eigenvalue Index\", fontsize=26)\n",
    "\n",
    "plt.legend(['First', 'Second', 'Third', 'Last'], fontsize=22)\n",
    "# plt.savefig('plots/eigens_100.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G_2_6.compute_laplacian(lap_type='normalized')\n",
    "G_9_13.compute_laplacian(lap_type='normalized')\n",
    "G_16_20.compute_laplacian(lap_type='normalized')\n",
    "G_23_27.compute_laplacian(lap_type='normalized')\n",
    "\n",
    "G_2_6.compute_fourier_basis(recompute=True)\n",
    "G_9_13.compute_fourier_basis(recompute=True)\n",
    "G_16_20.compute_fourier_basis(recompute=True)\n",
    "G_23_27.compute_fourier_basis(recompute=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,10))\n",
    "plt.plot(G_2_6.e[:5], color='blue', linewidth=2)\n",
    "_ = plt.plot(G_9_13.e[:5], color='gray', linewidth=2)\n",
    "_ = plt.plot(G_16_20.e[:5], color='red', linewidth=2)\n",
    "_ = plt.plot(G_23_27.e[:5], color='green', linewidth=2)\n",
    "\n",
    "plt.grid(True)\n",
    "_ = plt.xticks(fontsize=22)\n",
    "_ = plt.yticks(fontsize=22)\n",
    "_ = plt.xlabel(\"Eigenvalue Index\", fontsize=26)\n",
    "\n",
    "plt.legend(['First', 'Second', 'Third', 'Last'], fontsize=22)\n",
    "# plt.savefig('plots/eigens_100.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.max(G_2_6.e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.ones((5,5)) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a @ a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(list(nx.triangles(g_2_6).values())) / np.sum(adj_2_6 @ adj_2_6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(list(nx.triangles(g_9_13).values())) / np.sum(adj_9_13 @ adj_9_13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(list(nx.triangles(g_16_20).values())) / np.sum(adj_16_20 @ adj_16_20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(list(nx.triangles(g_23_27).values())) / np.sum(adj_23_27 @ adj_23_27)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gt = nx.complete_graph(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_html('data/view-source_https___www.newyork-demographics.com_zip_codes_by_population.html')[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# start at 70 \n",
    "re.split('>|<','<a href=\\\"11385-demographics\">11385</a>')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop_zip = {}\n",
    "i = 5\n",
    "for k,v in data.iloc[70:].iterrows():\n",
    "    i += 1\n",
    "    if type(v[1]) is str and \"-demographics\\\">\" in v[1]:\n",
    "        zip_code = re.split('>|<',v[1])[2]\n",
    "        i = 0\n",
    "    elif i == 3:\n",
    "        pop_zip[zip_code] = int(v[1])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop_zip_id = {}\n",
    "for k, v in pop_zip.items():\n",
    "    if k in set_zip:\n",
    "        pop_zip_id[zip_index[k]] = pop_zip[k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_covid, bins_pop = calc_color(list(pop_zip_id.values()), num_bins=12, name='jet')\n",
    "ax = plot_map_fill(list(pop_zip_id.keys()), sf, color_covid, x_lim=x_lim, y_lim=y_lim, figsize=(8,6))\n",
    "\n",
    "plt.savefig('plots/covid_cases.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = plt.get_cmap('jet')\n",
    "plt.figure(figsize=(12,10))\n",
    "bins_int = [int(x) for x in bins_pop]\n",
    "mat = bins_int.copy()\n",
    "mat = np.expand_dims(mat, 0)\n",
    "plt.imshow(mat, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9, 12], (mat[0, [0, 3, 6, 9, 12]]), fontsize=18)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(covid_zip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(pop_zip_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "perc_zip ={}\n",
    "for k, v in covid_zip.items():\n",
    "    perc_zip[k] = v / int(pop_zip_id[k] / 1000 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_covid, bins = calc_color(list(perc_zip.values()), num_bins=12, name='jet')\n",
    "ax = plot_map_fill(list(perc_zip.keys()), sf, color_covid, x_lim=x_lim, y_lim=y_lim, figsize=(8,6))\n",
    "plt.xlabel(\"Longitude\", fontsize=22)\n",
    "plt.ylabel(\"Latitude\",  fontsize=22)\n",
    "\n",
    "plt.savefig('plots/covid_perc.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = plt.get_cmap('jet')\n",
    "plt.figure(figsize=(12,10))\n",
    "bins_int = [int(x) for x in bins]\n",
    "mat = bins_int.copy()\n",
    "mat = np.expand_dims(mat, 0)\n",
    "plt.imshow(mat, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9, 12], (mat[0, [0, 3, 6, 9, 12]]), fontsize=18)\n",
    "\n",
    "plt.savefig('plots/cmap_covid_perc.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_covid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "death_zip = {}\n",
    "for k, row in df_covid.iterrows():\n",
    "    if str(row['MODIFIED_ZCTA']) in set_zip:\n",
    "        death_zip[zip_index[str(row['MODIFIED_ZCTA'])]] = row['COVID_DEATH_COUNT']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "death_perc_zip ={}\n",
    "for k, v in death_zip.items():\n",
    "    perc_zip[k] = v / int(pop_zip_id[k] / 1000 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = [0, 1, 2, 3, 4, 5, 6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_covid, bins = calc_color(list(perc_zip.values()), ref_bins=bins, name='jet')\n",
    "ax = plot_map_fill(list(perc_zip.keys()), sf, color_covid, x_lim=[-74.05, -73.8], y_lim=[40.6, 40.9], figsize=(8,6))\n",
    "\n",
    "plt.savefig('plots/covid_death_perc.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = plt.get_cmap('jet')\n",
    "plt.figure(figsize=(12,10))\n",
    "bins_int = [int(x) for x in bins]\n",
    "mat = bins_int.copy()\n",
    "mat = np.expand_dims(mat, 0)\n",
    "plt.imshow(mat, cmap=cmap, interpolation='nearest')\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 2, 4, 5], (mat[0, [0, 2, 4, 5]]), fontsize=18)\n",
    "\n",
    "plt.savefig('plots/cmap_death_perc.eps')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# April #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_4 = pd.read_csv('./data/202004-citibike-tripdata.csv')\n",
    "n_start_station = len(df_4['start station id'].unique())\n",
    "n_end_station = len(df_4['end station id'].unique())\n",
    "\n",
    "\n",
    "total_station_id = set(df_4['start station id']).union(set(df_4['end station id']))\n",
    "n_tot_station = len(total_station_id)\n",
    "\n",
    "id_index_4 = dict(zip(sorted(total_station_id), np.arange(n_tot_station)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locations_4 = dict()\n",
    "for e in id_index_4.keys():\n",
    "    if df_4[df_4['start station id'] == e]['start station latitude'].shape[0]:\n",
    "        locations_4[id_index_4[e]] = (df_4[df_4['start station id'] == e]['start station longitude'].iloc[0],\n",
    "                                  df_4[df_4['start station id'] == e]['start station latitude'].iloc[0])\n",
    "    else:\n",
    "        locations_4[id_index_4[e]] = (df_4[df_4['end station id'] == e]['end station longitude'].iloc[0],\n",
    "                                  df_4[df_4['end station id'] == e]['end station latitude'].iloc[0])\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dictionary with days starting indexes\n",
    "\"\"\"d_4 = days_index(df_4)\n",
    "with open('./data/python_variables/april_index.pickle', 'wb') as file:\n",
    "    pickle.dump(d_4, file)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load day indexes for March\n",
    "with open('./data/python_variables/april_index.pickle', 'rb') as file:\n",
    "    d_4 = pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_a = [e for e in d_4.values()]\n",
    "d_a.append(df_4.shape[0])\n",
    "d_a = np.array(d_a)\n",
    "rides_day_4 = d_a[1:] - d_a[:-1]\n",
    "\n",
    "# Rides Per Day\n",
    "plt.figure(figsize=(12,8))\n",
    "_ = plt.plot(np.arange(1, 31), rides_day_4, color='blue', marker='*')\n",
    "# _ = plt.xticks(np.arange(1, 8), [1, 5, 10, 15, 20, 25, 31], fontsize=18)\n",
    "_ = plt.xlabel(\"Day of April\", fontsize=24, color='black')\n",
    "_ = plt.ylabel(\"Rides per day\", fontsize=24, color='black')\n",
    "_ = plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)\n",
    "plt.grid(False)\n",
    "# plt.savefig('./plots/days_aggregated.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_4.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_4['usertype'].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# May # "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_5 = pd.read_csv('./data/202005-citibike-tripdata.csv')\n",
    "n_start_station = len(df_5['start station id'].unique())\n",
    "n_end_station = len(df_5['end station id'].unique())\n",
    "\n",
    "\n",
    "total_station_id = set(df_5['start station id']).union(set(df_5['end station id']))\n",
    "n_tot_station = len(total_station_id)\n",
    "\n",
    "id_index_5 = dict(zip(sorted(total_station_id), np.arange(n_tot_station)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locations_5 = dict()\n",
    "for e in id_index_5.keys():\n",
    "    if df_5[df_5['start station id'] == e]['start station latitude'].shape[0]:\n",
    "        locations_5[id_index_5[e]] = (df_5[df_5['start station id'] == e]['start station longitude'].iloc[0],\n",
    "                                  df_5[df_5['start station id'] == e]['start station latitude'].iloc[0])\n",
    "    else:\n",
    "        locations_5[id_index_5[e]] = (df_5[df_5['end station id'] == e]['end station longitude'].iloc[0],\n",
    "                                  df_5[df_5['end station id'] == e]['end station latitude'].iloc[0])\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dictionary with days starting indexes\n",
    "\"\"\"d_5 = days_index(df_5)\n",
    "with open('./data/python_variables/may_index.pickle', 'wb') as file:\n",
    "    pickle.dump(d_5, file)\"\"\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load day indexes for March\n",
    "with open('./data/python_variables/april_index.pickle', 'rb') as file:\n",
    "    d_4 = pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_a = [e for e in d_5.values()]\n",
    "d_a.append(df_5.shape[0])\n",
    "d_a = np.array(d_a)\n",
    "rides_day_5 = d_a[1:] - d_a[:-1]\n",
    "\n",
    "# Rides Per Day\n",
    "plt.figure(figsize=(12,8))\n",
    "_ = plt.plot(np.arange(1, 32), rides_day_5, color='blue', marker='*')\n",
    "# _ = plt.xticks(np.arange(1, 8), [1, 5, 10, 15, 20, 25, 31], fontsize=18)\n",
    "_ = plt.xlabel(\"Day of April\", fontsize=24, color='black')\n",
    "_ = plt.ylabel(\"Rides per day\", fontsize=24, color='black')\n",
    "_ = plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)\n",
    "plt.grid(False)\n",
    "# plt.savefig('./plots/days_aggregated.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(rides_day_5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_5.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# June # "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_6 = pd.read_csv('./data/CitiBike_Covid_Padua.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = [ int(e) for e in df_6.logicalTerminal.unique()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(df_6['logicalTerminal'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_6.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Workforce #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_wf = pd.read_csv('data/CitiBike_2020_Padua.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_wf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_log_coor = pd.read_csv('data/logical_coordinates.csv')\n",
    "df_log_coor.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logi_id = dict()\n",
    "for i, e in df_log_coor.iterrows():\n",
    "    for k,v in locations.items():\n",
    "        if np.round(e.latitude, 4) == np.round(v[1], 4) and np.round(e.longitude, 4) == np.round(v[0], 4):\n",
    "            logi_id[e.logicalTerminal] = k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_wf.iloc[0]['start_time'].split('/')[1][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "id_wf = []\n",
    "for i,e in df_wf.iterrows():\n",
    "    if e.CW == True:\n",
    "        id_wf.append(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "id_not_wf = []\n",
    "for i,e in df_wf.iterrows():\n",
    "    if e.CW == False:\n",
    "        id_not_wf.append(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wf_month = dict()\n",
    "for e in id_wf:\n",
    "    m = df_wf.iloc[e]['start_time'].split('/')[1][1]\n",
    "    if wf_month.get(m) is not None:\n",
    "        wf_month[m].append(e)\n",
    "    else:\n",
    "        wf_month[m] = [e]\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "not_wf_month = dict()\n",
    "for e in id_not_wf:\n",
    "    m = df_wf.iloc[e]['start_time'].split('/')[1][1]\n",
    "    if not_wf_month.get(m) is not None:\n",
    "        not_wf_month[m].append(e)\n",
    "    else:\n",
    "        not_wf_month[m] = [e]\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(wf_month['3']) / len(not_wf_month['3'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(wf_month['4']) / (len(not_wf_month['4']) + len(wf_month['4']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(wf_month['5']) / (len(not_wf_month['5']) + len(wf_month['5']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(wf_month['6']) / (len(not_wf_month['6']) + len(wf_month['6']) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = 0\n",
    "for e in wf_month:\n",
    "    s += len(wf_month[e])\n",
    "print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = 0\n",
    "for e in not_wf_month:\n",
    "    s += len(not_wf_month[e])\n",
    "print(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# May WF #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_wf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First Week\n",
    "days = {4, 5, 6, 7, 8}\n",
    "\n",
    "\n",
    "# Find all connections with weights\n",
    "adj_4_8_tot = np.zeros((n_tot_station, n_tot_station))\n",
    "\n",
    "for index, row in df_wf.iloc[wf_month['5']].iterrows():\n",
    "    if int(row.start_time.split('/')[2].split()[0]) in days:\n",
    "        s_s = str(row['logicalTerminal'])\n",
    "        e_s = str(row['logicalTerminal.1'])\n",
    "        if logi_id.get(s_s) is not None and logi_id.get(e_s) is not None:\n",
    "            adj_4_8_tot[logi_id[s_s], logi_id[e_s]] += 1\n",
    "            adj_4_8_tot[logi_id[e_s], logi_id[s_s]] += 1\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adj_4_8 = adj_4_8_tot.copy()\n",
    "adj_4_8[adj_4_8 <= ths] = 0\n",
    "\"\"\"adj_23_27[adj_23_27 >= ths] = 1\"\"\"\n",
    "print(\"Weighted number of rides : {}\".format(np.sum(adj_4_8) / 2))\n",
    "\n",
    "# Remove Self-loops\n",
    "for i in range(n_tot_station):\n",
    "    adj_4_8[i, i] = 0\n",
    "    \n",
    "g_4_8 = nx.from_numpy_matrix(adj_4_8)\n",
    "\n",
    "g_4_8.remove_nodes_from(out_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_4_8 = np.sum(adj_4_8, axis=0)\n",
    "f_4_8 = np.delete(f_4_8, list(out_nodes))\n",
    "f_4_8 = f_4_8 / np.mean(f_4_8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(f_4_8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = matplotlib.cm.get_cmap('jet')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14,15))\n",
    "mat = np.arange(0, 300, 30)\n",
    "label = np.array([round(x, 2) for x in mat])\n",
    "mat = np.expand_dims(mat, axis=0)\n",
    "plt.imshow(mat, origin=\"lower\", cmap=cmap, interpolation='nearest')\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.yticks([0], [''])\n",
    "_ = plt.xticks([0, 3, 6, 9], label[[0, 3, 6, 9]], fontsize=22)\n",
    "plt.set_cmap(cmap=cmap)\n",
    "\n",
    "plt.savefig('plots/means_cmap_wf_may.jpeg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hottest = np.argsort(f_4_8)[-10:]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_4_8 = dict(zip(np.arange(len(locations)), list(g_4_8.nodes)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 12))\n",
    "ax = plt.gca()\n",
    "\n",
    "\n",
    "nx.draw_networkx_nodes(g_4_8, locations,  node_size=65, node_color=f_4_8, cmap=cmap, vmin=0, vmax=5)\n",
    "nx.draw_networkx_edges(g_4_8, locations, alpha=0.8, width=0.2, edge_color='black', ax=ax)\n",
    "\n",
    "plt.grid(False)\n",
    "_ = plt.axis('off')\n",
    "\n",
    "circle = plt.Circle(locations[map_4_8[hottest[-1]]], 0.008, color='red', fill=False, zorder=10)\n",
    "ax.set_aspect(1)\n",
    "ax.add_patch(circle)\n",
    "plt.text(locations[map_4_8[hottest[-1]]][0] - 0.015, locations[map_4_8[hottest[-1]]][1] + 0.01,\n",
    "         'Presbyterian/Cornell', fontsize=15, zorder=100, color='green', weight='bold')\n",
    "\n",
    "circle = plt.Circle(locations[map_4_8[hottest[-2]]], 0.008, color='red', fill=False, zorder=10)\n",
    "ax.set_aspect(1)\n",
    "ax.add_patch(circle)\n",
    "plt.text(locations[map_4_8[hottest[-2]]][0] - 0.015, locations[map_4_8[hottest[-2]]][1] - 0.01,\n",
    "         'Bellevue/Tisch', fontsize=15, zorder=100, color='green', weight='bold')\n",
    "\n",
    "circle = plt.Circle(locations[map_4_8[hottest[-6]]], 0.008, color='red', fill=False, zorder=10)\n",
    "ax.set_aspect(1)\n",
    "ax.add_patch(circle)\n",
    "plt.text(locations[map_4_8[hottest[-6]]][0] - 0.015, locations[map_4_8[hottest[-6]]][1] + 0.01,\n",
    "         'Mt.Sinai', fontsize=15, zorder=100, color='green', weight='bold')\n",
    "\n",
    "circle = plt.Circle(locations[map_4_8[hottest[-7]]], 0.008, color='red', fill=False, zorder=10)\n",
    "ax.set_aspect(1)\n",
    "ax.add_patch(circle)\n",
    "plt.text(locations[map_4_8[hottest[-7]]][0] - 0.035, locations[map_4_8[hottest[-7]]][1] + 0.01,\n",
    "         'Mt.Sinai West', fontsize=15, zorder=100, color='green', weight='bold')\n",
    "\n",
    "circle = plt.Circle(locations[map_4_8[hottest[-9]]], 0.008, color='red', fill=False, zorder=10)\n",
    "ax.set_aspect(1)\n",
    "ax.add_patch(circle)\n",
    "plt.text(locations[map_4_8[hottest[-9]]][0] - 0.03, locations[map_4_8[hottest[-9]]][1] - 0.01,\n",
    "         'Presbyterian Lower Manhattan', fontsize=15, zorder=100, color='green', weight='bold')\n",
    "\n",
    "plt.show()\n",
    "fig.savefig('plots/workforce_may_first_week.jpeg')\n",
    "\n",
    "print(ax.dataLim)\n",
    "\n",
    "#mplleaflet.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First Week\n",
    "days = {4, 5, 6, 7, 8}\n",
    "\n",
    "\n",
    "# Find all connections with weights\n",
    "not_adj_4_8_tot = np.zeros((n_tot_station, n_tot_station))\n",
    "\n",
    "for index, row in df_wf.iloc[not_wf_month['5']].iterrows():\n",
    "    if int(row.start_time.split('/')[2].split()[0]) in days:\n",
    "        s_s = str(row['logicalTerminal'])\n",
    "        e_s = str(row['logicalTerminal.1'])\n",
    "        if logi_id.get(s_s) is not None and logi_id.get(e_s) is not None:\n",
    "            not_adj_4_8_tot[logi_id[s_s], logi_id[e_s]] += 1\n",
    "            not_adj_4_8_tot[logi_id[e_s], logi_id[s_s]] += 1\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "not_adj_4_8 = not_adj_4_8_tot.copy()\n",
    "not_adj_4_8[not_adj_4_8 <= ths] = 0\n",
    "\"\"\"adj_23_27[adj_23_27 >= ths] = 1\"\"\"\n",
    "print(\"Weighted number of rides : {}\".format(np.sum(not_adj_4_8) / 2))\n",
    "\n",
    "# Remove Self-loops\n",
    "for i in range(n_tot_station):\n",
    "    not_adj_4_8[i, i] = 0\n",
    "    \n",
    "not_g_4_8 = nx.from_numpy_matrix(not_adj_4_8)\n",
    "\n",
    "not_g_4_8.remove_nodes_from(out_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(wf_month['5']) / len(not_wf_month['5']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = 0\n",
    "for i in np.arange(adj_2_6.shape[0]):\n",
    "    s += adj_2_6_tot[i,i]\n",
    "\n",
    "print(s / np.sum(adj_2_6_tot))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = 0\n",
    "for i in np.arange(not_adj_4_8_tot.shape[0]):\n",
    "    s += not_adj_4_8_tot[i,i]\n",
    "\n",
    "print(s / np.sum(not_adj_4_8_tot))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = 0\n",
    "for i in np.arange(adj_4_8_tot.shape[0]):\n",
    "    s += adj_4_8_tot[i,i]\n",
    "\n",
    "print(s / np.sum(adj_4_8_tot))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "not_wf_month"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "thesis",
   "language": "python",
   "name": "thesis"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}