def __init__(self, y, x, yend, q,
                 w=None,
                 robust=None, gwk=None, sig2n_k=False,
                 spat_diag=False,
                 vm=False, name_y=None, name_x=None,
                 name_yend=None, name_q=None,
                 name_w=None, name_gwk=None, name_ds=None):

        n = USER.check_arrays(y, x, yend, q)
        USER.check_y(y, n)
        USER.check_weights(w, y)
        USER.check_robust(robust, gwk)
        USER.check_spat_diag(spat_diag, w)
        x_constant = USER.check_constant(x)
        BaseTSLS.__init__(self, y=y, x=x_constant, yend=yend, q=q,
                          robust=robust, gwk=gwk, sig2n_k=sig2n_k)
        self.title = "TWO STAGE LEAST SQUARES"
        self.name_ds = USER.set_name_ds(name_ds)
        self.name_y = USER.set_name_y(name_y)
        self.name_x = USER.set_name_x(name_x, x)
        self.name_yend = USER.set_name_yend(name_yend, yend)
        self.name_z = self.name_x + self.name_yend
        self.name_q = USER.set_name_q(name_q, q)
        self.name_h = USER.set_name_h(self.name_x, self.name_q)
        self.robust = USER.set_robust(robust)
        self.name_w = USER.set_name_w(name_w, w)
        self.name_gwk = USER.set_name_w(name_gwk, gwk)
        SUMMARY.TSLS(reg=self, vm=vm, w=w, spat_diag=spat_diag)

def _work(y,x,regi_ids,r,yend,q,w_r,w_lags,lag_q,robust,sig2n_k,name_ds,name_y,name_x,name_yend,name_q,name_w,name_regimes):
    y_r = y[regi_ids[r]]
    x_r = x[regi_ids[r]]
    if yend != None:
        yend_r = yend[regi_ids[r]]
    else:
        yend_r = yend
    if q != None:
        q_r = q[regi_ids[r]]
    else:
        q_r = q
    yend_r, q_r = set_endog_sparse(y_r, x_r, w_r, yend_r, q_r, w_lags, lag_q)
    x_constant = USER.check_constant(x_r)
    if robust == 'hac':
        robust = None
    model = BaseTSLS(y_r, x_constant, yend_r, q_r, robust=robust, sig2n_k=sig2n_k)
    model.title = "SPATIAL TWO STAGE LEAST SQUARES ESTIMATION - REGIME %s" %r
    model.robust = USER.set_robust(robust)
    model.name_ds = name_ds
    model.name_y = '%s_%s'%(str(r), name_y)
    model.name_x = ['%s_%s'%(str(r), i) for i in name_x]
    model.name_yend = ['%s_%s'%(str(r), i) for i in name_yend]
    model.name_z = model.name_x + model.name_yend
    model.name_q = ['%s_%s'%(str(r), i) for i in name_q]
    model.name_h = model.name_x + model.name_q
    model.name_w = name_w
    model.name_regimes = name_regimes
    return model

