Merge 0309207 into e97dadf

pandana/network.py

@@ -18,9 +18,15 @@
     "SUM": 0,
     "AVE": 1,
     "AVERAGE": 1,
-    "STD": 5,
-    "STDDEV": 5,
-    "COUNT": 6
+    "MIN": 2,
+    "25PCT": 3,
+    "MEDIAN": 4,
+    "MED": 4,
+    "75PCT": 5,
+    "MAX": 6,
+    "STD": 7,
+    "STDDEV": 7,
+    "COUNT": 8
 }
 
 DECAYS = {
@@ -282,7 +288,9 @@ def aggregate(self, distance, type="sum", decay="linear", imp_name=None,
             The maximum distance to aggregate data within
         type : string
             The type of aggregation, can be one of "ave", "sum", "std",
-            and "count"
+            "count", and now "min", "25pct", "median", "75pct", and "max" will
+            compute the associated quantiles.  (Quantiles are computed by
+            sorting so might be slower than the others.)
         decay : string
             The type of decay to apply, which makes things that are further
             away count less in the aggregation - must be one of "linear",

pandana/tests/test_pandana.py

@@ -1,6 +1,7 @@
 import os.path
 
 import numpy as np
+from numpy.testing import assert_allclose
 import pandas as pd
 import pytest
 from pandas.util import testing as pdt
@@ -52,6 +53,63 @@ def test_create_network(sample_osm):
     pass
 
 
+def test_agg_variables_accuracy(sample_osm):
+    net = sample_osm
+
+    # test accuracy compared to pandas functions
+
+    net.set(pd.Series(net.node_ids))
+    s = net.aggregate(10000, type="COUNT")
+    # not all the nodes in the sample network are connected
+    # get the nodes in the largest connected subgraph
+    # from printing the result out I know the largest subgraph has
+    # 477 nodes in the sample data
+    connected_nodes = s[s == 477].index.values
+
+    ssize = 50
+    r = random_data(ssize)
+    nodes = pd.Series(np.random.choice(connected_nodes, ssize))
+    net.set(nodes, variable=r)
+
+    s = net.aggregate(100000, type="COUNT").loc[connected_nodes]
+    assert s.unique().size == 1
+    assert s.iloc[0] == 50
+
+    s = net.aggregate(100000, type="AVE").loc[connected_nodes]
+    assert s.describe()['std'] < .01  # assert almost equal
+    assert_allclose(s.mean(), r.mean(), atol=1e-3)
+
+    s = net.aggregate(100000, type="MIN").loc[connected_nodes]
+    assert s.describe()['std'] < .01  # assert almost equal
+    assert_allclose(s.mean(), r.min(), atol=1e-3)
+
+    s = net.aggregate(100000, type="MAX").loc[connected_nodes]
+    assert s.describe()['std'] < .01  # assert almost equal
+    assert_allclose(s.mean(), r.max(), atol=1e-3)
+
+    r.sort()
+
+    s = net.aggregate(100000, type="MEDIAN").loc[connected_nodes]
+    assert s.describe()['std'] < .01  # assert almost equal
+    assert_allclose(s.mean(), r.iloc[25], atol=1e-2)
+
+    s = net.aggregate(100000, type="25PCT").loc[connected_nodes]
+    assert s.describe()['std'] < .01  # assert almost equal
+    assert_allclose(s.mean(), r.iloc[12], atol=1e-2)
+
+    s = net.aggregate(100000, type="75PCT").loc[connected_nodes]
+    assert s.describe()['std'] < .01  # assert almost equal
+    assert_allclose(s.mean(), r.iloc[37], atol=1e-2)
+
+    s = net.aggregate(100000, type="SUM").loc[connected_nodes]
+    assert s.describe()['std'] < .05  # assert almost equal
+    assert_allclose(s.mean(), r.sum(), atol=1e-2)
+
+    s = net.aggregate(100000, type="STD").loc[connected_nodes]
+    assert s.describe()['std'] < .01  # assert almost equal
+    assert_allclose(s.mean(), r.std(), atol=1e-2)
+
+
 def test_agg_variables(sample_osm):
     net = sample_osm
 

src/accessibility.cpp

@@ -1,3 +1,4 @@
+#include <algorithm>
 #include "accessibility.h"
 
 namespace MTC {
@@ -9,8 +10,7 @@ namespace MTC {
 		}
 
 		double Accessibility::compute_centrality(
-										int srcnode, DistanceVec &distances,
-									    int gno) {
+			int srcnode, DistanceVec &distances, int gno) {
 			if(distances.size() < 3) return 0.0;
 			int cnt = 0;
 			for(int i = 0 ; i < distances.size() ; i++) {
@@ -20,7 +20,7 @@ namespace MTC {
 
 					if(target <= source) continue;
 					std::vector<NodeID> path = ga[gno]->Route(
-											source,target,omp_get_thread_num());
+						source, target, omp_get_thread_num());
 
 					for(int i = 0 ; i < path.size() ; i++) {
 						if(path[i] == srcnode) cnt++;
@@ -74,9 +74,8 @@ namespace MTC {
 
 			assert(cat >= 0 && cat < POI_MAXVAL);
 
-			DistanceMap distances = ga[gno]->NearestPOI(cat,srcnode,maxradius,
-			                                            number,
-														omp_get_thread_num());
+			DistanceMap distances = ga[gno]->NearestPOI(cat, srcnode,
+				maxradius, number, omp_get_thread_num());
 			std::vector<float> ret;
 
 			accessibility_vars_t &vars = accessibilityVarsForPOIs[cat];
@@ -107,12 +106,12 @@ namespace MTC {
 			            std::vector<float> ( number ));
 			#pragma omp parallel for
 			for(int i = 0 ; i < numnodes ; i++) {
-				std::vector<float> d = findNearestPOIs(i, maxradius, number,
-				                                       cat, gno);
+				std::vector<float> d = findNearestPOIs(i, maxradius,
+					number, cat, gno);
 				for(int j = 0 ; j < number ; j++) {
-                    if(j < d.size()) dists[i][j] = d[j];
-                    else dists[i][j] = -1;
-                }
+                    			if(j < d.size()) dists[i][j] = d[j];
+                    			else dists[i][j] = -1;
+                		}
 			}
 			return dists;
 		}
@@ -128,7 +127,7 @@ namespace MTC {
 			{
 			#pragma omp for schedule(guided)
 			for(int i = 0 ; i < numnodes ; i++) {
-                for(int j = 0 ; j < ga.size() ; j++) {
+                		for(int j = 0 ; j < ga.size() ; j++) {
 				    ga[j]->Range(i,radius,omp_get_thread_num(),dms[j][i]);
 				}
 			}
@@ -172,6 +171,50 @@ namespace MTC {
 			}
         }
 
+        double
+		Accessibility::quantileAccessibilityVariable(DistanceVec &distances,
+			accessibility_vars_t &vars, float quantile, float radius) {
+
+			// first iterate through nodes in order to get count of items
+			int cnt = 0;
+
+			for (int i = 0 ; i < distances.size() ; i++) {
+
+				int nodeid = distances[i].first;
+				double distance = distances[i].second;
+
+				if(distance > radius) continue;
+
+				cnt += vars[nodeid].size();
+			}
+
+			if(cnt == 0) return -1;
+
+			std::vector<float> vals(cnt);
+
+			// make a second pass to put items in a single array for sorting
+			for (int i = 0, cnt = 0 ; i < distances.size() ; i++) {
+
+				int nodeid = distances[i].first;
+				double distance = distances[i].second;
+
+				if(distance > radius) continue;
+
+				// and then iterate through all items at the node
+				for(int j = 0 ; j < vars[nodeid].size() ; j++)
+					vals[cnt++] = vars[nodeid][j];
+			}
+
+			std::sort(vals.begin(), vals.end());
+
+			int ind = (int)(vals.size() * quantile);
+
+			if(quantile <= 0.0) ind = 0;
+			if(quantile >= 1.0) ind = vals.size()-1;
+
+			return vals[ind];
+		}
+
 		double
 		Accessibility::aggregateAccessibilityVariable(int srcnode, float radius,
 									accessibility_vars_t &vars,
@@ -194,6 +237,19 @@ namespace MTC {
 			}
 
 			if(distances.size() == 0) return -1;
+
+			if(aggtyp == AGG_MIN) {
+				return this->quantileAccessibilityVariable(distances, vars, 0.0, radius);
+			} else if (aggtyp == AGG_25PERCENTILE) {
+				return this->quantileAccessibilityVariable(distances, vars, 0.25, radius);
+			} else if (aggtyp == AGG_MEDIAN) {
+				return this->quantileAccessibilityVariable(distances, vars, 0.5, radius);
+			} else if (aggtyp == AGG_75PERCENTILE) {
+				return this->quantileAccessibilityVariable(distances, vars, 0.75, radius);
+			} else if (aggtyp == AGG_MAX) {
+				return this->quantileAccessibilityVariable(distances, vars, 1.0, radius);
+			}
+
 		    if(aggtyp == AGG_STDDEV) decay = DECAY_FLAT;
 
 			int cnt = 0;

src/accessibility.h

@@ -18,8 +18,10 @@ namespace MTC {
 			AGG_SUM,
 			AGG_AVE,
 			AGG_MIN,
-			AGG_MAX,
+			AGG_25PERCENTILE,
 			AGG_MEDIAN,
+			AGG_75PERCENTILE,
+			AGG_MAX,
 			AGG_STDDEV,
 			AGG_COUNT,
 			AGG_MAXVAL
@@ -49,6 +51,9 @@ namespace MTC {
 		        accessibility_vars_t &vars, aggregation_types_t aggtyp,
 		        decay_func_t gravity_func, int graphno=0);
 
+			double quantileAccessibilityVariable(DistanceVec &distances,
+				accessibility_vars_t &vars, float quantile, float radius);
+
 			// computes the accessibility for every node in the network
 			std::vector<double>
 			getAllAggregateAccessibilityVariables(float radius, int index,