def inverse_scg(w, data, scalar, transpose=False, symmetric=False,\
               threshold=0.00001,\
               max_iterations=None):
    
    #multiplier = SP.identity(w.n) - (scalar*w.sparse)       # A      n x n   
    count = 0                                               # k      scalar (step 1)
    run_tot = copy.copy(data)                               # z_k    n x 1  (step 1)
    #residuals = data - run_tot * multiplier                 # r_k    n x 1  (step 2)
    residuals = data - pysal.lag_spatial(w, scalar*data)
    #test1 = la.norm(residuals)                              # G_k    scalar (step 3)
    test1 = norm(residuals)
    directions = copy.copy(residuals)                       # d_k    n x 1  (step 6)
    while test1 > threshold:                                #               (step 4)
        count += 1                                          #               (step 5)
        #changes = multiplier * directions                   # t      n x 1  (step 7)
        changes = directions - pysal.lag_spatial(w, scalar*directions)
        intensity = test1 / (np.dot(directions.T, changes)) # alpha  scalar (step 8)
        #int_dir = intensity * directions                    #               (step 8)
        run_tot += intensity * directions
        #run_tot += int_dir                                  #               (step 8)
        #residuals -= int_dir                                #               (step 8)
        residuals -= intensity * changes
        #test2 = la.norm(residuals)                          #               (step 3)
        test2 = norm(residuals)
        directions = residuals + ((test2/test1)*directions) #               (step 6)
        test1 = test2
    return run_tot

def set_endog(y, x, w, yend, q, w_lags, lag_q):
    # Create spatial lag of y
    yl = lag_spatial(w, y)
    # spatial and non-spatial instruments
    if issubclass(type(yend), np.ndarray):
        if lag_q:
            lag_vars = sphstack(x, q)
        else:
            lag_vars = x
        spatial_inst = get_lags(w, lag_vars, w_lags)
        q = sphstack(q, spatial_inst)
        yend = sphstack(yend, yl)
    elif yend == None:  # spatial instruments only
        q = get_lags(w, x, w_lags)
        yend = yl
    else:
        raise Exception, "invalid value passed to yend"
    return yend, q

    lag = lag_spatial(w, x)
    spat_lags = lag
    for i in range(w_lags - 1):
        lag = lag_spatial(w, lag)
        spat_lags = sphstack(spat_lags, lag)
    return spat_lags

def get_lags(w, x, w_lags):
    '''
    Calculates a given order of spatial lags and all the smaller orders

    Parameters
    ----------
    w       : weight
              PySAL weights instance
    x       : array
              nxk arrays with the variables to be lagged  
    w_lags  : integer
              Maximum order of spatial lag

    Returns
    --------
    rs      : array
              nxk*(w_lags+1) array with original and spatially lagged variables

    '''
    lag = lag_spatial(w, x)
    spat_lags = lag
    for i in range(w_lags-1):
        lag = lag_spatial(w, lag)
        spat_lags = sphstack(spat_lags, lag)
    return spat_lags

 def __calc(self, z):
     zl = pysal.lag_spatial(self.w, z)
     bb = sum(z * zl) / 2.0
     zw = 1 - z
     zl = pysal.lag_spatial(self.w, zw)
     ww = sum(zw * zl) / 2.0
     bw = self.J - (bb + ww)
     return (bb, ww, bw)

 def get_x_lag(self, w, regimes_att):
     if regimes_att:
         xlag = ps.lag_spatial(w, regimes_att['x'])
         xlag = REGI.Regimes_Frame.__init__(self, xlag,
                                            regimes_att['regimes'], constant_regi=None, cols2regi=regimes_att['cols2regi'])[0]
         xlag = xlag.toarray()
     else:
         xlag = ps.lag_spatial(w, self.x)
     return xlag

