Add comments to cSNMF.py and cSNMF.ipynb

.gitignore

-.ipynb_checkpoints/*
-MultiplicativeAlgorithm/.ipynb_checkpoints/*
+*/.ipynb_checkpoints/
+*/__pycache__/
 MultiplicativeAlgorithm/D_trips.txt
 MultiplicativeAlgorithm/D_speeds.txt
 MultiplicativeAlgorithm/*.png

MultiplicativeAlgorithm/Extreme_Events.ipynb → MultiplicativeAlgorithm/ExtremeEvents.ipynb

@@ -163,7 +163,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.5"
+   "version": "3.6.6"
   }
  },
  "nbformat": 4,

MultiplicativeAlgorithm/Phase2.ipynb

@@ -2518,7 +2518,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.5"
+   "version": "3.6.6"
   }
  },
  "nbformat": 4,

MultiplicativeAlgorithm/Visualizations.py

@@ -315,3 +315,14 @@ def spikeyness_vs_error():
     plt.savefig('Spikeyness_vs_Error')
     return None
 
+def spikeyness2(W):
+    def spikeyness_column(column_number):
+        trend = W[:52*24*7, column_number].reshape(52,24*7)
+        if np.mean(W[:52*24*7, column_number]) == 0:
+            return np.nan
+        else:
+            return np.mean(np.nanstd(trend, axis=0))/np.mean(W[:, column_number])
+        
+    spikey_mean = np.mean([spikeyness_column(sig_number) for sig_number in range(W.shape[1])])
+    spikey_std = np.std([spikeyness(sig_number) for sig_number in range(W.shape[1])])
+    return spikey_mean, spikey_std
\ No newline at end of file

MultiplicativeAlgorithm/__init__.py

@@ -22,24 +22,4 @@ formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(messag
 fh.setFormatter(formatter)
 
 logger.addHandler(fh)
-logger.addHandler(sh)
-
-
-
-'''Select path name for files depending on whether you are working on your PC or cluster'''
-# path_for_data is an optional variable that defaults to 'path'. This is included to facilate a seperate path
-# address for large data files
-
-#path = './Data_Files/'         # Use this if working on PC
-#path_for_data = '../../'       # Use this if working on PC and if data files are in a different folder
-
-#path = './scratch/'           # Use this if working on cluster
-#path_for_data = './scratch/'  # Use this if working on cluster and if data files are in a different folder
-
-#filenames = config.generate_filenames(path, path_for_data)
-
-#logger.info('Filenames have been loaded')
-#[logger.debug(i) for i in filenames]
-
-#logger.info('Constants have been loaded')
-#[logger.debug(i) for i in dir(globals()['config'])];
\ No newline at end of file
+logger.addHandler(sh)
\ No newline at end of file

MultiplicativeAlgorithm/cSNMF.ipynb

@@ -4,12 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Run cSNMF on Trips data using multiplicative update rules\n",
-    "\n",
-    "**Constraint:** L1 norm of columns of W should be 1\n",
-    "\n",
-    "Get a copy of the data matrices in your local machine from the following links:\n",
-    " - https://uofi.box.com/s/yo60oe084d68obohgraek4weuaqtpgsp"
+    "# Run cSNMF on Trips data using multiplicative update rules"
    ]
   },
   {
@@ -18,13 +13,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from __init__ import *\n",
     "import numpy as np\n",
-    "import pandas as pd\n",
-    "import config\n",
-    "import cSNMF\n",
     "import matplotlib.pyplot as plt\n",
-    "%matplotlib inline"
+    "%matplotlib inline\n",
+    "\n",
+    "from __init__ import *  # initializes logger\n",
+    "import config           # contains all the global variables\n",
+    "import cSNMF            # contains update rules and factorization functions"
    ]
   },
   {
@@ -33,17 +28,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "## Read Full-Link data and prep for running NMF.\n",
-    "D = np.loadtxt('D_trips.txt')\n",
-    "logger.info('Full_link data has been read')\n",
+    "## Read Full-Link data and prep for running cSNMF.\n",
+    "if config.TRIPS:\n",
+    "    D = np.loadtxt('D_trips.txt')\n",
+    "    logger.info('Full_link data for trips has been read')\n",
+    "else:\n",
+    "    D = np.loadtxt('D_traveltimes.txt')\n",
+    "    logger.info('Full_link data for traveltimes has been read')\n",
     "\n",
-    "if config.SEEDED == 1:\n",
-    "    seed_W = 0; seed_H = 1\n",
-    "elif config.SEEDED == 0:\n",
-    "    seed_W = None; seed_H = None\n",
+    "# Seeded W and H will always return the same result. Turn off seeded for different results with each run.    \n",
+    "if config.SEEDED:\n",
+    "    seed_W, seed_H = 0, 1\n",
     "else:\n",
-    "    logger.critical('Seed value invalid. Needs to be 0 or 1. Check config.py!')\n",
-    "    quit()"
+    "    seed_W, seed_H = None, None"
    ]
   },
   {
@@ -52,10 +49,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# Factorizes D and imposes contraints.\n",
+    "# Returns numpy arrays W and H, and pandas dataframe results.\n",
     "W, H, results = cSNMF.factorize(D,\n",
-    "                                beta = 5000,\n",
+    "                                beta = config.BETA,\n",
     "                                rank = config.RANK,\n",
-    "                                max_iter = 600,\n",
     "                                seed_W = seed_W,\n",
     "                                seed_H = seed_H,\n",
     "                                debug = True,\n",
@@ -68,6 +66,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# Displays runtime statistics as a dataframe\n",
     "results"
    ]
   },
@@ -77,8 +76,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#np.savetxt('W_(seed_W = 10,seed_H = 21).txt', W)\n",
-    "#np.savetxt('H_(seed_W = 10,seed_H = 21).txt', H)"
+    "np.savetxt('W_trips.txt', W)\n",
+    "np.savetxt('H_trips.txt', H)"
    ]
   }
  ],
@@ -98,7 +97,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.3"
+   "version": "3.6.6"
   }
  },
  "nbformat": 4,

MultiplicativeAlgorithm/cSNMF.py

-''' This contains functions for cSNMF.'''
+''' This contains functions called by cSNMF.ipynb for constrained Sparse Nonnegative Matrix Factorization'''
 
 import numpy as np
 import pandas as pd
@@ -11,119 +11,100 @@ logger = logging.getLogger(__name__)
 
 
 ## Useful matrix operations:
-def m(A, B): return np.multiply(A,B)
 def div(A, B): return np.array([[A[i,j]/B[i,j] if A[i,j]>0 else 0 for j in range(A.shape[1])] for i in range(A.shape[0])])
-def d(A): return np.diag(np.diag(A))  # d(A) replaces all off-diagonal entries of A with 0.
-def N(A): return m(A, NONZEROS)       # NONZEROS is a global variable.
-
-
+def d(A):      return np.diag(np.diag(A))  # d(A) replaces all off-diagonal entries of A with 0.
+def N(A):      return A*NONZEROS           # NONZEROS is a global variable.
 
 
+# Calculates sparsity of H as in Low_Rank_Manhattan_Traffic.pdf
 def sparsity_metric(H):
     def sparsity_column(column):
         n = len(column)
-        if sum(column) == 0:
-            ratio = 1
-            logger.info('Column of zeros in H')
-        else:
-            ratio = np.linalg.norm(column, 1)/np.linalg.norm(column, 2)
+        assert sum(column) != 0, 'Column of zeros in H. Reduce rank'
+        ratio = np.linalg.norm(column, 1)/np.linalg.norm(column, 2)
         return (math.sqrt(n) - ratio)/(math.sqrt(n) - 1)
     return np.mean([sparsity_column(column) for column in H.T])
 
 
 def calculate_error(D, W, H):
-    norm_D = np.linalg.norm(D)
-    return np.linalg.norm(D - N(W@H))/norm_D*100
-
-def spikeyness(W):
-    def spikeyness_column(column_number):
-        trend = W[:52*24*7, column_number].reshape(52,24*7)
-        if np.mean(W[:52*24*7, column_number]) == 0:
-            return np.nan
-        else:
-            return np.mean(np.nanstd(trend, axis=0))/np.mean(W[:, column_number])
-        
-    spikey_mean = np.mean([spikeyness_column(sig_number) for sig_number in range(W.shape[1])])
-    spikey_std = np.std([spikeyness(sig_number) for sig_number in range(W.shape[1])])
-    return spikey_mean, spikey_std
-
-def zeros_in_H(H):
-    # Percentage of zero entries in H
-    # Use this on H_axed 
-    return sum(sum(H==0))/H.size
+    return np.linalg.norm(D - N(W@H))/NORM_D*100
 
 
+# Constraint imposed on the columns of W
 def impose_L1_constraint(W, H):
     W2 = W/sum(W)
     H2 = np.array([H[i]*sum(W)[i] for i in range(len(H))])
     return W2, H2
 
-
-def axe_H(H, tail_cutoff = 0.4):
+# Sets the bottom tail_cutoff fraction of entries in each column of H to 0.
+def axe_H(H, tail_cutoff = config.AXING_CUTOFF):
     H2 = H.copy()
     for i, column in enumerate(H.T):
-        sorted_column = sorted(column)
+        sorted_column  = sorted(column)
         cumulative_sum = it.accumulate(sorted_column)
-        kill_count = sum(list(cumulative_sum) < tail_cutoff*sum(column))
-        H2.T[i] = np.array([entry if entry >= sorted_column[kill_count] else 0 for entry in column])
+        kill_count     = sum(list(cumulative_sum) < tail_cutoff*sum(column))
+        H2.T[i]        = np.array([entry if entry >= sorted_column[kill_count] else 0 for entry in column])
     return H2
 
 
+'''Sort columns of W and rows of H according to the following rules --
+   - Give each link a total weight of 1.
+   - Distribute that weight across all of its signatures in proportion to the corresponding entry in H.
+   - For each signature, calculate popularity to be sum of weights across all links.
+   - Sort columns of W in decreasing order of popularity.
+   - Sort rows of H based on the same permutaion in order to preserve W@H'''
+# Sorting is done so we have a rough intuition for which signatures get used the most.
 def sort_WH(W, H):
-    H2 = axe_H(H)
-    weights = H2/sum(H2)
+    weights = H/sum(H)
     signature_popularities = sum(weights.T)
-    signature_popularities = sorted(enumerate(signature_popularities), key = lambda x : x[1], reverse = True)
+    signature_popularities = sorted(enumerate(signature_popularities),\
+                                    key = lambda x : x[1], reverse = True)
     W2 = W.copy()
-    H3 = H2.copy()
+    H2 = H.copy()
     for i in range(len(signature_popularities)):
         W2[:, i] = W[:, signature_popularities[i][0]]
-        H3[i] = H2[signature_popularities[i][0]]
-    try:
-        assert ( np.linalg.norm(W@H2 - W2@H3) / np.linalg.norm(W@H2) ) < 0.01
-    except AssertionError:
-        logger.critical('Sorting is not working correctly')
-        return W2, H3
-    
-    return W2, H3
+        H2[i] = H[signature_popularities[i][0]]
+        
+    assert ( np.linalg.norm(W@H - W2@H2) / np.linalg.norm(W@H) ) < 0.01,\
+        'Sorting is not working correctly'
+    return W2, H2
 
 
 
 def factorize(data_array,
-              init_W = None,
-              rank = config.RANK,
-              beta = None,
-              threshold = 0.5,
-              max_iter = 600,
-              seed_W = None,
-              seed_H = None,
-              log = logger,
-              debug = False,
-              axing = True,
-              update_W = True):
+              rank      = 50,
+              beta      = None,
+              threshold = 0.5,    # algorithm stops if change in W and H is less than threshold % 
+              max_iter  = 600,
+              seed_W    = None,
+              seed_H    = None,
+              log       = logger,
+              debug     = False,  # if debug, then return runtime-statistics 
+              axing     = True,   # if axing, then set small entries in H to zero
+              init_W    = None):  # optionally specify W as a numpy array and keep it unchanged in every iteration
 
     log.info('Rank= %s, Threshold= %s', rank, threshold)
 
     D = np.nan_to_num(data_array)
     eta = D.max()
     if beta is None:
-        beta = len(D.T)//2
-    global NONZEROS
-    NONZEROS = ~np.isnan(data_array)
+        beta = len(D.T)/2
     
-    JNN = np.ones((rank, rank))
+    global NONZEROS, NORM_D
+    NONZEROS = ~ np.isnan(data_array)
+    NORM_D   = np.linalg.norm(D)
     
+    JNN = np.ones((rank, rank))
     def update_W(W, H):
         Term_1 = D @ H.T 
         Term_2 = N(W @ H) @ H.T + eta*W
-        W = div(m(W, Term_1), Term_2)
+        W = div(W * Term_1, Term_2)
         return W
 
-    
     def update_H(W, H):
         Term_3 = W.T @ D
         Term_4 = W.T @ N(W @ H) + beta*(JNN @ H)
-        H = div(m(H, Term_3), Term_4)
+        H = div(H * Term_3, Term_4)
         return H
     
 
@@ -134,32 +115,30 @@ def factorize(data_array,
 
     
     logger.info('Initializing W and H...')
-    
     W_shape = (D.shape[0], rank)
     H_shape = (rank, D.shape[1])
     
     if init_W is None:
+        update_W = True
         np.random.seed(seed_W)
         W = np.random.uniform(low = 0.01, high = 1., size = W_shape)    
     else:
+        update_W = False
         W = init_W
         
     np.random.seed(seed_H)
     H = np.random.uniform(low = 0.01, high = 1., size = H_shape)    
     log.info('W, H chosen')
 
-    iterations = 0
+    
+    iterations   = 0
     column_names = ['error', 'sparsity', 'diff_W', 'diff_H', 'W_minmax',
                     'H_0th', 'H_25th', 'H_50th', 'H_75th', 'H_100th'] # For results DataFrame
     results = pd.DataFrame(columns = column_names)
-
-
-    diff_W = 100
-    diff_H = 100
-    
-    global norm_D
-    norm_D = np.linalg.norm(N(D))
+    diff_W  = 100
+    diff_H  = 100
     
+    # Run till stopping condition is met
     while abs(diff_W) + abs(diff_H) > threshold or iterations < 100:
         if iterations > max_iter:
             break
@@ -194,16 +173,14 @@ def factorize(data_array,
         if debug is True:
             H_0th, H_25th, H_50th, H_75th, H_100th = quartiles(H)
             W_minmax = (W.min(), W.max())
-        
             column_values = [error, sparsity, diff_W, diff_H, W_minmax,\
                              H_0th, H_25th, H_50th, H_75th, H_100th]
-
             results = results.append(dict(zip(column_names, column_values)), ignore_index = True)
-            
     
-    if axing:
-        W, H = impose_L1_constraint(W, H)
-        W, H = sort_WH(W, H)
+    
+    W, H = impose_L1_constraint(W, H)
+    if axing:  H = axe_H(H)
+    W, H = sort_WH(W, H)
     
     error = calculate_error(D, W, H)
     sparsity = sparsity_metric(H)

MultiplicativeAlgorithm/config.py

-''' * This file defines addresses to all data files, input files and output files.
-    * It also lists all the Global variables used in Non-Negative Matrix Factorization.
+''' * This file lists all the Global variables used in Non-Negative Matrix Factorization.
     * In case any of the Factorization paramaters need to be changed, they should be changed in this file instead of locally within each program.'''
 
+# Global variables
+RANK  = 50           # Rank for Matrix Factorization in Phase1.py
+TRIPS = True         # False reads travel_traveltimes dataset, True reads trips dataset.
+                     # Currently False doesn't work due to bad hyperparameters.    
+SEEDED = True        # Determines whether cSNMF.py uses seeded random matrices for W,H or not.
+BETA   = 5000        # Higher value penalizes lack of sparsity in H. (See Low_Rank_Manhattan_Traffic.pdf).
+                     # If unsure, then set BETA = None
+AXING_CUTOFF = 0.4   # Set between 0 and 1. Sets the bottom AXING_CUTOFF fraction of entries in each column of H to 0.
 
-def generate_filenames(path, path_for_data=None):
-    if path_for_data is None:
-        path_for_data = path
-    filenames = {# Raw data files
-                 'links':                path + 'links.csv',
-                 'nodes':                path + 'nodes.csv',
-                 'raw_data':             path_for_data + 'travel_times_2011.csv',
-                 'data_coo_form':        path_for_data + 'data_coo_form.txt',
-                 'data_trips':           path_for_data + 'data_trips.csv',
-                 'data_traveltimes':     path_for_data + 'data_travel_times.csv',
-                 'data_trips_transpose': path_for_data + 'data_trips_transpose.csv',
-    
-                 # Inputs and Outputs of Read_data.py
-                 'full_link_ids':         path + 'full_link_ids.txt',
-                 'empty_link_ids':        path + 'empty_link_ids.txt',
-                 'full_link_trips':       path + 'full_link_trips.json',
-                 'full_link_traveltimes': path + 'full_link_travel_times.json',
-                 'full_link_speeds':      path + 'full_link_speeds.json',
 
-                 # Inputs and Outputs of Phase1.py
-                 'random':             path + 'Random_numbers.txt',
-                 'W_trips':            path + 'W_trips.txt',
-                 'W_speeds':           path + 'W_speeds.csv',
-                 'W_trips_seed':       path + 'Seeded0,1/W_trips.csv',
-                 'W_speeds_seed':      path + 'Seeded0,1/W_speeds.csv',
-                 'H_trips':           path + 'H_trips.txt',
-                 'H_speeds':          path + 'H_speeds.csv',
-                 'H_trips_seed':      path + 'Seeded0,1/H_trips.csv',
-                 'H_speeds_seed':     path + 'Seeded0,1/H_speeds.csv',
-                 'H_trips_axed':      path + 'H_trips_axed.csv',
-                 'H_speeds_axed':     path + 'H_speeds_axed.csv',
-                 'H_trips_axed_seed': path + 'Seeded0,1/H_trips_axed.csv',
-                 'W_speeds_axed_seed': path + 'Seeded0,1/H_speeds_axed.csv'}
-    return filenames
+# Fixed global parameters for 2011 traffic dataset
+HOURS_IN_YEAR = 8760   # 24*365
+TOTAL_LINKS   = 260855 # Total Links aka Road-segments in NYC.
+FULL_LINKS    = 2302   # Links with less than 720 hours worth of missing data in a year.
+EMPTY_LINKS   = 234892 # Links with less that 720 hours worth of data in a year.
+MID_LINKS     = 23661  # TOTAL_LINKS - FULL_LINKS - EMPTY_LINKS
 
-
-'''Global variables'''
-
-RANK = 50            # Rank for Matrix Factorization in Phase1.py
-TRIPS = 1            # Boolean variable: 0 reads travel_speeds data, 1 reads number_of_trips data. 
-HOURS_IN_YEAR = 8760 # 24*365
-SEEDED = 1           # Boolean variable: 1 or 0 means Phase1.py is seeded or not respectively.
-
-TOTAL_LINKS = 260855 # Total Links aka Road-segments in NYC.
-FULL_LINKS = 2302    # Links with less than 720 hours worth of missing data in a year.
-EMPTY_LINKS = 234892 # Links with less that 720 hours worth of data in a year.
-MID_LINKS = 23661  # TOTAL_LINKS - FULL_LINKS - EMPTY_LINKS
 assert MID_LINKS == TOTAL_LINKS - FULL_LINKS - EMPTY_LINKS
\ No newline at end of file

ReadData/ReadData.ipynb

@@ -115,7 +115,7 @@
    "source": [
     "## Read from full_link_ids.txt\n",
     "\n",
-    "with open('../Multiplicative Algorithm/full_link_ids.txt', 'r') as readfile:\n",
+    "with open('../MultiplicativeAlgorithm/full_link_ids.txt', 'r') as readfile:\n",
     "    full_links = [int(line.strip()) for line in readfile]\n",
     "    print(len(full_links))"
    ]