updating to run in crankshaft. Output still wrong but getting there

CONTRIBUTING.md

@@ -45,8 +45,8 @@ source envs/dev/bin/activate
 
 Update extension in a working database with:
 
-* `ALTER EXTENSION crankshaft UPDATE TO 'current';`
-  `ALTER EXTENSION crankshaft UPDATE TO 'dev';`
+* `ALTER EXTENSION crankshaft VERSION TO 'current';`
+  `ALTER EXTENSION crankshaft VERSION TO 'dev';`
 
 Note: we keep the current development version install as 'dev' always;
 we update through the 'current' alias to allow changing the extension
@@ -58,10 +58,7 @@ should be dropped manually before the update.
 If the extension has not previously been installed in a database,
 it can be installed directly with:
 
-* `CREATE EXTENSION IF NOT EXISTS plpythonu;`
-  `CREATE EXTENSION IF NOT EXISTS postgis;`
-  `CREATE EXTENSION IF NOT EXISTS cartodb;`
-  `CREATE EXTENSION crankshaft WITH VERSION 'dev';`
+* `CREATE EXTENSION crankshaft WITH VERSION 'dev';`
 
 Note: the development extension uses the development python virtual
 environment automatically.

doc/02_moran.md

@@ -1,102 +1,4 @@
-## Name
-
-CDB_AreasOfInterest -- returns a table with a cluster/outlier classification, the significance of a classification, an autocorrelation statistic (Local Moran's I), and the geometry id for each geometry in the original dataset.
-
-## Synopsis
-
-```sql
-table(numeric moran_val, text quadrant, numeric significance, int ids, numeric column_values) CDB_AreasOfInterest(text query, text column_name)
-
-table(numeric moran_val, text quadrant, numeric significance, int ids, numeric column_values) CDB_AreasOfInterest(text query, text column_name, int permutations, text geom_column, text id_column, text weight_type, int num_ngbrs)
-```
-
-## Description
-
-CDB_AreasOfInterest is a table-returning function that classifies the geometries in a table by an attribute and gives a significance for that classification. This information can be used to find "Areas of Interest" by using the correlation of a geometry's attribute with that of its neighbors. Areas can be clusters, outliers, or neither (depending on which significance value is used).
-
-Inputs:
-
-* `query` (required): an arbitrary query against tables you have access to (e.g., in your account, shared in your organization, or through the Data Observatory). This string must contain the following columns: an id `INT` (e.g., `cartodb_id`), geometry (e.g., `the_geom`), and the numeric attribute which is specified in `column_name`
-* `column_name` (required): column to perform the area of interest analysis tool on. The data must be numeric (e.g., `float`, `int`, etc.)
-* `permutations` (optional): used to calculate the significance of a classification. Defaults to 99, which is sufficient in most situations.
-* `geom_column` (optional): the name of the geometry column. Data must be of type `geometry`.
-* `id_column` (optional): the name of the id column (e.g., `cartodb_id`). Data must be of type `int` or `bigint` and have a unique condition on the data.
-* `weight_type` (optional): the type of weight used for determining what defines a neighborhood. Options are `knn` or `queen`.
-* `num_ngbrs` (optional): the number of neighbors in a neighborhood around a geometry. Only used if `knn` is chosen above.
-
-Outputs:
-
-* `moran_val`: underlying correlation statistic used in analysis
-* `quadrant`: human-readable interpretation of classification
-* `significance`: significance of classification (closer to 0 is more significant)
-* `ids`: id of original geometry (used for joining against original table if desired -- see examples)
-* `column_values`: original column values from `column_name`
-
-Availability: crankshaft v0.0.1 and above
-
-## Examples
-
-```sql
-SELECT
-  t.the_geom_webmercator,
-  t.cartodb_id,
-  aoi.significance,
-  aoi.quadrant As aoi_quadrant
-FROM
-  observatory.acs2013 As t
-JOIN
-  crankshaft.CDB_AreasOfInterest('SELECT * FROM observatory.acs2013',
-                                 'gini_index')
-```
-
-## API Usage
-
-Example
-
-```text
-http://eschbacher.cartodb.com/api/v2/sql?q=SELECT * FROM crankshaft.CDB_AreasOfInterest('SELECT * FROM observatory.acs2013','gini_index')
-```
-
-Result
-```json
-{
-  time: 0.120,
-  total_rows: 100,
-  rows: [{
-    moran_vals: 0.7213,
-    quadrant: 'High area',
-    significance: 0.03,
-    ids: 1,
-    column_value: 0.22
-  },
-  {
-    moran_vals: -0.7213,
-    quadrant: 'Low outlier',
-    significance: 0.13,
-    ids: 2,
-    column_value: 0.03
-  },
-  ...
-  ]
-}
-```
-
-## See Also
-
-crankshaft's areas of interest functions:
-
-* [CDB_AreasOfInterest_Global]()
-* [CDB_AreasOfInterest_Rate_Local]()
-* [CDB_AreasOfInterest_Rate_Global]()
-
-
-PostGIS clustering functions:
-
-* [ST_ClusterIntersecting](http://postgis.net/docs/manual-2.2/ST_ClusterIntersecting.html)
-* [ST_ClusterWithin](http://postgis.net/docs/manual-2.2/ST_ClusterWithin.html)
-
-
--- removing below, working into above
+### Moran's I
 
 #### What is Moran's I and why is it significant for CartoDB?
 

doc/docs_template.md

@@ -1,24 +0,0 @@
-
-## Name
-
-## Synopsis
-
-## Description
-
-Availability: v...
-
-## Examples
-
-```SQL
--- example of the function in use
-SELECT cdb_awesome_function(the_geom, 'total_pop')
-FROM table_name
-```
-
-## API Usage
-
-_asdf_
-
-## See Also
-
-_Other function pages_

src/pg/sql/50_contours.sql

@@ -0,0 +1,32 @@
+
+CREATE OR REPLACE FUNCTION
+  _CDB_Contours (
+      subquery TEXT,
+      grid_size NUMERIC DEFAULT 100,
+      bandwidth NUMERIC DEFAULT 0.0001,
+      levels NUMERIC[] DEFAULT null
+      )
+RETURNS table (level Numeric, geom_text text )
+AS $$
+  from crankshaft.contours import cdb_generate_contours
+  # TODO: use named parameters or a dictionary
+  return cdb_generate_contours(subquery, grid_size, bandwidth, levels)
+$$ LANGUAGE plpythonu;
+
+
+CREATE OR REPLACE FUNCTION
+  CDB_Contours (
+    subquery TEXT,
+    grid_size NUMERIC DEFAULT 100,
+    bandwidth NUMERIC DEFAULT 0.0001,
+    levels NUMERIC[] DEFAULT null
+    )
+RETURNS table (level Numeric, geom geometry )
+AS $$
+BEGIN
+
+  RETURN QUERY
+    select cont.level as level, ST_GeomFromText(cont.geom_text, 4326)::geometry as geom from _CDB_Contours(subquery,grid_size,bandwidth,levels) as cont;
+END;
+$$ LANGUAGE plpgsql;
+

src/pg/sql/80_similarity_rank.sql

@@ -1,15 +0,0 @@
-CREATE OR REPLACE FUNCTION cdb_SimilarityRank(cartodb_id numeric, query text)
-returns TABLE (cartodb_id NUMERIC, similarity NUMERIC)
-as $$
-  plpy.execute('SELECT cdb_crankshaft._cdb_crankshaft_activate_py()')
-  from crankshaft.similarity import similarity_rank
-  return similarity_rank(cartodb_id, query)
-$$ LANGUAGE plpythonu;
-
-CREATE OR REPLACE FUNCTION cdb_MostSimilar(cartodb_id numeric, query text ,matches numeric)
-returns TABLE (cartodb_id NUMERIC, similarity NUMERIC)
-as $$
-  plpy.execute('SELECT cdb_crankshaft._cdb_crankshaft_activate_py()')
-  from crankshaft.similarity import most_similar
-  return most_similar(matches, query)
-$$ LANGUAGE plpythonu;

src/py/crankshaft/crankshaft/__init__.py

@@ -1,3 +1,3 @@
 import random_seeds
 import clustering
-import similarity
+import contours

src/py/crankshaft/crankshaft/contours/__init__.py

@@ -0,0 +1 @@
+from contours import *

src/py/crankshaft/crankshaft/contours/contours.py

@@ -0,0 +1,58 @@
+from scipy.stats import gaussian_kde
+from scipy.interpolate import griddata
+import numpy as np 
+from sklearn.neighbors import KernelDensity
+from skimage.measure import find_contours
+import plpy
+
+def cdb_generate_contours(query, grid_size, bandwidth, levels):
+    plpy.notice('one')
+    data   = plpy.execute( 'select ST_X(the_geom) as x , ST_Y(the_geom) as y from ({0}) as a '.format(query))
+    plpy.notice('two')
+
+    xs = [d['x'] for d in data]
+    ys = [d['y'] for d in data]
+    plpy.notice('three')
+    return generate_contours(xs,ys,grid_size,bandwidth,levels)
+  
+def scale_coord(coord, x_range,y_range,grid_size):
+    plpy.notice('ranges %,  % ', x_range, y_range)
+    return [coord[0]*(x_range[1]-x_range[0])/float(grid_size)+x_range[0],
+            coord[1]*(y_range[1]-y_range[0])/float(grid_size)+y_range[0]]
+    
+def make_wkt(data,x_range, y_range, grid_size):
+    joined = ','.join([' '.join(map(str,scale_coord(coord_pair, x_range, y_range, grid_size))) for coord_pair in data])
+    return '({0})'.format(joined)
+    
+def make_multi_line(data,x_range,y_range, grid_size):
+    joined = ','.join([ make_wkt(ring,x_range,y_range,grid_size)  for ring in data ])
+    return 'MULTILINESTRING({0})'.format(joined)
+
+def generate_contours(xs,ys, grid_res=100, bandwidth=0.001, levels=None):
+    plpy.notice("HERE")
+    max_y, min_y = np.max(ys), np.min(ys)
+    max_x, min_x = np.max(xs), np.min(xs)
+    positions = np.vstack([ys,xs]).T
+    grid_x,grid_y = np.meshgrid(np.linspace(min_x, max_x , grid_res), np.linspace(min_y, max_y, grid_res))
+    xy = np.vstack([grid_y.ravel(), grid_x.ravel()]).T
+    xy *= np.pi / 180.
+
+    plpy.notice(" Generating kernel density")
+    kde = KernelDensity(bandwidth=bandwidth, metric='haversine',
+                        kernel='gaussian', algorithm='ball_tree')
+    kde.fit(positions*np.pi/180.)
+    results = np.exp(kde.score_samples(xy))
+    results = results.reshape((grid_x.shape[0], grid_y.shape[0]))
+    
+    if not levels:
+        levels = np.linspace(results.min(), results.max(),60)
+    plpy.notice(' finding contours')
+    CS = [find_contours(results, level) for level in levels]
+    
+    vertices = []
+    for contours,level in zip(CS,levels):
+        if len(contours)>0:
+            multiline = make_multi_line(contours, (min_x,max_x), (min_y, max_y), grid_res)
+            vertices.append([level, multiline ])
+    plpy.notice('generated vertices retunring ?')
+    return vertices

src/py/crankshaft/crankshaft/similarity/__init__.py

@@ -1 +0,0 @@
-from similarity import * 

src/py/crankshaft/crankshaft/similarity/similarity.py

@@ -1,91 +0,0 @@
-from sklearn.neighbors import NearestNeighbors
-import  scipy.stats as stats
-import numpy as np
-import plpy
-import time
-import cPickle
-
-
-def query_to_dictionary(result):
-    return [ dict(zip(r.keys(), r.values())) for r in result ]
-
-def drop_all_nan_columns(data):
-    return data[ :, ~np.isnan(data).all(axis=0)]
-    
-def fill_missing_na(data,val=None):
-    inds = np.where(np.isnan(data))
-    if val==None:
-        col_mean = stats.nanmean(data,axis=0)
-        data[inds]=np.take(col_mean,inds[1])
-    else:
-        data[inds]=np.take(val, inds[1])
-    return data
-    
-def similarity_rank(target_cartodb_id, query):
-    start_time  = time.time() 
-    #plpy.notice('converting to dictionary ', start_time) 
-    #data = query_to_dictionary(plpy.execute(query))  
-    plpy.notice('coverted , running query ', time.time() - start_time) 
-    
-    data = plpy.execute(query_only_values(query))
-    plpy.notice('run query  , getting cartodb_idsi', time.time() - start_time)
-    cartodb_ids = plpy.execute(query_cartodb_id(query))[0]['a']
-    target_id  = cartodb_ids.index(target_cartodb_id)
-    plpy.notice('run query  , extracting ', time.time() - start_time)
-    features, target = extract_features_target(data,target_id)
-    plpy.notice('extracted  , cleaning ', time.time() - start_time)
-    features = fill_missing_na(drop_all_nan_columns(features))
-    plpy.notice('cleaned , normalizing', start_time - time.time())
-    
-    normed_features, normed_target  = normalize_features(features,target)
-    plpy.notice('normalized , training ', time.time() - start_time )
-    tree = train(normed_features)
-    plpy.notice('normalized , pickling ', time.time() - start_time )
-    #plpy.notice('tree_dump ',  len(cPickle.dumps(tree, protocol=cPickle.HIGHEST_PROTOCOL)))
-    plpy.notice('pickles, querying ', time.time() - start_time)
-    dist, ind  = tree.kneighbors(normed_target)
-    plpy.notice('queried , rectifying', time.time() - start_time)
-    return zip(cartodb_ids, dist[0])
-
-def query_cartodb_id(query):
-    return 'select array_agg(cartodb_id) a from ({0}) b'.format(query)
-
-def query_only_values(query):
-    first_row = plpy.execute('select * from ({query}) a limit 1'.format(query=query))
-    just_values = ','.join([ key for key in  first_row[0].keys()  if key not in ['the_geom', 'the_geom_webmercator','cartodb_id']])
-    return 'select Array[{0}] a from ({1}) b '.format(just_values, query)
-
-
-def most_similar(matches,query):
-    data = plpy.execute(query)    
-    features, _ = extract_features_target(data)
-    results = []
-    for i in features:
-        target = features
-        dist,ind = tree.query(target, k=matches)
-        cartodb_ids  = [ dist[ind]['cartodb_id'] for index in ind ]
-        results.append(cartodb_ids)
-    return cartodb_ids, results
-    
-    
-def train(features):
-    tree = NearestNeighbors( n_neighbors=len(features), algorithm='auto').fit(features)
-    return tree
-    
-def normalize_features(features, target):
-    maxes = features.max(axis=0)
-    mins  = features.min(axis=0)
-    return (features - mins)/(maxes-mins), (target-mins)/(maxes-mins)
- 
-def extract_row(row):
-    keys = row.keys()
-    values = row.values()
-    del values[ keys.index('cartodb_id')]
-    return values
-
-def extract_features_target(data, target_index=None):
-    target   = None
-    features = [row['a'] for row in data]
-    target   = features[target_index]
-    return np.array(features, dtype=float), np.array(target, dtype=float)
-    

src/py/crankshaft/setup.py

@@ -40,9 +40,9 @@ setup(
 
     # The choice of component versions is dictated by what's
     # provisioned in the production servers.
-    install_requires=['pysal==1.9.1', 'scikit-learn==0.17.1'],
+    install_requires=['pysal==1.9.1'],
 
-    requires=['pysal', 'numpy','sklearn'],
+    requires=['pysal', 'numpy', 'sklearn', 'scikit-image'],
 
     test_suite='test'
 )