def moran_plot(IM):

    import pandas as pd
    import seaborn as sns
    import matplotlib.pyplot as plt
    import numpy as np
    import pysal as ps

    y_norm = normalize(IM.y)
    y_lag = ps.lag_spatial(IM.w, IM.y)
    y_lag_norm = normalize(y_lag)
    dados = pd.DataFrame({'y':IM.y, 'y_norm':y_norm,
                          'y_lag':y_lag, 'y_lag_norm':y_lag_norm})

    f, ax = plt.subplots(1, figsize=(7, 5))
    sns.regplot('y_norm', 'y_lag_norm', data=dados, ci=None,
                color='black', line_kws={'color':'red'})
    plt.axvline(0, c='gray', alpha=0.7)
    plt.axhline(0, c='gray', alpha=0.7)

    limits = np.array([y_norm.min(), y_norm.max(), y_lag_norm.min(), y_lag_norm.max()])
    limits = np.abs(limits).max()
    border = 0.02
    ax.set_xlim(- limits - border, limits + border)
    ax.set_ylim(- limits - border, limits + border)

    plt.show();

 def __calc(self, z, op):
     if op == 'c':     # cross-product
         zl = pysal.lag_spatial(self.w, z)
         g = (z * zl).sum()
     elif op == 's':   # squared difference
         zs = np.zeros(z.shape)
         z2 = z ** 2
         for i, i0 in enumerate(self.w.id_order):
             neighbors = self.w.neighbor_offsets[i0]
             wijs = self.w.weights[i0]
             zw = zip(neighbors, wijs)
             zs[i] = sum([wij * (z2[i] - 2.0 * z[i] * z[
                 j] + z2[j]) for j, wij in zw])
         g = zs.sum()
     elif op == 'a':    # absolute difference
         zs = np.zeros(z.shape)
         for i, i0 in enumerate(self.w.id_order):
             neighbors = self.w.neighbor_offsets[i0]
             wijs = self.w.weights[i0]
             zw = zip(neighbors, wijs)
             zs[i] = sum([wij * abs(z[i] - z[j]) for j, wij in zw])
         g = zs.sum()
     else:              # any previously defined function op
         zs = np.zeros(z.shape)
         for i, i0 in enumerate(self.w.id_order):
             neighbors = self.w.neighbor_offsets[i0]
             wijs = self.w.weights[i0]
             zw = zip(neighbors, wijs)
             zs[i] = sum([wij * op(z, i, j) for j, wij in zw])
         g = zs.sum()
     return g

 def test_lag_spatial(self):
     yl = pysal.lag_spatial(self.w, self.y)
     np.testing.assert_array_almost_equal(yl, [1., 2., 1.])
     self.w.id_order = ['b', 'c', 'a']
     y = np.array([1, 2, 0])
     yl = pysal.lag_spatial(self.w, y)
     np.testing.assert_array_almost_equal(yl, [2., 1., 1.])
     w = pysal.lat2W(3, 3)
     y = np.arange(9)
     yl = pysal.lag_spatial(w, y)
     ylc = np.array([4., 6., 6., 10., 16., 14., 10., 18., 12.])
     np.testing.assert_array_almost_equal(yl, ylc)
     w.transform = 'r'
     yl = pysal.lag_spatial(w, y)
     ylc = np.array(
         [2., 2., 3., 3.33333333, 4.,
          4.66666667, 5., 6., 6.])
     np.testing.assert_array_almost_equal(yl, ylc)

 def rose(self, Y, w, k=8):
     sw = 2*np.pi/k 
     cuts = np.arange(0.0,2*np.pi+sw,sw)
     wY = ps.lag_spatial(w,Y)
     dx = Y[:,-1]-Y[:,0]
     dy = wY[:,-1]-wY[:,0]
     theta = np.arctan2(dy,dx)
     neg = theta < 0.0
     utheta = theta*(1-neg) + neg * (2*np.pi+theta)         
     return cuts, utheta, dx, dy

def drawW(rel,w,k):
    r=np.random.permutation(rel)
    wy=pysal.lag_spatial(w,r)
    y=wy[:,-1]-wy[:,0]
    x=r[:,-1]-r[:,0]
    theta=np.arctan2(y,x)
    neg=theta < 0.0
    utheta=theta*(1-neg) + neg * (2*np.pi+theta)
    k=8
    width=2*np.pi/k
    cuts=np.arange(0.0,2*np.pi+width,width)
    counts,bin=np.histogram(utheta,cuts)
    return counts

