def contiguity_from_shapefile(shapefile, criteria='rook'):
    """
    Create a spatial weights object based on a contiguity criteria from a
    shapefile.

    Produces a "*.gal" file in the directory of the shapefile


    Argument
    --------
    shapefile: string with full path to shapefile

    criteria: string for type of contiguity ['rook'|'queen']

    Returns
    -------
    cards: nx1 numpy array with the number of neighbors for each element based
    on criterion

    """

    if criteria == 'rook':
        PS.rook_from_shapefile(shapefile)
        abb = 'r'
    else:
        PS.queen_from_shapefile(shapefile)
        abb = 'q'
    cards = NP.array(w.cardinalities.values())
    cards.shape = (len(cards),1)
    galfile = shapefile.split(".")[0] + "_" + abb + ".gal"
    gal = PS.open(galfile,'w')
    gal.write(w)
    gal.close()

    return cards

    def processAlgorithm(self, progress):
        field = self.getParameterValue(self.FIELD)
        field = field[0:10] # try to handle Shapefile field length limit
        filename = self.getParameterValue(self.INPUT)
        layer = dataobjects.getObjectFromUri(filename)
        filename = dataobjects.exportVectorLayer(layer)        
        
        contiguity = self.getParameterValue(self.CONTIGUITY)
        if contiguity == 0: # queen
            print 'INFO: Moran\'s using queen contiguity'
            w=pysal.queen_from_shapefile(filename)
        else: # 1 for rook
            print 'INFO: Moran\'s using rook contiguity'
            w=pysal.rook_from_shapefile(filename)
    
        f = pysal.open(filename.replace('.shp','.dbf'))
        y=np.array(f.by_col[str(field)])
        m = pysal.Moran(y,w,transformation = "r", permutations = 999)

        self.setOutputValue(self.I,m.I)
        
        print "Moran's I: %f" % (m.I)
        print "INFO: Moran's I values range from -1 (indicating perfect dispersion) to +1 (perfect correlation). Values close to -1/(n-1) indicate a random spatial pattern."
        print "p_norm: %f" % (m.p_norm)
        print "p_rand: %f" % (m.p_rand)
        print "p_sim: %f" % (m.p_sim)
        print "INFO: p values smaller than 0.05 indicate spatial autocorrelation that is significant at the 5% level."
        print "z_norm: %f" % (m.z_norm)
        print "z_rand: %f" % (m.z_rand)
        print "z_sim: %f" % (m.z_sim)
        print "INFO: z values greater than 1.96 or smaller than -1.96 indicate spatial autocorrelation that is significant at the 5% level."

 def setUp(self):
     db=pysal.open(pysal.examples.get_path("columbus.dbf"),"r")
     y = np.array(db.by_col("CRIME"))
     self.y = np.reshape(y, (49,1))
     X = []
     X.append(db.by_col("HOVAL"))
     X.append(db.by_col("INC"))
     self.X = np.array(X).T
     X2 = []
     X2.append(db.by_col("INC"))
     self.X2 = np.array(X2).T
     yd = []
     yd.append(db.by_col("HOVAL"))
     self.yd = np.array(yd).T
     q = []
     q.append(db.by_col("DISCBD"))
     self.q = np.array(q).T
     self.w = pysal.queen_from_shapefile(pysal.examples.get_path("columbus.shp"))
     self.w.transform = 'r'
     self.r_var = 'NSA'
     self.regimes = db.by_col(self.r_var)
     #Artficial:
     n = 256
     self.n2 = n/2
     self.x_a1 = np.random.uniform(-10,10,(n,1))
     self.x_a2 = np.random.uniform(1,5,(n,1))
     self.q_a = self.x_a2 + np.random.normal(0,1,(n,1))
     self.x_a = np.hstack((self.x_a1,self.x_a2))
     self.y_a = np.dot(np.hstack((np.ones((n,1)),self.x_a)),np.array([[1],[0.5],[2]])) + np.random.normal(0,1,(n,1))
     latt = int(np.sqrt(n))
     self.w_a = pysal.lat2W(latt,latt)
     self.w_a.transform='r'
     self.regi_a = [0]*(n/2) + [1]*(n/2)
     self.w_a1 = pysal.lat2W(latt/2,latt)
     self.w_a1.transform='r'

 def get(self,shpName='',width=0,height=0):
     print shpName,width,height
     if not shpName in self.SHPS:
         return self.index()
     shp = pysal.open(self.SHPS[shpName])
     W = None
     if 'w' in self.request.GET:
         wtype = self.request.GET['w']
         if wtype.lower() == 'rook':
             W = pysal.rook_from_shapefile(self.SHPS[shpName])
         elif wtype.lower() == 'queen':
             W = pysal.queen_from_shapefile(self.SHPS[shpName])
         else:
             try:
                 k = int(wtype)
                 W = pysal.knnW_from_shapefile(self.SHPS[shpName],k)
             except:
                 print "No valid W"
     print shp
     if width and height:
         width=int(width)
         height=int(height)
         if W:
             return self.write({'len':len(shp), 'polygons':shift_scale_shp(shp,width,height),'width':width,'height':height,'W':W.neighbors})
         else:
             return self.write({'len':len(shp), 'polygons':shift_scale_shp(shp,width,height),'width':width,'height':height,'W':'null'})
     return self.write({'len':len(shp)})