def _work(y, x, w, regi_ids, r, yend, q, robust, sig2n_k, name_ds, name_y, name_x, name_yend, name_q, name_w, name_regimes):
    y_r = y[regi_ids[r]]
    x_r = x[regi_ids[r]]
    yend_r = yend[regi_ids[r]]
    q_r = q[regi_ids[r]]
    x_constant = USER.check_constant(x_r)
    if robust == 'hac' or robust == 'ogmm':
        robust2 = None
    else:
        robust2 = robust
    model = BaseTSLS(y_r, x_constant, yend_r, q_r,
                     robust=robust2, sig2n_k=sig2n_k)
    model.title = "TWO STAGE LEAST SQUARES ESTIMATION - REGIME %s" % r
    if robust == 'ogmm':
        _optimal_weight(model, sig2n_k, warn=False)
    model.robust = USER.set_robust(robust)
    model.name_ds = name_ds
    model.name_y = '%s_%s' % (str(r), name_y)
    model.name_x = ['%s_%s' % (str(r), i) for i in name_x]
    model.name_yend = ['%s_%s' % (str(r), i) for i in name_yend]
    model.name_z = model.name_x + model.name_yend
    model.name_q = ['%s_%s' % (str(r), i) for i in name_q]
    model.name_h = model.name_x + model.name_q
    model.name_w = name_w
    model.name_regimes = name_regimes
    if w:
        w_r, warn = REGI.w_regime(w, regi_ids[r], r, transform=True)
        set_warn(model, warn)
        model.w = w_r
    return model

    def __init__(self, y, x, yend=None, q=None,\
                 w=None, w_lags=1, lag_q=True,\
                 robust=None, gwk=None, sig2n_k=False,\
                 spat_diag=False,\
                 vm=False, name_y=None, name_x=None,\
                 name_yend=None, name_q=None,\
                 name_w=None, name_gwk=None, name_ds=None):

        n = USER.check_arrays(x, yend, q)
        USER.check_y(y, n)
        USER.check_weights(w, y, w_required=True)
        USER.check_robust(robust, gwk)
        yend2, q2 = set_endog(y, x, w, yend, q, w_lags, lag_q)
        x_constant = USER.check_constant(x)
        BaseGM_Lag.__init__(self, y=y, x=x_constant, w=w.sparse, yend=yend2, q=q2,\
                            w_lags=w_lags, robust=robust, gwk=gwk,\
                            lag_q=lag_q, sig2n_k=sig2n_k)
        self.predy_e, self.e_pred, warn = sp_att(w,self.y,self.predy,\
                      yend2[:,-1].reshape(self.n,1),self.betas[-1])
        set_warn(self,warn)
        self.title = "SPATIAL TWO STAGE LEAST SQUARES"        
        self.name_ds = USER.set_name_ds(name_ds)
        self.name_y = USER.set_name_y(name_y)
        self.name_x = USER.set_name_x(name_x, x)
        self.name_yend = USER.set_name_yend(name_yend, yend)
        self.name_yend.append(USER.set_name_yend_sp(self.name_y))
        self.name_z = self.name_x + self.name_yend
        self.name_q = USER.set_name_q(name_q, q)
        self.name_q.extend(USER.set_name_q_sp(self.name_x, w_lags, self.name_q, lag_q))
        self.name_h = USER.set_name_h(self.name_x, self.name_q)
        self.robust = USER.set_robust(robust)
        self.name_w = USER.set_name_w(name_w, w)
        self.name_gwk = USER.set_name_w(name_gwk, gwk)
        SUMMARY.GM_Lag(reg=self, w=w, vm=vm, spat_diag=spat_diag)

 def __init__(self, y, x, w, method='full', epsilon=0.0000001,
              spat_diag=False, vm=False, name_y=None, name_x=None,
              name_w=None, name_ds=None):
     n = USER.check_arrays(y, x)
     USER.check_y(y, n)
     USER.check_weights(w, y, w_required=True)
     x_constant = USER.check_constant(x)
     method = method.upper()
     if method in ['FULL', 'ORD']:
         BaseML_Lag.__init__(self, y=y, x=x_constant,
                             w=w, method=method, epsilon=epsilon)
         # increase by 1 to have correct aic and sc, include rho in count
         self.k += 1
         self.title = "MAXIMUM LIKELIHOOD SPATIAL LAG" + \
             " (METHOD = " + method + ")"
         self.name_ds = USER.set_name_ds(name_ds)
         self.name_y = USER.set_name_y(name_y)
         self.name_x = USER.set_name_x(name_x, x)
         name_ylag = USER.set_name_yend_sp(self.name_y)
         self.name_x.append(name_ylag)  # rho changed to last position
         self.name_w = USER.set_name_w(name_w, w)
         self.aic = DIAG.akaike(reg=self)
         self.schwarz = DIAG.schwarz(reg=self)
         SUMMARY.ML_Lag(reg=self, w=w, vm=vm, spat_diag=spat_diag)
     else:
         raise Exception, "{0} is an unsupported method".format(method)

    def __init__(self, y, x, yend=None, q=None,\
                 w=None, w_lags=1, lag_q=True,\
                 vm=False, name_y=None, name_x=None,\
                 name_yend=None, name_q=None,\
                 name_w=None, name_ds=None):

        n = USER.check_arrays(y, x, yend, q)
        USER.check_y(y, n)
        USER.check_weights(w, y, w_required=True)
        yend2, q2 = set_endog(y, x, w, yend, q, w_lags, lag_q)
        x_constant = USER.check_constant(x)
        BaseGM_Combo.__init__(self, y=y, x=x_constant, w=w.sparse, yend=yend2, q=q2,\
                                w_lags=w_lags, lag_q=lag_q)
        self.rho = self.betas[-2]
        self.predy_e, self.e_pred, warn = sp_att(w,self.y,\
                   self.predy,yend2[:,-1].reshape(self.n,1),self.rho)
        set_warn(self, warn)
        self.title = "SPATIALLY WEIGHTED TWO STAGE LEAST SQUARES"        
        self.name_ds = USER.set_name_ds(name_ds)
        self.name_y = USER.set_name_y(name_y)
        self.name_x = USER.set_name_x(name_x, x)
        self.name_yend = USER.set_name_yend(name_yend, yend)
        self.name_yend.append(USER.set_name_yend_sp(self.name_y))
        self.name_z = self.name_x + self.name_yend
        self.name_z.append('lambda')
        self.name_q = USER.set_name_q(name_q, q)
        self.name_q.extend(USER.set_name_q_sp(self.name_x, w_lags, self.name_q, lag_q))
        self.name_h = USER.set_name_h(self.name_x, self.name_q)
        self.name_w = USER.set_name_w(name_w, w)
        SUMMARY.GM_Combo(reg=self, w=w, vm=vm)   