def AddLagVars(dataframe, shp, idfield, SumVars, AvgVars, Adj=None):
    df=dataframe.copy()
    if idfield not in df.columns:
        print 'Dataframe missing id field'
    if idfield not in pysal.open(shp[:-3]+'dbf').header:
        print 'Shp missing id field'

    if Adj==None:
        Adj=pysal.queen_from_shapefile(shp, idVariable=idfield)
    else:
        pass
   
    df.set_index(idfield, inplace=True)
    df=df.reindex(Adj.id_order)
    
    Adj.transform='o'
    for Var in SumVars:
        df[Var+'_LAG_SUM']=pysal.lag_spatial(Adj, np.array(df[Var]))
    Adj.transform=('r')    
    for Var in AvgVars:
        df[Var+'_LAG_AVG']=pysal.lag_spatial(Adj, np.array(df[Var])) 

    df.reset_index(inplace=True) 
    return df

def log_lik_lag(ldet, w, b, X, y):
    n = w.n
    r = b[0]    # ml estimate of rho
    b = b[1:]   # ml for betas
    yl = ps.lag_spatial(w,y)
    ys = y - r * yl
    XX = np.dot(X.T, X)
    iXX = np.linalg.inv(XX)
    b = np.dot(iXX, np.dot(X.T,ys))
    yhat = r * yl + np.dot(X,b)
    e = y - yhat
    e2 = (e**2).sum()
    sig2 = e2 / n
    ln2pi = np.log(2*np.pi)
    return ldet - n/2. * ln2pi - n/2. * np.log(sig2) - e2/(2 * sig2)

def mplot(m, xlabel='', ylabel='', title='', custom=(7,7)):
    '''
    Produce basic Moran Plot 
    ...
    Parameters
    ---------
    m            : array
                   values of Moran's I 
    xlabel       : str
                   label for x axis
    ylabel       : str
                   label for y axis                
    title        : str
                   title of plot
    custom       : tuple
                   dimensions of figure size

    Returns
    ---------
    plot         : png 
                    image file showing plot
            
    '''
    
    lag = ps.lag_spatial(m.w, m.z)
    fit = ps.spreg.OLS(m.z[:, None], lag[:,None])

    ## Customize plot
    fig = plt.figure(figsize=custom)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.suptitle(title)

    plt.scatter(m.z, lag, s=60, color='k', alpha=.6)
    plt.plot(lag, fit.predy, color='r')

    plt.axvline(0, alpha=0.5)
    plt.axhline(0, alpha=0.5)
    plt.show()

    return None

def hac_multi(reg, gwk, constant=False):
    """
    HAC robust estimation of the variance-covariance matrix for multi-regression object 

    Parameters
    ----------

    reg             : Regression object (OLS or TSLS)
                      output instance from a regression model

    gwk             : PySAL weights object
                      Spatial weights based on kernel functions

    Returns
    --------

    psi             : kxk array
                      Robust estimation of the variance-covariance

    """
    if not constant:
        reg.hac_var = check_constant(reg.hac_var)
    xu = spbroadcast(reg.hac_var, reg.u)
    gwkxu = lag_spatial(gwk, xu)
    psi0 = spdot(xu.T, gwkxu)
    counter = 0
    for m in reg.multi:
        reg.multi[m].robust = 'hac'
        reg.multi[m].name_gwk = reg.name_gwk
        try:
            psi1 = spdot(reg.multi[m].varb, reg.multi[m].zthhthi)
            reg.multi[m].vm = spdot(psi1, np.dot(psi0, psi1.T))
        except:
            reg.multi[m].vm = spdot(
                reg.multi[m].xtxi, np.dot(psi0, reg.multi[m].xtxi))
        reg.vm[(counter * reg.kr):((counter + 1) * reg.kr),
               (counter * reg.kr):((counter + 1) * reg.kr)] = reg.multi[m].vm
        counter += 1

    def _calc(self, y, w, classes, k):
        # lag markov
        ly = pysal.lag_spatial(w, y)
        npm = np.matrix
        npa = np.array
        if self.fixed:
            l_classes = pysal.Quantiles(ly.flatten(), k=k).yb
            l_classes.shape = ly.shape
        else:
            l_classes = npa([pysal.Quantiles(
                ly[:, i], k=k).yb for i in np.arange(self.cols)])
            l_classes = l_classes.transpose()
        l_classic = Markov(l_classes)
        T = np.zeros((k, k, k))
        n, t = y.shape
        for t1 in range(t - 1):
            t2 = t1 + 1
            for i in range(n):
                T[l_classes[i, t1], classes[i, t1], classes[i, t2]] += 1

        P = np.zeros_like(T)
        F = np.zeros_like(T)  # fmpt
        ss = np.zeros_like(T[0])
        for i, mat in enumerate(T):
            row_sum = mat.sum(axis=1)
            row_sum = row_sum + (row_sum == 0)
            p_i = np.matrix(np.diag(1. / row_sum) * np.matrix(mat))
            #print i
            #print mat
            #print p_i
            ss[i] = steady_state(p_i).transpose()
            try:
                F[i] = fmpt(p_i)
            except:
                #pylint; "No exception type(s) specified"
                print "Singlular fmpt matrix for class ", i
            P[i] = p_i
        return T, P, ss, F

