Pro1_CQ_W6_e1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Network construction & Community detection--Twitter data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We are given the data from tweeter, recording users' ID and the time and location when they post anything. Our basic idea is to construct a network where we observe the relation among places by analyzing people's activities near each location, trying to detect grouping communities and distinguish weak links from strong ones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1,  load the data from twitter and categorize it by weekdays/weekend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "u'/Users/CQ/Documents/Project1'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Failed `CDLL(/Library/Frameworks/GEOS.framework/Versions/Current/GEOS)`\n"
     ]
    }
   ],
   "source": [
    "__author__ = \"Stephen\"\n",
    "__date__ = \"2016_03_27\"\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import networkx as nx #library supporting networks\n",
    "#make sure plots are embedded into the notebook\n",
    "#%pylab inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gp\n",
    "from datetime import datetime\n",
    "from datetime import timedelta, date\n",
    "import time\n",
    "import csv\n",
    "import community #Louvain method\n",
    "import operator\n",
    "import os\n",
    "\n",
    "PWD = '/Users/CQ/Documents/Project1'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Some information about these files:\n",
    "* twitter_file is the main database here, recording all tweeting activities from *Nov.4th 2015 to Feb.12th 2016*\n",
    "* USzipcode is the dataframe containing zipcodes in U.S. but there are some zipcodes missing, and we will only provide analysis based on the zipcodes we have here in the database\n",
    "* geo_NY is the shapefile for plotting, required by geopandas\n",
    "* NOTE: For the majority of missing zipcodes, personal experience tells me that those zipcodes are usually parks or forests, in other words, the places where postal service is not needed much"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "twitter_file = pd.read_csv(PWD + '/2016_02_18_sintetic.csv')\n",
    "USzipcode = pd.read_csv(PWD + '/USzipcode_XY.csv')\n",
    "zipPath = PWD + '/nyc-zip-code-tabulation-areas-polygons.geojson'\n",
    "geo_NY = gp.read_file(zipPath)[['geometry', 'postalCode']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Convert time stamp into time format and then categorize our data into four different datasets:\n",
    "1. All data, labeled  All days with **date_index 0**\n",
    "2. Data from weekdays, labeled  Weekdays with **date_index 1**\n",
    "3. Data from Saturdays, labeled  Saturdays with **date_index 2**\n",
    "4. Data from Sundays, labeled  Sundays with **date_index 3**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "TimeList = range(4)\n",
    "TimeList[0] = [time.localtime(x) for x in twitter_file.timestamp]\n",
    "TimeList[1] = [x for x in TimeList[0] if x.tm_wday < 5]\n",
    "TimeList[2] = [x for x in TimeList[0] if x.tm_wday == 5]\n",
    "TimeList[3] = [x for x in TimeList[0] if x.tm_wday == 6]\n",
    "twitter_file.iloc[:,0] = TimeList[0]\n",
    "\n",
    "twitter_file.columns = ['DateTime','Stamp','ZipCode','User','ID']\n",
    "USzipcode.columns = ['ZipCode','lat','lon']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Change the columns' names as above, compute and filter the results, drop all zeros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Clean_zipcode(dataset):\n",
    "    zipcode_all = list(sorted(set(dataset.ZipCode)))\n",
    "    zip_exist = [i for i in zipcode_all if i in USzipcode.ZipCode.values]\n",
    "    return dataset[dataset.ZipCode.isin(zip_exist)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Data = range(4)\n",
    "for i in range(4):\n",
    "    Data[i] = twitter_file[twitter_file.DateTime.isin(TimeList[i])]\n",
    "    Data[i] = Clean_zipcode(Data[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2, construct a network with the probability, or equivalently, the weight of link(a,b) defined as\n",
    "### $$link(a,b) = \\Sigma_c\\frac{t(c,a)*(t(c,b)-\\delta(a,b))}{T*(t(c)-1)}$$\n",
    "### where\n",
    "* *t(c,a)* denotes the total number of tweets that user *c* has posted at location (in our case, zipcode) *a*\n",
    "* $t(c) = \\Sigma_a t(c,a)$\n",
    "* $T = \\Sigma_c t(c)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Network_prob_build(dataset):\n",
    "    graph = nx.DiGraph()\n",
    "    places = [int(x[-5:]) for x in list(dataset.columns)[1:]]\n",
    "    graph.add_nodes_from(places)\n",
    "    T = dataset.iloc[:,1:].sum().sum()\n",
    "    tc = dataset.iloc[:,1:].sum(axis = 1)\n",
    "    \n",
    "    #ix_tc is the index for tc > 1, so we can avoid having zero terms in denominator: (t(c) - 1) == 0\n",
    "    ix_tc = tc > 1\n",
    "    #the formula is provided above\n",
    "    for a in places:\n",
    "        ix_a = places.index(a)\n",
    "        for b in places:\n",
    "            ix_b = places.index(b)\n",
    "            delta = (ix_a == ix_b)+0\n",
    "            w = 1.0 * dataset.iloc[:,ix_a+1][ix_tc] * (dataset.iloc[:,ix_b+1][ix_tc] - delta) / (tc[ix_tc]-1)\n",
    "            graph.add_edge(a,b, weight = sum(w)/T )\n",
    "    return graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Network_naive_build(date_index = 0):\n",
    "    Graph = nx.Graph()\n",
    "    Graph.add_nodes_from(Data[date_index].ZipCode.unique())\n",
    "    for i in range(len(Net1[date_index])):\n",
    "        w = Net1[date_index].iloc[i,:]\n",
    "        Graph.add_edge(w[0],w[1],weight = w[2])\n",
    "    return Graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compared with our original network, here we label it as Network2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#add unit weights to unweighted network nodes\n",
    "def make_directed_and_weighted(G):\n",
    "    DG = nx.DiGraph()\n",
    "    DG.add_nodes_from(G.nodes())\n",
    "    \n",
    "    for e in G.edges():\n",
    "        v1 = min(e[0],e[1])\n",
    "        v2 = max(e[0],e[1])\n",
    "        w = G[v1][v2]['weight']\n",
    "        #if it's a loop edge, we double its weight\n",
    "        if v1==v2:\n",
    "            DG.add_edge(v1,v2, weight = 2 * w)\n",
    "        #otherwise we set the weight of A->B equal to the weight of B->A\n",
    "        else:\n",
    "            DG.add_edge(v1,v2, weight = w)\n",
    "            DG.add_edge(v2,v1, weight = w)\n",
    "    return DG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# WARNING: NEXT CELL MAY TAKE MORE THAN 22 HOURS TO RUN, WE'VE EXTRACTED THOSE DATA TO EXTERNAL FILES FOR FUTURE USAGE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "for i in range(4):\n",
    "    Net2[i] = Normalized_weight(Data[i])\n",
    "    Net2[i].to_csv('/Users/CQ/Documents/Project1/Network2 for '+LABEL[i]+'days.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now we load the data for further analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Net1, Net2 = range(4), range(4)\n",
    "LABEL = ['All ','Week','Satur','Sun']\n",
    "for i in range(4):\n",
    "    Net1[i] = pd.read_csv(PWD + '/Network1/Network of '+LABEL[i]+'days.csv').iloc[:,1:]\n",
    "    Net2[i] = pd.read_csv(PWD + '/Network2/Network2 for '+LABEL[i]+'days.csv')\n",
    "    Net2[i].rename(columns={'Unnamed: 0':'UserID'}, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3, proceed to community detection, first construct the directed graphs based on the probabilistic methods, then plot the graphs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Convert these dataframes into directed graphs, and then perform combo partition and community plotting to see the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Graph, Graph1, Graph2 = range(4), range(4), range(4)\n",
    "for i in range(4):\n",
    "    Graph[i] = Network_naive_build(i)\n",
    "    Graph1[i] = make_directed_and_weighted(Graph[i])\n",
    "    Graph2[i] = Network_prob_build(Net2[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start with Louvain method if curious but not necessary, next try combo algorithm with Inf, 2, and 4 communities parititions.\n",
    "* We will run tests for graphs both with and without loop edges\n",
    "* Please note that you may have to downgrade the version of package networkx from 1.11 to 1.9.1 in order to successfully plot the graph  (--week 3, March 3rd 2016)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "workfolder = PWD + '/CommunityDetection'\n",
    "os.chdir(workfolder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#interface for running compiled combo over the network G given a maximal number of communities maxcom\n",
    "def getComboPartition(G, maxcom):\n",
    "    #save network in net format\n",
    "    nodes={}\n",
    "    nodenum={}\n",
    "    i=0\n",
    "    #create a dictionary transforming nodes to unique numbers\n",
    "    for n in G.nodes():\n",
    "        nodenum[n]=i\n",
    "        nodes[i]=n\n",
    "        i+=1\n",
    "    f = open('combo/temp.net', 'w')\n",
    "    f.write('*Arcs\\n')\n",
    "    for e in G.edges(data=True):\n",
    "        f.write('{0} {1} {2}\\n'.format(nodenum[e[0]],nodenum[e[1]],e[2]['weight']))\n",
    "    f.close()\n",
    "    #run combo\n",
    "    command=workfolder+'/combo/comboCPP combo/temp.net'\n",
    "    if maxcom<float('Inf'):\n",
    "        command=command+' {0}'.format(maxcom)  \n",
    "    os.system(command)\n",
    "    #read resulting partition\n",
    "    f = open('combo/temp_comm_comboC++.txt', 'r')\n",
    "    i=0\n",
    "    partition={}\n",
    "    for line in f:\n",
    "        partition[nodes[i]]=int(line)\n",
    "        i+=1\n",
    "    return partition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def modularity(G, partition):\n",
    "    #compute network modularity according to the given partitioning\n",
    "    nodes=G.nodes()\n",
    "    #compute node weights and total network weight\n",
    "    if G.is_directed():\n",
    "        w1=G.out_degree(weight='weight')\n",
    "        w2=G.in_degree(weight='weight')\n",
    "        T=1.0*sum([e[2]['weight'] for e in G.edges(data=True)])\n",
    "    else:\n",
    "        w1=G.degree(weight='weight')\n",
    "        w2=G.degree(weight='weight')\n",
    "        T=1.0*sum([(1+(e[0]!=e[1]))*e[2]['weight'] for e in G.edges(data=True)])\n",
    "    M=0.0 #start accumulating modularity score\n",
    "    for a in nodes:\n",
    "        for b in nodes:\n",
    "            #if (G.is_directed())|(b>=a):\n",
    "            if partition[a]==partition[b]: #if nodes belong to the same community\n",
    "                    #get edge weight\n",
    "                if G.has_edge(a,b):\n",
    "                    e=G[a][b]['weight']\n",
    "                else:\n",
    "                    e=0\n",
    "                M+=e/T-w1[a]*w2[b]/(T**2) #add modularity score for the considered edge \n",
    "    return M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def PartitionSorting(partition):\n",
    "    ClusterList = [ (j, partition.values().count(j)*1.0/len(partition)) for j in set(partition.values())]\n",
    "    ClusterList.sort(key = operator.itemgetter(1), reverse = 1)\n",
    "    ClusterRanking = [i[0] for i in ClusterList]\n",
    "    sorted_partition = {i:0 for i in partition.keys()}\n",
    "    for i in sorted_partition.keys():\n",
    "        sorted_partition[i] = ClusterRanking.index(partition[i]) + 1\n",
    "    return sorted_partition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N_CP4, P_CP4 = range(4), range(4)\n",
    "SN_CP4, SP_CP4 = range(4), range(4)\n",
    "for i in range(4):\n",
    "    N_CP4[i] = getComboPartition(Graph1[i], 4)\n",
    "    P_CP4[i] = getComboPartition(Graph2[i], 4)\n",
    "    SN_CP4[i] = PartitionSorting(N_CP4[i])\n",
    "    SP_CP4[i] = PartitionSorting(P_CP4[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def PlotMapPart(fig, partition, date_index=0, zips=geo_NY, key='postalCode', cmap='spectral', title='Community Detection',size=221, alpha=.7):\n",
    "    ax = fig.add_subplot(size + date_index)\n",
    "    y = {i:'Community %d'%(partition[i]+1) for i in partition.keys()}\n",
    "    p = pd.Series(y).reset_index().rename(columns={'index':'postalCode',0:'part'})\n",
    "    zips.postalCode = zips.postalCode.astype(int)\n",
    "    z = zips.merge(p, on=key, how='left')\n",
    "    #z.loc[z.part.notnull(),'part'] = z[z.part.notnull()].part.astype(int)\n",
    "    level = len(set(partition.values()))\n",
    "    Legend = ['Community %d'%(i+1) for i in range(level)]\n",
    "    z[z.part.notnull()].plot(column = 'part', categorical=1, ax=ax, alpha=alpha, cmap = cmap, legend=True)\n",
    "    plt.title(title, fontweight=\"bold\", size=15)\n",
    "    ax.axis('off')\n",
    "    #ax.legend(Legend, loc = 2)\n",
    "    plt.show()\n",
    "    plt.subplots_adjust(hspace = 0.1, wspace = 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interesting observation here:\n",
    "### The majority of the map is consistent regardless of the date, and both airports JFK&LaGuardia are labeled as the same color as midtown Manhattan for almost everyday, indicating high amount of Twitter activities. However, LaGuardia is more quite during weekends, while JFK has even higher activities on Saturdays due to the color consistency with lower Manhattan. I'm thinking maybe this results from 3 different possibilities:\n",
    "1. International air traffic level on Saturdays is relatively higher than usual, while domestic level drops\n",
    "2. Activity level in the city is relatively lower than usual on Saturdays\n",
    "3. Our data is incomplete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig1 = plt.figure(num = 1, figsize = (20,15))\n",
    "Title = ['All Data','Weekdays', 'Saturdays', 'Sundays']\n",
    "plt.suptitle('Figure 1(Naive Network)\\nCommunity Detection by Combo method--4', fontweight=\"bold\", size=18)\n",
    "for i in range(4):\n",
    "    PlotMapPart(fig1, SN_CP4[i], date_index=i, title = Title[i]+', modularity %f' %modularity(Graph1[i],SN_CP4[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some even more interesting observations here:\n",
    "* During work days, some places in Brooklyn, e.g. Brooklyn Heights have almost same level of activity as Manhattan, but it drops significantly during weekends\n",
    "* JFK's activity level is almost same as the Queens on Sundays, making me doubting even more about the aforementioned possibility in previous plot\n",
    "* Lower Manhatten is always ranking highest, while Bronx and Queens are always at the bottom, with Staten Island and lower Brooklyn being second-to-last"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig2 = plt.figure(num = 2, figsize = (20,15))\n",
    "plt.suptitle('Figure 2(Probabilistic network)\\nCommunity Detection by Combo method--4', fontweight=\"bold\", size=18)\n",
    "for i in range(4):\n",
    "    PlotMapPart(fig2, SP_CP4[i], date_index=i, title = Title[i]+', modularity %f' %modularity(Graph2[i],SP_CP4[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4, some data visualization and demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A) First glance, twitter users' daily activities, with weekends highlighted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def DailyActivity():\n",
    "    start_date = datetime.fromtimestamp(time.mktime(TimeList[0][0])).date()\n",
    "    end_date = datetime.fromtimestamp(time.mktime(TimeList[0][-1])).date()\n",
    "    day_count = (end_date-start_date).days\n",
    "    DateList = [start_date + timedelta(n) for n in range(day_count)]\n",
    "    T_ACT = {i:0 for i in DateList}\n",
    "    for single_day in DateList:\n",
    "        day_list = [x for x in Data[0].DateTime if datetime.fromtimestamp(time.mktime(x)).date() == single_day]\n",
    "        T_ACT[single_day] = len(day_list)\n",
    "    Thanksgiving = date(2015, 11, 26)\n",
    "    Christmas = date(2015, 12, 25)\n",
    "    NewYear = date(2016, 1, 1)\n",
    "    a = Thanksgiving - start_date\n",
    "    b = date(2015,12,1) - start_date\n",
    "    c = Christmas - start_date\n",
    "    d = NewYear - start_date\n",
    "    e = date(2016,2,1) - start_date\n",
    "    fig = plt.figure(figsize=(20,8))\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.set_title(\"Daily Twitter Activities\", fontweight=\"bold\", size=18)\n",
    "    bar = ax.bar(range(day_count), T_ACT.values(), color = ['g','g','g','r','r','g','g'] , align='center', width = 1.5)\n",
    "    ax.legend((bar[0],bar[4]), ('Weekdays', 'Weekends'), loc = 2, prop={'size':18})\n",
    "    ax.set_xticks([0,a.days, b.days, c.days, d.days, e.days])\n",
    "    ax.set_xticklabels(['Nov/4/2015','Thanksgiving','Dec/1/2015','Christmas','New Year','Feb/1/2016'],\n",
    "                       rotation=90, fontsize = 14)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This confirms that we do miss some data here as many days having \"zero\" activity which is clearly against common sense."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "DailyActivity()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### B) Top places (zip codes) with highest amount of tweets, marked in map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def TopZipCode(fig, date_index, top = 10):\n",
    "    ziplist = sort(Data[date_index].ZipCode.unique())\n",
    "    percentage_ziplist = {i:0 for i in ziplist}\n",
    "    T = len(Data[date_index])\n",
    "    for i in ziplist:\n",
    "        percentage_ziplist[i] = 100.0 * len(Data[date_index][Data[date_index].ZipCode == i])/T\n",
    "    sorted_ziplist = sorted(percentage_ziplist.items(), key=operator.itemgetter(1), reverse = 1)\n",
    "    new_ziplist = [x[0] for x in sorted_ziplist]\n",
    "    percentage = [x[1] for x in sorted_ziplist]\n",
    "    \n",
    "    title = ['All days','Weekdays', 'Saturdays', 'Sundays']\n",
    "    ax = fig.add_subplot(221 + date_index)\n",
    "    plt.title(title[date_index], fontweight=\"bold\", size=15)\n",
    "    ax.bar(range(top), percentage[:top], color='g')\n",
    "    ax.set_xticks(range(top))\n",
    "    ax.set_xticklabels(new_ziplist[:top], fontsize = 13-top/4, fontweight=\"bold\", rotation = 30)\n",
    "    plt.subplots_adjust(hspace = 0.2, wspace = 0.1)\n",
    "    return new_ziplist[:top]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### You can change the number of \"top\" as you want to see the difference (but make sure you've loaded all the prerequisite data in previous cells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "top = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "TopZip = range(4)\n",
    "fig4 = plt.figure(num = 4, figsize = (15,10))\n",
    "plt.suptitle('Figure 4\\nTop %d zipcodes with their percentage'%top, fontweight=\"bold\", size=18)\n",
    "for i in range(4):\n",
    "    TopZip[i] = TopZipCode(fig4, i, top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### C) Some demographic data of New York gathered in 2014"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Demographic = pd.read_csv(PWD + \"/Demographic_2014.csv\", low_memory=False)\n",
    "variable_names = list(Demographic.iloc[0,:])\n",
    "#first row is variable name\n",
    "Demographic = Demographic.iloc[1:,:]\n",
    "Demographic.index = range(len(Demographic))\n",
    "#get zipcode from the strange string\n",
    "zipcode = pd.DataFrame([int(Demographic.loc[i,'Name of Area'].split(\" \")[0]) for i in range(len(Demographic))])\n",
    "Demographic.columns = variable_names\n",
    "Demographic = pd.concat((zipcode,Demographic),axis=1)\n",
    "Demographic = Demographic.rename(columns={0:'zipcode'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Slice and simplify the dataframe according to our zipcode database so we can focus on New York instead of the whole country"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "Demo_NY = Demographic[Demographic.zipcode.isin(Data[0].ZipCode.unique())]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First plotting the sample base, which is the \"All days\" category, then we observe if there's any interesting features showing on each community based on our partition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "fig4 = plt.figure(num = 4, figsize = (20,15))\n",
    "plt.suptitle('Figure 4(Loop edges included)\\nCommunity Detection by Combo method--4', fontweight=\"bold\", size=18)\n",
    "for i in range(4):\n",
    "    CommPartShape(fig4, D_Graph2[i], CP4[i],i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig5 = plt.figure(figsize = (15,12))\n",
    "PlotMapPart(fig5, Graph2[0], SP_CP4[0], size = 111, date_index=0, alpha = 0.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def ShowFeature(feature_name, sorted_partition, LEGEND = None, title = 'Community Data', date_index = 0):\n",
    "    '''date_index has the similar function here: controlling the sample base\n",
    "    0--All days\n",
    "    1--Weekdays\n",
    "    2--Saturdays\n",
    "    3--Sundays\n",
    "    with 0 being the default'''\n",
    "    #please make sure to use sorted_partition so that ClusterRanking is no longer needed\n",
    "    partition = sorted_partition[date_index]\n",
    "    feature_len = len(feature_name)\n",
    "    data = Demo_NY.loc[:,feature_name]\n",
    "    data.fillna(value = 0, inplace=True)\n",
    "    data = data.applymap(int)\n",
    "    data['ZipCode'] = Demo_NY.zipcode\n",
    "    #zip_intersection = set(data.ZipCode) & set(partition.keys())\n",
    "    data = data[data.ZipCode.isin(set(partition.keys()))]\n",
    "    data['Part'] = [partition[i] for i in data.ZipCode]\n",
    "    grouped = data.iloc[:,:-2].groupby(data.Part)\n",
    "    SUM = grouped.sum()*1.0\n",
    "    RATIO = SUM.divide(sum(SUM, axis = 1), axis = 0)\n",
    "    STD = grouped.std()\n",
    "    MEAN = grouped.mean()\n",
    "    level = len(set(partition.values()))\n",
    "    Width = 3.0/feature_len/level\n",
    "    fig = plt.figure(figsize = (20, 9))\n",
    "    x = range(level)\n",
    "    LABEL = ['Community %d'%(i+1) for i in x]\n",
    "    bar_step = arange(0, feature_len)\n",
    "    Color = matplotlib.cm.spectral(np.linspace(0,1,feature_len+1))\n",
    "    ax = fig.add_subplot(111)\n",
    "    for i in range(feature_len):\n",
    "        plt.bar(x+bar_step[i]*Width, MEAN.iloc[:,i] ,width = Width, yerr = list(STD.iloc[:,i]), color = Color[i+1], align='center')\n",
    "    Title = ['ALL days','Weekdays', 'Saturdays', 'Sundays']\n",
    "    plt.title(title+' data on each community (Partition based on '+Title[date_index]+')', fontweight=\"bold\", size = 18)\n",
    "    y = [i+(feature_len-1)*0.5*Width for i in x]\n",
    "    ax.set_xticks(y)\n",
    "    ax.set_xticklabels(LABEL, fontsize = 14, fontweight=\"bold\")\n",
    "    if LEGEND != None:\n",
    "        plt.legend(LEGEND, fontsize = 18-feature_len/2, loc =2)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "population = ['SE_T009_002','SE_T009_003','SE_T009_004','SE_T009_005']\n",
    "legend_pop = ['Under 18', '18-34', '35-64', '65+']\n",
    "#'SE_T001_001' is the Total population\n",
    "ShowFeature(population,SP_CP4,legend_pop,'Population')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "race = ['SE_T013_002','SE_T013_003','SE_T013_004','SE_T013_005','SE_T013_006','SE_T013_007','SE_T013_008']\n",
    "legend_race = ['White', 'Black or African', 'American Indian and Alaska Native', 'Asian',\n",
    "            'Native Hawaiian and Other Pacific Islander', 'Some Other', 'Two or More races']\n",
    "#'SE_T013_001' is the Total population\n",
    "ShowFeature(race,SP_CP4,legend_race,'Race')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "households = ['SE_T017_002','SE_T017_003','SE_T017_004','SE_T017_007','SE_T018_002']\n",
    "legend_households = ['Family', 'Married-couple', 'Single-parent', 'Nofamily',\n",
    "                     'Household with one or more people under 18']\n",
    "#'SE_T017_001' is the Total Households\n",
    "ShowFeature(households,SP_CP4,legend_households,'Household')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "education = ['SE_T025_002','SE_T025_003','SE_T025_004','SE_T025_005','SE_T025_006','SE_T025_007','SE_T025_008']\n",
    "legend_edu = ['< High school','High school', 'Some college', 'Bechelor', 'Master', 'Professional school', 'Doctorate']\n",
    "# 'SE_T025_001'is the Total Population over 25\n",
    "ShowFeature(education,SP_CP4,legend_edu,'Education')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "insurance=['SE_T145_002','SE_T145_003','SE_T145_004','SE_T145_005']\n",
    "#'SE_T145_001' is the Total number, and \n",
    "legend_insur = ['No Coverage','With Some Coverage', 'Public Insurance', 'Private Insurance']\n",
    "ShowFeature(insurance,SP_CP4,legend_insur,'Insurance')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "born_place = ['SE_T133_002','SE_T133_004','SE_T133_005']\n",
    "legend_born = ['Native Born','Foreign Born: Naturalized Citizen', 'Foreign Born: Not a citizen']\n",
    "#'SE_T133_001' Total Population; 'SE_T133_003',Foreign Born\n",
    "ShowFeature(born_place, SP_CP4,legend_born,'Born Places')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "commute_time = ['SE_T129_003','SE_T129_004','SE_T129_005',\n",
    "                     'SE_T129_006','SE_T129_007','SE_T129_008','SE_T129_009','SE_T129_010']\n",
    "# 'SE_T129_001', Workers 16 Years and over; 'SE_T129_002', did not work at home\n",
    "legend_commute_time = ['Communte time < 10 minutes','10 to 19 minutes','20 to 29 minutes','30 to 39 minutes',\n",
    "                       '40 to 59 minutes','60 to 89 minutes','> 90 minutes', 'Work at home']\n",
    "ShowFeature(commute_time,SP_CP4,legend_commute_time,'Communte Time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "income=['SE_T056_002','SE_T056_003','SE_T056_004','SE_T056_005','SE_T056_006','SE_T056_007',\n",
    "       'SE_T056_008','SE_T056_009','SE_T056_010','SE_T056_011','SE_T056_012','SE_T056_013',\n",
    "       'SE_T056_014','SE_T056_015','SE_T056_016','SE_T056_017']\n",
    "#'SE_T056_001' is the number of total Households\n",
    "legend_income = ['< $10,000','\\$10,000 - $14,999', '\\$15,000 - $19,999', '\\$20,000 - $24,999', '\\$25,000 - $29,999',\n",
    "                 '\\$30,000 - $34,999', '\\$35,000 - $39,999', '\\$40,000 - $44,999', '\\$45,000 - $49,999',\n",
    "                 '\\$50,000 - $59,999','\\$60,000 - $74,999','\\$75,000 - $99,999', '\\$100,000 - $124,999',\n",
    "                 '\\$125,000 - $149,999','\\$150,000 - $199,999', '>= $200,000']\n",
    "ShowFeature(income, SP_CP4,legend_income,'Income')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "House_price =['SE_T100_002','SE_T100_003','SE_T100_004','SE_T100_005',\n",
    "                  'SE_T100_006','SE_T100_007','SE_T100_008','SE_T100_009','SE_T100_010']\n",
    "legend_house_price = ['House Price < $20,000','\\$20,000 to $49,999','\\$50,000 to $99,999','\\$100,000 to $149,999',\n",
    "                      '\\$150,000 to $299,999','\\$300,000 to $499,999','\\$500,000 to $749,999','\\$750,000 to $999,999',\n",
    "                      '> $1,000,000']\n",
    "#'SE_T100_001', Owner-occupied housing units\n",
    "ShowFeature(House_price,SP_CP4,legend_house_price,'Housing Price')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "Rent_price=['SE_T102_002','SE_T102_003','SE_T102_004','SE_T102_005',\n",
    "            'SE_T102_006','SE_T102_007','SE_T102_008','SE_T102_009']\n",
    "#'SE_T102_001', Renter-occupied housing units with cash rent\n",
    "legen_rent_price = ['Rent < $300','\\$300 to $599','\\$600 to $799','\\$800 to $999','\\$1,000 to $1,249',\n",
    "                    '\\$1,250 to $1,499','\\$1,500 to $1,999','Rent > $2,000']\n",
    "ShowFeature(Rent_price, SP_CP4,legen_rent_price,'Rent Price')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "employment = ['SE_T037_002','SE_T037_003']\n",
    "legend_emp = ['Employed', 'Unemployed']\n",
    "# 'SE_T037_001' is the total Civilian Population In Labor Force 16 Years And Over\n",
    "ShowFeature(employment,SP_CP4,legend_emp,'Employment')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "type_house = ['SE_T095_003','SE_T094_002','SE_T094_003']\n",
    "legend_housetype = ['Vacant Houses','Owner Occupied','Renter Occupied']\n",
    "#'SE_T094_001','Occupied houses'\n",
    "ShowFeature(type_house, SP_CP4,legend_housetype,'House Type')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "income_per_capital=['SE_T083_001']\n",
    "#legend_income_per_cap = ['Income Per Capital']\n",
    "#Per capita income (In 2014 Inflation adjusted dollars)\n",
    "ShowFeature(income_per_capital,SP_CP4,title='Income Per Capital')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Next we try page rank algorithm to see if it can provide anything interesting, we will compare our probabilistic network with naive networks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#produce symmetrized undirected version of a directed network\n",
    "def as_undirected(G):\n",
    "    UG = nx.Graph()\n",
    "    UG.add_nodes_from(G.nodes())\n",
    "    for e in G.edges():\n",
    "        UG.add_edge(e[0],e[1], weight = G[e[0]][e[1]]['weight'])\n",
    "    return UG  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "PR1, PR2, UGraph1_PR, UGraph2_PR = range(4), range(4), range(4), range(4)\n",
    "for i in range(4):\n",
    "    PR2[i] = as_undirected(Graph2[i])\n",
    "    UGraph1_PR[i] = list(nx.connected_component_subgraphs(Graph[i], copy=True))[0]\n",
    "    PR1[i] = nx.pagerank(UGraph1_PR[i],0.85)\n",
    "    UGraph2_PR[i] = list(nx.connected_component_subgraphs(PR2[i], copy=True))[0]\n",
    "    PR2[i] = nx.pagerank(UGraph2_PR[i],0.85)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def PlotPageRankPart(fig, partition, cluster = 4, date_index=0, zips=geo_NY, key='postalCode', cmap='spectral', title='Community Detection',size=221, alpha=.7):\n",
    "    a = np.arange(0.,100.,step = 100.0/cluster)\n",
    "    x = np.percentile(partition.values(),a)\n",
    "    #v2 = max(x.values())\n",
    "    #v1 = min(x.values())\n",
    "    y = {i:0 for i in partition.keys()}\n",
    "    for i in partition.keys():\n",
    "        for j in range(cluster):\n",
    "            if partition[i] >= x[j]:\n",
    "                y[i] = 'Community %d'%(j+1)\n",
    "\n",
    "    ax = fig.add_subplot(size + date_index)\n",
    "    p = pd.Series(y).reset_index().rename(columns={'index':'postalCode',0:'part'})\n",
    "    zips.postalCode = zips.postalCode.astype(int)\n",
    "    z = zips.merge(p, on=key, how='left')\n",
    "    z[z.part.notnull()].plot(column = 'part', ax=ax, alpha=alpha, cmap = cmap, categorical=1, legend = True)\n",
    "    plt.title(title, fontweight=\"bold\", size=15)\n",
    "    ax.axis('off')\n",
    "    plt.subplots_adjust(hspace = 0.1, wspace = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig20 = plt.figure(num = 20, figsize = (20,15))\n",
    "plt.suptitle('Naive network\\nCommunity Detection by Page Rank--4', fontweight=\"bold\", size=18)\n",
    "for i in range(4):\n",
    "    PlotPageRankPart(fig20, PR1[i], date_index=i, title = Title[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig21 = plt.figure(num = 21, figsize = (20,15))\n",
    "plt.suptitle('Probabilistic network\\nCommunity Detection by Page Rank--4', fontweight=\"bold\", size=18)\n",
    "for i in range(4):\n",
    "    PlotPageRankPart(fig21, PR2[i], date_index=i, title = Title[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig22 = plt.figure(num = 22, figsize = (20,15))\n",
    "plt.suptitle('Probabilistic network\\nCommunity Detection by Page Rank--10', fontweight=\"bold\", size=18)\n",
    "for i in range(4):\n",
    "    PlotPageRankPart(fig21, PR2[i], cluster=10, date_index=i, title = Title[i])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}