def _run_stp1(y,x,eq_ids,r):
    y_r = y[eq_ids[r]]
    x_r = x[eq_ids[r]]     
    x_constant = USER.check_constant(x_r)
    model = BaseOLS(y_r, x_constant)
    #model.logll = diagnostics.log_likelihood(model) 
    return model            

    def __init__(self, y, x, regimes, w=None, constant_regi='many',
                 cols2regi='all', method='full', epsilon=0.0000001,
                 regime_err_sep=False, cores=None, spat_diag=False,
                 vm=False, name_y=None, name_x=None,
                 name_w=None, name_ds=None, name_regimes=None):

        n = USER.check_arrays(y, x)
        USER.check_y(y, n)
        USER.check_weights(w, y, w_required=True)
        self.constant_regi = constant_regi
        self.cols2regi = cols2regi
        self.regime_err_sep = regime_err_sep
        self.name_ds = USER.set_name_ds(name_ds)
        self.name_y = USER.set_name_y(name_y)
        self.name_w = USER.set_name_w(name_w, w)
        self.name_regimes = USER.set_name_ds(name_regimes)
        self.n = n
        self.y = y

        x_constant = USER.check_constant(x)
        name_x = USER.set_name_x(name_x, x)
        self.name_x_r = name_x

        cols2regi = REGI.check_cols2regi(constant_regi, cols2regi, x)
        self.regimes_set = REGI._get_regimes_set(regimes)
        self.regimes = regimes
        USER.check_regimes(self.regimes_set, self.n, x.shape[1])
        self.regime_err_sep = regime_err_sep

        if regime_err_sep == True:
            if set(cols2regi) == set([True]):
                self._error_regimes_multi(y, x, regimes, w, cores,
                                          method, epsilon, cols2regi, vm, name_x, spat_diag)
            else:
                raise Exception, "All coefficients must vary accross regimes if regime_err_sep = True."
        else:
            regimes_att = {}
            regimes_att['x'] = x_constant
            regimes_att['regimes'] = regimes
            regimes_att['cols2regi'] = cols2regi
            x, name_x = REGI.Regimes_Frame.__init__(self, x_constant,
                                                    regimes, constant_regi=None, cols2regi=cols2regi,
                                                    names=name_x)

            BaseML_Error.__init__(
                self, y=y, x=x, w=w, method=method, epsilon=epsilon, regimes_att=regimes_att)

            self.title = "MAXIMUM LIKELIHOOD SPATIAL ERROR - REGIMES" + \
                " (METHOD = " + method + ")"
            self.name_x = USER.set_name_x(name_x, x, constant=True)
            self.name_x.append('lambda')
            self.kf += 1  # Adding a fixed k to account for lambda.
            self.chow = REGI.Chow(self)
            self.aic = DIAG.akaike(reg=self)
            self.schwarz = DIAG.schwarz(reg=self)
            self._cache = {}
            SUMMARY.ML_Error(reg=self, w=w, vm=vm,
                             spat_diag=spat_diag, regimes=True)

 def _get_spat_diag_props(self,y, x, w, yend, q, w_lags, lag_q):
     self._cache = {}
     yend, q = set_endog(y, x, w, yend, q, w_lags, lag_q)
     x = USER.check_constant(x)
     x = REGI.regimeX_setup(x, self.regimes, [True] * x.shape[1], self.regimes_set)
     self.z = sphstack(x,REGI.regimeX_setup(yend, self.regimes, [True] * (yend.shape[1]-1)+[False], self.regimes_set))
     self.h = sphstack(x,REGI.regimeX_setup(q, self.regimes, [True] * q.shape[1], self.regimes_set))
     hthi = np.linalg.inv(spdot(self.h.T,self.h))
     zth = spdot(self.z.T,self.h)     
     self.varb = np.linalg.inv(spdot(spdot(zth,hthi),zth.T))

 def _get_spat_diag_props(self, results, regi_ids, x, yend, q):
     self._cache = {}
     x = USER.check_constant(x)
     x = REGI.regimeX_setup(
         x, self.regimes, [True] * x.shape[1], self.regimes_set)
     self.z = sphstack(x, REGI.regimeX_setup(
         yend, self.regimes, [True] * yend.shape[1], self.regimes_set))
     self.h = sphstack(
         x, REGI.regimeX_setup(q, self.regimes, [True] * q.shape[1], self.regimes_set))
     hthi = np.linalg.inv(spdot(self.h.T, self.h))
     zth = spdot(self.z.T, self.h)
     self.varb = np.linalg.inv(spdot(spdot(zth, hthi), zth.T))