def moran_scatter_plot(shp, dbf, var, w):
    y = np.array(dbf.by_col[var])
    y_lag = pysal.lag_spatial(w, y)
   
    y_z = (y - y.mean()) / y.std()
    y_lag_z = (y_lag - y_lag.mean()) / y_lag.std()
    
    global SHP_DICT 
    uuid = SHP_DICT[shp]
    
    global WS_SERVER 
    ws = create_connection(WS_SERVER)
    msg = {
        "command": "moran_scatter_plot",
        "uuid":  uuid,
        "title": "Moran Scatter plot for variable [%s]" % var,
        "data": { "x": y_z.tolist(), "y" : y_lag_z.tolist() },
        "fields": [var, "lagged %s" % var]
    }
    str_msg = json.dumps(msg)
    ws.send(str_msg)
    #print "send:", str_msg
    ws.close()

def moran_dispersao(IM, title='', xlabel='', ylabel=''):

    y_norm = normalizar(IM.y)
    y_lag = ps.lag_spatial(IM.w, IM.y)
    y_lag_norm = normalizar(y_lag)
    dados = pd.DataFrame({'y':IM.y, 'y_norm':y_norm,
                          'y_lag':y_lag, 'y_lag_norm':y_lag_norm})

    f, ax = plt.subplots(1, figsize=(7, 5))
    sns.regplot('y_norm', 'y_lag_norm', data=dados, ci=None,
                color='black', line_kws={'color':'red'})
    plt.axvline(0, c='gray', alpha=0.7)
    plt.axhline(0, c='gray', alpha=0.7)

    limits = np.array([y_norm.min(), y_norm.max(), y_lag_norm.min(), y_lag_norm.max()])
    limits = np.abs(limits).max()
    border = 0.02
    ax.set_xlim(- limits - border, limits + border)
    ax.set_ylim(- limits - border, limits + border)
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.show();

    def grafico(self):

        """Grafico de dispersion"""
        w=self.w
        y=self.y
        ystd=(y-np.mean(y))/np.std(y)
        w.transform = 'r'
        yl = pysal.lag_spatial(w,ystd)

        colors = np.random.rand(100)
        """area = np.pi * (15 * np.random.rand(100))**2 # 0 to 15 point radiuses"""

        fig, ax = plt.subplots()
        m, fit = np.polyfit(ystd, yl, deg=1)
        """ax.plot(self.y, fit[0] * self.y + fit[1], color='blue', alpha=0.4, linewidth=0.3, linestyle='dotted')"""
        ax.scatter(ystd, yl, c=colors, alpha=0.5)
        ax.plot(ystd, m*ystd + fit, color='blue', alpha=0.4, linewidth=0.3)

        fig.suptitle('Moran`s I: '+str(round(self.mi.I,5)))
        plt.xlabel(self.fieldName)
        plt.ylabel('Spatial Lag '+self.fieldName)

        ax.set_yticks([np.mean(yl)],minor=False)
        ax.yaxis.set_major_locator(FixedLocator([round(np.mean(yl),5)]))
        ax.yaxis.grid(True)

        ax.set_xticks([np.mean(ystd), np.amax(ystd)],minor=True)
        ax.xaxis.set_major_locator(FixedLocator([round(np.mean(ystd),5)]))
        ax.xaxis.grid(True)

        """ax.spines['left'].set_position((y,np.mean(y)))"""
        """ax.spines['bottom'].set_position((yl,np.mean(yl)))"""
        """ax.set_yticklabels(['Bill', 'Jim'])"""
        """plt.grid(True)"""
        """fig.savefig('test.jpg')"""
        fig.show()

    def run(self, path):
        if self.verify():
            print "running"
            # print self.data['wtFiles'][self.data['wtFile']]
            newVars = [var.run() for var in self.newVars]
            names = [v[0] for v in newVars]
            vars = [v[1] for v in newVars]
            db = self.db()
            xid = [db.header.index(i) for i in vars]
            X = [db[:, i] for i in xid]
            W = self.loadWeights()
            lag = [pysal.lag_spatial(W, y) for y in X]
            lag = zip(*lag)  # transpose
            lag = map(list, lag)
            new_header = db.header + names

            if path.endswith('.dbf'):
                new_spec = db.field_spec + [('N', 20, 10) for n in names]
                data = db.read()
                db.close()
                newdb = pysal.open(path, 'w')
                newdb.header = new_header
                newdb.field_spec = new_spec
                for i, row in enumerate(data):
                    newdb.write(row + lag[i])
                newdb.close()

            elif path.endswith('.csv'):
                data = db.read()
                db.close()
                newdb = pysal.open(path, 'wb')
                writer = csv.writer(newdb)
                writer.writerow(new_header)
                for i, row in enumerate(data):
                    writer.writerow(row + lag[i])
                newdb.close()

    def __init__(self, y, w, permutations=0,
                 significance_level=0.05):
        y = y.transpose()
        pml = pysal.Moran_Local

        #################################################################
        # have to optimize conditional spatial permutations over a
        # time series - this is a place holder for the foreclosure paper
        ml = [pml(yi, w, permutations=permutations) for yi in y]
        #################################################################

        q = np.array([mli.q for mli in ml]).transpose()
        classes = np.arange(1, 5)  # no guarantee all 4 quadrants are visited
        Markov.__init__(self, q, classes)
        self.q = q
        self.w = w
        n, k = q.shape
        k -= 1
        self.significance_level = significance_level
        move_types = np.zeros((n, k), int)
        sm = np.zeros((n, k), int)
        self.significance_level = significance_level
        if permutations > 0:
            p = np.array([mli.p_z_sim for mli in ml]).transpose()
            self.p_values = p
            pb = p <= significance_level
        else:
            pb = np.zeros_like(y.T)
        for t in range(k):
            origin = q[:, t]
            dest = q[:, t + 1]
            p_origin = pb[:, t]
            p_dest = pb[:, t]
            for r in range(n):
                move_types[r, t] = TT[origin[r], dest[r]]
                key = (origin[r], dest[r], p_origin[r], p_dest[r])
                sm[r, t] = MOVE_TYPES[key]
        if permutations > 0:
            self.significant_moves = sm
        self.move_types = move_types

        # null of own and lag moves being independent

        ybar = y.mean(axis=0)
        r = y / ybar
        ylag = np.array([pysal.lag_spatial(w, yt) for yt in y])
        rlag = ylag / ybar
        rc = r < 1.
        rlagc = rlag < 1.
        markov_y = pysal.Markov(rc)
        markov_ylag = pysal.Markov(rlagc)
        A = np.matrix([[1, 0, 0, 0],
                       [0, 0, 1, 0],
                       [0, 0, 0, 1],
                       [0, 1, 0, 0]])

        kp = A * np.kron(markov_y.p, markov_ylag.p) * A.T
        trans = self.transitions.sum(axis=1)
        t1 = np.diag(trans) * kp
        t2 = self.transitions
        t1 = t1.getA()
        self.chi_2 = pysal.spatial_dynamics.markov.chi2(t1, t2)
        self.expected_t = t1
        self.permutations = permutations