def read_files(filepath, **kwargs):
    """
    Reads a dbf/shapefile pair, squashing geometries into a "geometry" column.
    """
    #keyword arguments wrapper will strip all around dbf2df's required arguments
    geomcol = kwargs.pop('geomcol', 'geometry')
    weights = kwargs.pop('weights', '')
    
    dbf_path, shp_path = _pairpath(filepath)

    df = dbf2df(dbf_path, **kwargs)
    df[geomcol] = shp2series(shp_path)

    if weights != '' and isinstance(weights, str):
        if weights.lower() in ['rook', 'queen']:
            if weights.lower() == 'rook':
                df.W = ps.rook_from_shapefile(shp_path)
            else:
                df.W = ps.queen_from_shapefile(shp_path)
        else:
            try:
                W_path = os.path.splitext(dbf_path)[0] + '.' + weights
                df.W = ps.open(W_path).read()
            except IOError:
                print('Weights construction failed! Passing on weights')
    
    return df

    def contiguity_from_shapefile(shapefile, criteria='rook'):
        print shapefile
        if criteria == 'rook':
	    PS.rook_from_shapefile(shapefile)
            abb = 'r'
    	else:
            PS.queen_from_shapefile(shapefile)
            abb = 'q'
    	    cards = NP.array(w.cardinalities.values())
    	    cards.shape = (len(cards),1)
    	    galfile = shapefile.split(".")[0] + "_" + abb + ".gal"
    	    gal = PS.open(galfile,'w')
    	    gal.write(w)
    	    gal.close()

    	return cards

 def test_names(self):
     w = pysal.queen_from_shapefile(pysal.examples.get_path('columbus.shp'))
     gwk = pysal.kernelW_from_shapefile(pysal.examples.get_path('columbus.shp'),k=5,function='triangular', fixed=False)
     name_x = ['inc']
     name_y = 'crime'
     name_yend = ['hoval']
     name_q = ['discbd']
     name_w = 'queen'
     name_gwk = 'k=5'
     name_ds = 'columbus'
     reg = TSLS(self.y, self.X, self.yd, self.q,
             spat_diag=True, w=w, robust='hac', gwk=gwk,
             name_x=name_x, name_y=name_y, name_q=name_q, name_w=name_w,
             name_yend=name_yend, name_gwk=name_gwk, name_ds=name_ds)
     betas = np.array([[ 88.46579584], [  0.5200379 ], [ -1.58216593]])
     np.testing.assert_array_almost_equal(reg.betas, betas, 7)
     vm = np.array([[ 225.0795089 ,   17.11660041,  -12.22448566],
                    [  17.67097154,    2.47483461,   -1.4183641 ],
                    [ -12.45093722,   -1.40495464,    0.8700441 ]])
     np.testing.assert_array_almost_equal(reg.vm, vm, 7)
     self.assertListEqual(reg.name_x, ['CONSTANT']+name_x)
     self.assertListEqual(reg.name_yend, name_yend)
     self.assertListEqual(reg.name_q, name_q)
     self.assertEqual(reg.name_y, name_y)
     self.assertEqual(reg.name_w, name_w)
     self.assertEqual(reg.name_gwk, name_gwk)
     self.assertEqual(reg.name_ds, name_ds)

    def compute(self, vlayer, tfield, idvar,matType):
	vlayer=qgis.utils.iface.activeLayer()
	idvar=self.idVariable.currentText()
        # print type(idvar)
	tfield=self.inField.currentText()
        # print type(tfield)
	provider=vlayer.dataProvider()
	allAttrs=provider.attributeIndexes()
	caps=vlayer.dataProvider().capabilities()
	if caps & QgsVectorDataProvider.AddAttributes:
            TestField = idvar[:5]+"_qrr"
	    res = vlayer.dataProvider().addAttributes([QgsField(TestField, QVariant.Double)])
	wp=str(self.dic[str(self.inShape.currentText())])
        if matType == "Rook":
           w = py.rook_from_shapefile(wp, idVariable=unicode(idvar))
        else:
           w = py.queen_from_shapefile(wp, idVariable=unicode(idvar))

	w1=wp[:-3]+"dbf"
	db=py.open(w1)
	y=np.array(db.by_col[unicode(tfield)])
      	np.random.seed(12345)
	gc = py.Geary(y,w)
        #lm=py.Moran_Local(y,w)
	#l=lm.p_sim
	gg = gc.C
        self.SAresult.setText("The Global Geary's C index is " + str(gg))

