diff --git a/src/py/crankshaft/crankshaft/analysis_data_provider.py b/src/py/crankshaft/crankshaft/analysis_data_provider.py index cbc27bc..cd79000 100644 --- a/src/py/crankshaft/crankshaft/analysis_data_provider.py +++ b/src/py/crankshaft/crankshaft/analysis_data_provider.py @@ -65,3 +65,12 @@ class AnalysisDataProvider: return data except plpy.SPIError, err: plpy.error('Analysis failed: %s' % err) + + def get_gwr(self, params): + """fetch data for gwr analysis""" + query = pu.gwr_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/regression/gwr_cs.py b/src/py/crankshaft/crankshaft/regression/gwr_cs.py index b5b24d7..91d6d19 100644 --- a/src/py/crankshaft/crankshaft/regression/gwr_cs.py +++ b/src/py/crankshaft/crankshaft/regression/gwr_cs.py @@ -1,92 +1,96 @@ +""" + Geographically weighted regression +""" import numpy as np from gwr.base.gwr import GWR from gwr.base.sel_bw import Sel_BW import plpy import crankshaft.pysal_utils as pu import json +from crankshaft.analysis_data_provider import AnalysisDataProvider -def gwr(subquery, dep_var, ind_vars, bw=None, - fixed=False, kernel='bisquare'): - """ - 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' - """ +class GWR: + def __init__(self, analysis_provider=None): + if analysis_provider: + self.analysis_provider = analysis_provider + else: + self.analysis_provider = AnalysisDataProvider() - # query_result = subquery - params = {'geom_col': 'the_geom', - 'id_col': 'cartodb_id', - 'subquery': subquery, - 'dep_var': dep_var, - 'ind_vars': ind_vars} + def gwr(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' + """ - try: - query = pu.gwr_query(params) - plpy.notice(query) - query_result = plpy.execute(query) - except plpy.SPIError, err: - plpy.notice(query) - plpy.error('Analysis failed: %s' % err) + params = {'geom_col': geom_col, + 'id_col': id_col, + 'subquery': subquery, + 'dep_var': dep_var, + 'ind_vars': ind_vars} - # unique ids and variable names list - rowid = np.array(query_result[0]['rowid'], dtype=np.int) + # retrieve data + query_result = self.analysis_data_provider.get_gwr(params) - # TODO: should x, y be centroids? point on surface? - # lat, long coordinates - x = np.array(query_result[0]['x'], dtype=float) - y = np.array(query_result[0]['y'], dtype=float) - coords = zip(x, y) + # unique ids and variable names list + rowid = np.array(query_result[0]['rowid'], dtype=np.int) - # extract dependent variable - Y = np.array(query_result[0]['dep_var'], dtype=float).reshape((-1, 1)) + # TODO: should x, y be centroids? point on surface? + # lat, long coordinates + x = np.array(query_result[0]['x'], dtype=float) + y = np.array(query_result[0]['y'], dtype=float) + coords = zip(x, y) - n = Y.shape[0] - k = len(ind_vars) - X = np.zeros((n, k)) + # extract dependent variable + Y = np.array(query_result[0]['dep_var'], dtype=float).reshape((-1, 1)) - for attr in range(0, k): - attr_name = 'attr' + str(attr + 1) - X[:, attr] = np.array( - query_result[0][attr_name], dtype=float).flatten() + n = Y.shape[0] + k = len(ind_vars) + X = np.zeros((n, k)) - # add intercept variable name - ind_vars.insert(0, 'intercept') + # extract query result + for attr in range(0, k): + attr_name = 'attr' + str(attr + 1) + X[:, attr] = np.array( + query_result[0][attr_name], dtype=float).flatten() - # calculate bandwidth if none is supplied - plpy.notice(str(bw)) - if bw is None: - bw = Sel_BW(coords, Y, X, - fixed=fixed, kernel=kernel).search() - plpy.notice(str(bw)) - model = GWR(coords, Y, X, bw, - fixed=fixed, kernel=kernel).fit() + # add intercept variable name + ind_vars.insert(0, 'intercept') - # TODO: iterate from 0, n-1 and fill objects like this, for a - # column called coeffs: - # {'pctrural': ..., 'pctpov': ..., ...} - # Follow the same structure for other outputs + # calculate bandwidth if none is supplied + plpy.notice(str(bw)) + if bw is None: + bw = Sel_BW(coords, Y, X, + fixed=fixed, kernel=kernel).search() + plpy.notice(str(bw)) + model = GWR(coords, Y, X, bw, + fixed=fixed, kernel=kernel).fit() - coefficients = [] - stand_errs = [] - t_vals = [] - predicted = model.predy.flatten() - residuals = model.resid_response - r_squared = model.localR2.flatten() - bw = np.repeat(float(bw), n) + # containers for outputs + coeffs = [] + stand_errs = [] + t_vals = [] - for idx in xrange(n): - 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] + # extracted model information + predicted = model.predy.flatten() + residuals = model.resid_response + r_squared = model.localR2.flatten() + bw = np.repeat(float(bw), n) + + # create lists of json objs for model outputs + for idx in xrange(n): + coeffs.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)})) - t_vals.append(json.dumps({var: model.tvalues[idx, k] - for k, var in enumerate(ind_vars)})) - plpy.notice(str(zip(coefficients, stand_errs, t_vals, - predicted, residuals, r_squared, rowid, bw))) - return zip(coefficients, stand_errs, t_vals, - predicted, residuals, r_squared, rowid, bw) + return zip(coeffs, stand_errs, t_vals, + predicted, residuals, r_squared, rowid, bw)