From d95ca54cdc42df6bf4f02c9c1a6d8cf2db904546 Mon Sep 17 00:00:00 2001
From: Andy Eschbacher <andy.eschbacher@gmail.com>
Date: Thu, 5 Jan 2017 14:22:40 -0500
Subject: [PATCH] updating prediction to data provider

---
 .../crankshaft/analysis_data_provider.py      |  8 ++
 .../crankshaft/pysal_utils/pysal_utils.py     |  1 +
 .../crankshaft/regression/gwr_cs.py           | 94 ++++++++++++++++++-
 3 files changed, 101 insertions(+), 2 deletions(-)

diff --git a/src/py/crankshaft/crankshaft/analysis_data_provider.py b/src/py/crankshaft/crankshaft/analysis_data_provider.py
index cd79000..84b0bea 100644
--- a/src/py/crankshaft/crankshaft/analysis_data_provider.py
+++ b/src/py/crankshaft/crankshaft/analysis_data_provider.py
@@ -74,3 +74,11 @@ class AnalysisDataProvider:
             return query_result
         except plpy.SPIError, err:
             plpy.error('Analysis failed: %s' % err)
+    def get_gwr_predict(self, params):
+        """fetch data for gwr predict"""
+        query = pu.gwr_predict_query(params)
+        try:
+            query_result = plpy.execute(query)
+            return query_result
+        except plpy.SPIError, err:
+            plpy.error('Analysis failed: %s' % err)
diff --git a/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py b/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py
index 765beaa..073ab27 100644
--- a/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py
+++ b/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py
@@ -222,6 +222,7 @@ def gwr_query(params):
 
     return query.format(**params).strip()
 
+
 def gwr_predict_query(params):
     """
     GWR query
diff --git a/src/py/crankshaft/crankshaft/regression/gwr_cs.py b/src/py/crankshaft/crankshaft/regression/gwr_cs.py
index 3b812cf..2d29ed2 100644
--- a/src/py/crankshaft/crankshaft/regression/gwr_cs.py
+++ b/src/py/crankshaft/crankshaft/regression/gwr_cs.py
@@ -39,8 +39,7 @@ class GWR:
         # unique ids and variable names list
         rowid = np.array(query_result[0]['rowid'], dtype=np.int)
 
-        # TODO: should x, y be centroids? point on surface?
-        #       lat, long coordinates
+        # x, y are centroids of input geometries
         x = np.array(query_result[0]['x'], dtype=float)
         y = np.array(query_result[0]['y'], dtype=float)
         coords = zip(x, y)
@@ -90,3 +89,94 @@ class GWR:
 
         return zip(coeffs, stand_errs, t_vals,
                    predicted, residuals, r_squared, bw, rowid)
+
+    def gwr_predict(self, subquery, dep_var, ind_vars,
+                    bw=None, fixed=False, kernel='bisquare',
+                    geom_col='the_geom', id_col='cartodb_id'):
+        """
+        subquery: 'select * from demographics'
+        dep_var: 'pctbachelor'
+        ind_vars: ['intercept', 'pctpov', 'pctrural', 'pctblack']
+        bw: value of bandwidth, if None then select optimal
+        fixed: False (kNN) or True ('distance')
+        kernel: 'bisquare' (default), or 'exponential', 'gaussian'
+        """
+
+        params = {'geom_col': geom_col,
+                  'id_col': id_col,
+                  'subquery': subquery,
+                  'dep_var': dep_var,
+                  'ind_vars': ind_vars}
+
+        query_result = self.data_provider.get_gwr_predict(params)
+
+        # unique ids and variable names list
+        rowid = np.array(query_result[0]['rowid'], dtype=np.int)
+
+        x = np.array(query_result[0]['x'])
+        y = np.array(query_result[0]['y'])
+        coords = np.array(zip(x, y))
+
+        # extract dependent variable
+        Y = np.array(query_result[0]['dep_var']).reshape((-1, 1))
+
+        n = Y.shape[0]
+        k = len(ind_vars)
+        X = np.zeros((n, k))
+
+        for attr in range(0, k):
+            attr_name = 'attr' + str(attr + 1)
+            X[:, attr] = np.array(
+              query_result[0][attr_name]).flatten()
+
+        # add intercept variable name
+        ind_vars.insert(0, 'intercept')
+
+        # split data into "training" and "test" for predictions
+        # create index to split based on null y values
+        train = np.where(Y != np.array(None))[0]
+        test = np.where(Y == np.array(None))[0]
+
+        if len(test) < 1:
+            plpy.error('No rows flagged for prediction: verify that rows '
+                       'denoting prediction locations have a dependent '
+                       'variable value of null')
+
+        # split dependent variable (only need training which is non-Null's)
+        Y_train = Y[train].reshape((-1, 1))
+        Y_train = Y_train.astype(np.float)
+
+        # split coords
+        coords_train = coords[train]
+        coords_test = coords[test]
+
+        # split explanatory variables
+        X_train = X[train]
+        X_test = X[test]
+
+        # calculate bandwidth if none is supplied
+        if bw is None:
+            bw = Sel_BW(coords_train, Y_train, X_train,
+                        fixed=fixed, kernel=kernel).search()
+
+        # estimate model and predict at new locations
+        model = GWR(coords_train, Y_train, X_train, bw,
+                    fixed=fixed, kernel=kernel).predict(coords_test, X_test)
+
+        coefficients = []
+        stand_errs = []
+        t_vals = []
+        r_squared = model.localR2.flatten()
+        predicted = model.predy.flatten()
+
+        m = len(model.predy)
+        for idx in xrange(m):
+            coefficients.append(json.dumps({var: model.params[idx, k]
+                                            for k, var in enumerate(ind_vars)}))
+            stand_errs.append(json.dumps({var: model.bse[idx, k]
+                                          for k, var in enumerate(ind_vars)}))
+            t_vals.append(json.dumps({var: model.tvalues[idx, k]
+                                      for k, var in enumerate(ind_vars)}))
+
+        return zip(coefficients, stand_errs, t_vals,
+                   r_squared, predicted, rowid)