def test2():
    #d6cd52286e5d3e9e08a5a42489180df3.shp
    import pysal
    shp = pysal.open('../www/tmp/d6cd52286e5d3e9e08a5a42489180df3.shp')
    dbf = pysal.open('../www/tmp/d6cd52286e5d3e9e08a5a42489180df3.dbf')        
    
    w = pysal.queen_from_shapefile('../www/tmp/d6cd52286e5d3e9e08a5a42489180df3.shp')
    #d3viz.moran_scatter_plot(shp, dbf, "dog_cnt", w)    
    import numpy as np
    y = np.array(dbf.by_col["dog_cnt"])
    lm = pysal.Moran_Local(y, w)    
    
    for i,j in enumerate(lm.q):
        if lm.p_sim[i] >= 0.05:
            print 0
        else:
            print j
        
    import d3viz
    d3viz.setup()
    
    d3viz.show_map(shp)
    d3viz.scatter_plot(shp, ["dog_cnt","home_cnt"])
    
    d3viz.lisa_map(shp, "dog_cnt", lm)

def test_d3viz():
    import pysal
    shp = pysal.open(pysal.examples.get_path('NAT.shp'),'r')
    dbf = pysal.open(pysal.examples.get_path('NAT.dbf'),'r')    

    #shp = pysal.open('/Users/xun/github/PySAL-Viz/test_data/man_road.shp','r')
    #dbf = pysal.open('/Users/xun/github/PySAL-Viz/test_data/man_road.dbf','r')
    import d3viz
    d3viz.setup()
    
    d3viz.init_map(shp)
    d3viz.show_map(shp)
   
    d3viz.scatter_plot(shp, field_x="HR90", field_y="HC90") 
    w = pysal.queen_from_shapefile(pysal.examples.get_path('NAT.shp'))
    d3viz.moran_scatter_plot(shp, dbf, "HR90", w)
    
    d3viz.quantile_map(shp, "HR90", 5)
    d3viz.quantile_map(shp, "HR90", 5, basemap="leaflet")
   
    import numpy as np
    y = np.array(dbf.by_col["HR90"])
    lm = pysal.Moran_Local(y, w)
    
    d3viz.lisa_map(shp, "HR90", lm)

    def __init__(self, filepath, outname, namelist, idlist, nb="queen", factor=2):
        """
        Initiation of modules
        """
        f=Dbf(filepath+".dbf")
        #Create mapping of locations to row id
        self.locations = dict()
        i=0
        for row in f:
            uid=unicode("".join([row[k] for k in idlist]))
            locnames = unicode(", ".join([row[k] for k in namelist]),"ascii","ignore")
            self.locations[i] = {outname:locnames,"id":uid}
            i+=1
        self.__dict__[outname]= self.locations
        self.outname = outname
        #Get Neightbor weights by queen, rook, knn, distance
        if nb=="queen":
            self.wt = pysal.queen_from_shapefile(filepath+".shp")
        elif nb=="rook":
            self.wt = pysal.rook_from_shapefile(filepath+".shp")
        elif nb=="knn":
            self.wt = pysal.knnW_from_shapefile(filepath+".shp", k=factor)
        elif nb=="distance":
            self.wt = pysal.threshold_binaryW_from_shapefile(filepath+".shp",k)

        #Create dictionary of neighbors for each region
        self.neighbors ={}
        for i,j in enumerate(self.wt):
            self.neighbors[self.locations[i]["id"]] = {self.outname:self.locations[i][self.outname]
                            ,"neighbors":dict([[self.locations[k]["id"],self.locations[k][self.outname]] for k in j.keys()])}

    def accept(self):
        # look for open shapefile layers, if none 
        if len(self.comboBox.currentText()) == 0 and self.lineEdit.text() == "":
            QMessageBox.information(self, self.tr("Weights from Shapefile"), self.tr("Please select input polygon vector layer"))

        # elif self.outShape.text() == "":
        #    QMessageBox.information(self, self.tr("Sum Line Lengths In Polyons"), self.tr("Please specify output shapefile"))
        else:
            # run the PySAL logic
            if str(self.comboBox.currentText()) == "":
                shapefile = str(self.shapefile)
            else:
                shapefile = str(self.d[str(self.comboBox.currentText())])
                       
                if self.radioButton.isChecked():
                    w = PS.queen_from_shapefile(shapefile)
                    abb = 'q'
                else:
                    w = PS.rook_from_shapefile(shapefile)
                    abb = 'r'
    	        cards = NP.array(w.cardinalities.values())
    	        cards.shape = (len(cards),1)
    	        galfile = shapefile.split(".")[0] + "_" + abb + ".gal"
    	        gal = PS.open(galfile,'w')
    	        gal.write(w)
    	        gal.close()
            QDialog.accept(self)

 def setUp(self):
     self.w = pysal.queen_from_shapefile(pysal.examples.get_path("columbus.shp"))
     self.w.transform = "r"
     self.db = pysal.open(pysal.examples.get_path("columbus.dbf"), "r")
     y = np.array(self.db.by_col("CRIME"))
     self.y = np.reshape(y, (49, 1))
     self.r_var = "NSA"
     self.regimes = self.db.by_col(self.r_var)

    def get_weight_matrix(self, array, rook=False, shpfile=None):
        """Return the spatial weight matrix based on pysal functionalities

        Keyword arguments:
            array    Numpy array with inventory values.
            rook    Boolean to select spatial weights matrix as rook or
                    queen case.
            shpfile    Name of file used to setup weight matrix.
        """
        # Get case name.
        if rook:
            case = 'rook'
        else:
            case = 'queen'
        # Get grid dimension.
        dim = array.shape
        if self.sptype == 'vector':
            try:
                # Create weights based on shapefile topology using defined key.
                if shpfile is None:
                    shpfile = self.invfile
                # Differentiat between rook and queen's case.
                if rook:
                    w = pysal.rook_from_shapefile(shpfile, self.invcol)
                else:
                    w = pysal.queen_from_shapefile(shpfile, self.invcol)
            except:
                msg = "Couldn't build spatial weight matrix for vector "
                "inventory <%s>" % (self.name)
                raise RuntimeError(msg)

            # Match weight index to inventory array index.
            w.id_order = list(self.inv_index)

            logger.info("Weight matrix in %s's case successfully calculated "
                        "for vector dataset" % case)
        elif self.sptype == 'raster':
            try:
                # Construct weight matrix in input grid size.
                w = pysal.lat2W(*dim, rook=rook)
            except:
                msg = "Couldn't build spatial weight matrix for raster "
                "inventory <%s>" % (self.name)
                raise RuntimeError(msg)

            logger.info("Weight matrix in %s's case successfully calculated "
                        "for raster dataset" % case)

        # Print imported raster summary.
        print("[ WEIGHT NUMBER ] = ", w.n)
        print("[ MIN NEIGHBOR ] = ", w.min_neighbors)
        print("[ MAX NEIGHBOR ] = ", w.max_neighbors)
        print("[ ISLANDS ] = ", *w.islands)
        print("[ HISTOGRAM ] = ", *w.histogram)

        self._Inventory__modmtime()

        return(w)