def _work(
    y,
    x,
    regi_ids,
    r,
    yend,
    q,
    w_r,
    w_lags,
    lag_q,
    robust,
    sig2n_k,
    name_ds,
    name_y,
    name_x,
    name_yend,
    name_q,
    name_w,
    name_regimes,
):
    y_r = y[regi_ids[r]]
    x_r = x[regi_ids[r]]
    if yend != None:
        yend_r = yend[regi_ids[r]]
    else:
        yend_r = yend
    if q != None:
        q_r = q[regi_ids[r]]
    else:
        q_r = q
    yend_r, q_r = set_endog_sparse(y_r, x_r, w_r, yend_r, q_r, w_lags, lag_q)
    x_constant = USER.check_constant(x_r)
    if robust == "hac" or robust == "ogmm":
        robust2 = None
    else:
        robust2 = robust
    model = BaseTSLS(y_r, x_constant, yend_r, q_r, robust=robust2, sig2n_k=sig2n_k)
    model.title = "SPATIAL TWO STAGE LEAST SQUARES ESTIMATION - REGIME %s" % r
    if robust == "ogmm":
        _optimal_weight(model, sig2n_k, warn=False)
    model.rho = model.betas[-1]
    model.robust = USER.set_robust(robust)
    model.name_ds = name_ds
    model.name_y = "%s_%s" % (str(r), name_y)
    model.name_x = ["%s_%s" % (str(r), i) for i in name_x]
    model.name_yend = ["%s_%s" % (str(r), i) for i in name_yend]
    model.name_z = model.name_x + model.name_yend
    model.name_q = ["%s_%s" % (str(r), i) for i in name_q]
    model.name_h = model.name_x + model.name_q
    model.name_w = name_w
    model.name_regimes = name_regimes
    return model

def _work(y, x, regi_ids, r, w_r, method, epsilon, name_ds, name_y, name_x, name_w, name_regimes):
    y_r = y[regi_ids[r]]
    x_r = x[regi_ids[r]]
    x_constant = USER.check_constant(x_r)
    model = BaseML_Lag(y_r, x_constant, w_r, method=method, epsilon=epsilon)
    model.title = "MAXIMUM LIKELIHOOD SPATIAL LAG - REGIME " + str(r) + " (METHOD = " + method + ")"
    model.name_ds = name_ds
    model.name_y = "%s_%s" % (str(r), name_y)
    model.name_x = ["%s_%s" % (str(r), i) for i in name_x]
    model.name_w = name_w
    model.name_regimes = name_regimes
    model.aic = DIAG.akaike(reg=model)
    model.schwarz = DIAG.schwarz(reg=model)
    return model

def _work(y,x,regi_ids,r,robust,sig2n_k,name_ds,name_y,name_x,name_w,name_regimes):
    y_r = y[regi_ids[r]]
    x_r = x[regi_ids[r]]
    x_constant = USER.check_constant(x_r)
    if robust == 'hac':
        robust = None
    model = BaseOLS(y_r, x_constant, robust=robust, sig2n_k=sig2n_k)
    model.title = "ORDINARY LEAST SQUARES ESTIMATION - REGIME %s" %r
    model.robust = USER.set_robust(robust)
    model.name_ds = name_ds
    model.name_y = '%s_%s'%(str(r), name_y)
    model.name_x = ['%s_%s'%(str(r), i) for i in name_x]
    model.name_w = name_w
    model.name_regimes = name_regimes
    return model

    def __init__(self, y, x, w=None, optim='newton',scalem='phimean',maxiter=100,\
                 vm=False, name_y=None, name_x=None, name_w=None, name_ds=None, \
                 spat_diag=False):

        n = USER.check_arrays(y, x)
        USER.check_y(y, n)
        if w:
            USER.check_weights(w, y)
            spat_diag = True
        x_constant = USER.check_constant(x)
        BaseProbit.__init__(self,y=y,x=x_constant,w=w,optim=optim,scalem=scalem,maxiter=maxiter) 
        self.title = "CLASSIC PROBIT ESTIMATOR"        
        self.name_ds = USER.set_name_ds(name_ds)    
        self.name_y = USER.set_name_y(name_y)
        self.name_x = USER.set_name_x(name_x, x)
        self.name_w = USER.set_name_w(name_w, w)
        SUMMARY.Probit(reg=self, w=w, vm=vm, spat_diag=spat_diag)    