def calc_shp_adjacencies(shp, sort_id=None):

    if sort_id:
        adj=ps.queen_from_shapefile(shp,sort_id)
    else:
        adj=ps.queen_from_shapefile(shp)
    sort_order=sorted(adj.id_order)
    
    neigh_dict=adj.neighbors
    
    neigh_list_nest=[neigh_dict[i] for i in sort_order]
    neigh_list_flat=[k for sub in neigh_list_nest for k in sorted(sub)]
    
    neigh_id_list=[sort_order.index(X)+1 for X in neigh_list_flat]

    num_neigh_list=[len(x) for x in neigh_list_nest]
    
    return {'num_neigh':num_neigh_list, 'neigh_id':neigh_id_list,'shp':shp}

 def test_names(self):
     X = np.array(self.db.by_col("INC"))
     X = np.reshape(X, (49, 1))
     X = SP.csr_matrix(X)
     yd = np.array(self.db.by_col("CRIME"))
     yd = np.reshape(yd, (49, 1))
     q = np.array(self.db.by_col("DISCBD"))
     q = np.reshape(q, (49, 1))
     w = pysal.queen_from_shapefile(pysal.examples.get_path("columbus.shp"))
     gwk = pysal.kernelW_from_shapefile(
         pysal.examples.get_path("columbus.shp"), k=5, function="triangular", fixed=False
     )
     name_x = ["inc"]
     name_y = "crime"
     name_yend = ["crime"]
     name_q = ["discbd"]
     name_w = "queen"
     name_gwk = "k=5"
     name_ds = "columbus"
     reg = GM_Lag(
         self.y,
         X,
         yd,
         q,
         spat_diag=True,
         w=w,
         robust="hac",
         gwk=gwk,
         name_x=name_x,
         name_y=name_y,
         name_q=name_q,
         name_w=name_w,
         name_yend=name_yend,
         name_gwk=name_gwk,
         name_ds=name_ds,
     )
     betas = np.array([[5.46344924e01], [4.13301682e-01], [-5.92637442e-01], [-7.40490883e-03]])
     np.testing.assert_array_almost_equal(reg.betas, betas, 7)
     vm = np.array(
         [
             [5.70817052e02, -1.83655385e01, -8.36602575e00, 2.37538877e-02],
             [-1.85224661e01, 6.53311383e-01, 2.84209566e-01, -6.47694160e-03],
             [-8.31105622e00, 2.78772694e-01, 1.38144928e-01, -3.98175246e-03],
             [2.66662466e-02, -6.23783104e-03, -4.11092891e-03, 1.10936528e-03],
         ]
     )
     np.testing.assert_array_almost_equal(reg.vm, vm, 6)
     self.assertListEqual(reg.name_x, ["CONSTANT"] + name_x)
     name_yend.append("W_crime")
     self.assertListEqual(reg.name_yend, name_yend)
     name_q.extend(["W_inc", "W_discbd"])
     self.assertListEqual(reg.name_q, name_q)
     self.assertEqual(reg.name_y, name_y)
     self.assertEqual(reg.name_w, name_w)
     self.assertEqual(reg.name_gwk, name_gwk)
     self.assertEqual(reg.name_ds, name_ds)

    def estimate_elasticity(self, zones):
        dummies = pd.get_dummies(zones.county)
        zones = pd.concat([zones, dummies], axis=1)
        zones['avg_far'] = self.buildings_far.groupby('zone_id').far.mean() #use far_x because Xavier's code adds far to buildings

        #zones = zones[zones.residential_sqft_zone>0]

        #wrook = py.queen_from_shapefile('C:/users/jmartinez/documents/Test Zones/zones_prj_res2.shp')
        wqueen = py.queen_from_shapefile(os.path.join(misc.data_dir(),'shapefiles\\zones.shp'))
        w = py.weights.weights.W(wqueen.neighbors, wqueen.weights)
        x = zones[['zonal_pop','mean_income']]
        x = x.apply(np.log1p)

        x['ln_jobs_within_30min'] = zones['ln_jobs_within_30min']
        x['zone_contains_park'] = zones['zone_contains_park']
        x['Arapahoe'] = zones['Arapahoe']
        x['Boulder'] = zones['Boulder']
        x['Broomfield'] = zones['Broomfield']
        x['Clear Creek'] = zones['Clear Creek']
        x['Denver'] = zones['Denver']
        x['Douglas'] = zones['Douglas']
        x['Elbert'] = zones['Elbert']
        x['Gilpin'] = zones['Gilpin']
        x['Jefferson'] = zones['Jefferson']
        x['Weld'] = zones['Weld']
        x=x.fillna(0)
        x = x.as_matrix()

        imat = zones[['ln_avg_nonres_unit_price_zone','avg_far']]
        imat = imat.fillna(0)
        imat = imat.as_matrix()

        yend = zones['ln_avg_unit_price_zone']
        yend = yend.fillna(0)
        yend = yend.as_matrix()
        yend = np.reshape(yend,(zones.shape[0],1))

        y = zones['residential_sqft_zone']
        y = y.fillna(0)
        y = y.apply(np.log1p)
        y = y.as_matrix()
        y = np.reshape(y,(zones.shape[0],1))


        imat_names = ['non_res_price','avg_far']
        x_names = ['zonal_pop', 'mean_income', 'ln_jobs_within_30min', 'zone_contains_park','Arapahoe','Boulder','Broomfield','Clear Creek','Denver','Douglas','Elbert','Gilpin','Jefferson','Weld']
        yend_name = ['ln_avg_unit_price_zone']
        y_name = 'residential_sqft_zone'

        reg_2sls = py.spreg.twosls_sp.GM_Lag(y, x, yend=yend, q=imat, w=w, w_lags=2, robust ='white', name_x = x_names, name_q = imat_names, name_y = y_name, name_yend = yend_name)

        demand_elasticity = np.absolute(reg_2sls.betas[15])
        demand_elasticity = 1/demand_elasticity[0]
        #
        return demand_elasticity

 def test_spatial(self):
     w = pysal.queen_from_shapefile(pysal.examples.get_path('columbus.shp'))
     reg = TSLS(self.y, self.X, self.yd, self.q, spat_diag=True, w=w)
     betas = np.array([[ 88.46579584], [  0.5200379 ], [ -1.58216593]])
     np.testing.assert_array_almost_equal(reg.betas, betas, 7)
     vm = np.array([[ 229.05640809,   10.36945783,   -9.54463414],
                    [  10.36945783,    2.0013142 ,   -1.01826408],
                    [  -9.54463414,   -1.01826408,    0.62914915]]) 
     np.testing.assert_array_almost_equal(reg.vm, vm, 7)
     ak_test = np.array([ 1.16816972,  0.27977763])
     np.testing.assert_array_almost_equal(reg.ak_test, ak_test, 7)

def main(perms):
    np.random.seed(10)
    w = pysal.queen_from_shapefile(pysal.examples.get_path("NAT.shp"))
    f = pysal.open(pysal.examples.get_path('NAT.dbf'))
    y = np.array(f.by_col('POL90'))
    y = np.reshape(y, (3085,))
    t1 = time.time()
    lm = Moran_Local(y, w, transformation = "r", permutations = perms, cores=None)
    t2 = time.time()
    print "Local Morans (only) time was {}".format(t2-t1)
    print "I value: {}".format(lm.Is)
    print "p sum: {}".format(lm.p_z_sim)

def test_mplot():
    link = ps.examples.get_path('columbus.shp')

    db = read_files(link)
    y = db['HOVAL'].values
    w = ps.queen_from_shapefile(link)
    w.transform = 'R'

    m = ps.Moran(y, w)

    fig = mplot(m, xlabel='Response', ylabel='Spatial Lag',
                title='Moran Scatterplot', custom=(7,7))
    plt.close(fig)