def _work_error(y,x,regi_ids,r,w,method,epsilon,name_ds,name_y,name_x,name_w,name_regimes):
    w_r,warn = REGI.w_regime(w, regi_ids[r], r, transform=True)
    y_r = y[regi_ids[r]]
    x_r = x[regi_ids[r]]
    x_constant = USER.check_constant(x_r)
    model = BaseML_Error(y=y_r,x=x_constant,w=w_r,method=method,epsilon=epsilon)
    set_warn(model, warn)
    model.w = w_r
    model.title = "MAXIMUM LIKELIHOOD SPATIAL ERROR - REGIME "+str(r)+" (METHOD = "+method+")"
    model.name_ds = name_ds
    model.name_y = '%s_%s'%(str(r), name_y)
    model.name_x = ['%s_%s'%(str(r), i) for i in name_x]
    model.name_w = name_w
    model.name_regimes = name_regimes
    model.aic = DIAG.akaike(reg=model)
    model.schwarz = DIAG.schwarz(reg=model)
    return model

def _run_stp1(y,x,yend,q,eq_ids,r,sig2n_k,w,w_lags,lag_q):
    y_r = y[eq_ids[r]]
    x_r = x[eq_ids[r]]
    yend_r = None
    if yend != None:
        yend_r = yend[eq_ids[r]]
        q_r = q[eq_ids[r]]
    if w_lags:
        yend_r, q_r = set_endog_sparse(y_r, x_r, w, yend_r, q_r, w_lags, lag_q)        
    x_constant = USER.check_constant(x_r)
    if yend != None:
        model = BaseTSLS(y_r, x_constant, yend_r, q_r, sig2n_k=sig2n_k) 
        model.zhat = np.dot(model.h,np.dot(model.hthi,model.htz))
    else:
        model = BaseOLS(y_r, x_constant, sig2n_k=sig2n_k) 
        model.zhat = model.x
        model.z = model.x
    return model

 def __init__(self, y, x, w, method='full', epsilon=0.0000001,
              spat_diag=False, vm=False, name_y=None, name_x=None,
              name_w=None, name_ds=None):
     n = USER.check_arrays(y, x)
     USER.check_y(y, n)
     USER.check_weights(w, y, w_required=True)
     x_constant = USER.check_constant(x)
     method = method.upper()
     BaseML_Error.__init__(self, y=y, x=x_constant,
                           w=w, method=method, epsilon=epsilon)
     self.title = "MAXIMUM LIKELIHOOD SPATIAL ERROR" + \
         " (METHOD = " + method + ")"
     self.name_ds = USER.set_name_ds(name_ds)
     self.name_y = USER.set_name_y(name_y)
     self.name_x = USER.set_name_x(name_x, x)
     self.name_x.append('lambda')
     self.name_w = USER.set_name_w(name_w, w)
     self.aic = DIAG.akaike(reg=self)
     self.schwarz = DIAG.schwarz(reg=self)
     SUMMARY.ML_Error(reg=self, w=w, vm=vm, spat_diag=spat_diag)

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 __init__(self, y, x,\
                 w=None,\
                 robust=None, gwk=None, sig2n_k=True,\
                 nonspat_diag=True, spat_diag=False, moran=False,\
                 vm=False, name_y=None, name_x=None,\
                 name_w=None, name_gwk=None, name_ds=None):

        n = USER.check_arrays(y, x)
        USER.check_y(y, n)
        USER.check_weights(w, y)
        USER.check_robust(robust, gwk)
        USER.check_spat_diag(spat_diag, w)
        x_constant = USER.check_constant(x)
        BaseOLS.__init__(self, y=y, x=x_constant, robust=robust,\
                     gwk=gwk, sig2n_k=sig2n_k) 
        self.title = "ORDINARY LEAST SQUARES"
        self.name_ds = USER.set_name_ds(name_ds)
        self.name_y = USER.set_name_y(name_y)
        self.name_x = USER.set_name_x(name_x, x)
        self.robust = USER.set_robust(robust)
        self.name_w = USER.set_name_w(name_w, w)
        self.name_gwk = USER.set_name_w(name_gwk, gwk)
        SUMMARY.OLS(reg=self, vm=vm, w=w, nonspat_diag=nonspat_diag,\
                    spat_diag=spat_diag, moran=moran)