%matplotlib inline
import pandas as pd
import numpy as np
import pysal as ps
import matplotlib.pyplot as plt

db = pd.read_csv('../workshop_data/Houston_pop00.csv')
db.info()

db.head()

db.tail()

db.describe().T

downtown = db[db['dcbd'] < 15]
downtown.info()

l = [1, 2, 3, 4]
a = np.array(l)
a

print type(a[0])
l += 'a'
a = np.array(l)
print type(a[0])

l = [[1, 2], [3, 4], [5, 6]]
a = np.array(l)
print 'Array a has a dimension of: ', a.shape
print a

a = np.random.random((3, 2))
b = np.random.random((2, 3))
a, b

# Sum (note the transpose for dimensionality alignment)
a + b.T

# Difference (note the transpose for dimensionality alignment)
a - b.T

# Matrix product
c = np.dot(a, b)
c

# Matrix division by a scalar
c = a / 2.
c

db['pop_dens'] = db['POP00'] / db['ALAND']

db['pop_dens'].plot(kind='kde')

db['dcbd'].hist(bins=50, color='k', alpha=0.5, grid=False)

db.pop_dens = db['POP00'] / db['ALAND']
db['pop_dens'].describe()

dbf = ps.open('../workshop_data/houston_tract_pop_emp_wgs84.dbf')
dbf.header

df = pd.DataFrame({'emp_dens': dbf.by_col('emp_dens'), \
                   'pop_dens': dbf.by_col('pop_dens'), \
                   'dcbd': dbf.by_col('dcbd'), \
                   'downtown': dbf.by_col('downtown')})
df.info()

w = ps.queen_from_shapefile('../../sdar_mini_repo/data/houston_tract_pop_emp_wgs84.shp')

w

w.n

w[0]

w.transform = 'R'

w[0]

f = ps.open('houston_tract_queen.gal', 'w')
f.write(w)
f.close()

from IPython.display import HTML
HTML('<iframe src=http://pysal.readthedocs.org/en/v1.7/library/weights/Wsets.html#pysal.weights.Wsets width=100% height=350></iframe>')

# Downtown/suburbs weights
dt_sb = ps.weights.block_weights(dbf.by_col('downtown'))
dt_sb = ps.weights.regime_weights(dbf.by_col('downtown'))

# Queen example
queen = ps.queen_from_shapefile('../../sdar_mini_repo/data/houston_tract_pop_emp_wgs84.shp')

w = ps.weights.Wsets.w_intersection(queen, dt_sb)

w.transform = 'R'

wy = ps.lag_spatial(w, df['emp_dens'])

from pysal.contrib.viz import mapping as maps

shp_link = '../../sdar_mini_repo/data/houston_tract_pop_emp_wgs84.shp'
maps.plot_choropleth(shp_link, df['downtown'], 'unique_values', \
                     figsize=(12, 12))

shp_link = '../../sdar_mini_repo/data/houston_tract_pop_emp_wgs84.shp'
maps.plot_choropleth(shp_link, df['emp_dens'], 'quantiles', figsize=(12, 12))

mi = ps.Moran(df['emp_dens'], w)

mi.I

lmi = ps.Moran_Local(df['emp_dens'].values, w)

lmi.p_sim

ols_base = ps.spreg.OLS(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'dcbd']].values, \
                    name_x = ['pop_dens', 'dcbd'], name_y='emp_dens')

print ols_base.summary

ols_base.betas

ols_base.std_err

ols_white = ps.spreg.OLS(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'dcbd']].values, \
                    robust='white', \
                    name_x = ['pop_dens', 'dcbd'], name_y='emp_dens')

ols_white.std_err

ols_sp_diag = ps.spreg.OLS(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'dcbd']].values, w, \
                    name_x = ['pop_dens', 'dcbd'], name_y='emp_dens', \
                    spat_diag=True)

ols_sp_diag.lm_error

ols_sp_diag.lm_lag

print ols_sp_diag.summary

fes = pd.get_dummies(df['downtown'], prefix='downtown')
fes.head()

x = df[['pop_dens', 'dcbd']].join(fes.drop('downtown_0', axis=1))
ols_fe = ps.spreg.OLS(df['emp_dens'].values[:, None], x.values, w, \
                    name_x = list(x.columns), name_y='emp_dens', \
                    spat_diag=True)

print ols_fe.summary

downtown = df['downtown'].map({1: 'downtown', 0: 'suburbs'})

ols_regimes = ps.spreg.OLS_Regimes(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'dcbd']].values, downtown, w=w, \
                    name_x = ['pop_dens', 'dcbd'], name_y='emp_dens', \
                    spat_diag=True)

print ols_regimes.summary

df['w_pop_dens'] = ps.lag_spatial(w, df['pop_dens'])

ols_wx = ps.spreg.OLS(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'w_pop_dens', 'dcbd']].values, \
                    name_x = ['pop_dens', 'w_pop_dens', 'dcbd'], name_y='emp_dens')

print ols_wx.summary

lag_iv = ps.spreg.GM_Lag(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'dcbd']].values, w=w, \
                    name_x = ['pop_dens', 'dcbd'], name_y='emp_dens')

print lag_iv.summary

lag_ml = ps.spreg.ML_Lag(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'dcbd']].values, w=w, \
                    name_x = ['pop_dens', 'dcbd'], name_y='emp_dens', \
                    method='ord')

print lag_ml.summary

ps.spreg.GM_Endog_Error_Het?

error_gmm_arraiz = ps.spreg.GM_Error_Het(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'dcbd']].values, w=w, \
                    name_x = ['pop_dens', 'dcbd'], name_y='emp_dens')

print error_gmm_arraiz.summary

error_gmm_drucker = ps.spreg.GM_Error_Hom(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'dcbd']].values, w=w, \
                    name_x = ['pop_dens', 'dcbd'], name_y='emp_dens')

print error_gmm_drucker.summary

error_ml = ps.spreg.ML_Error(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'dcbd']].values, w=w, \
                    name_x = ['pop_dens', 'dcbd'], name_y='emp_dens')

print error_ml.summary

wk = ps.kernelW_from_shapefile('../../sdar_mini_repo/data/houston_tract_pop_emp_wgs84.shp', \
                                  k=15,function='triangular', fixed=False)

ols_hac = ps.spreg.OLS(df['emp_dens'].values[:, None], \
                    df[['pop_dens', 'dcbd']].values, \
                    name_x = ['pop_dens', 'dcbd'], name_y='emp_dens', \
                    robust='hac', gwk=wk)

print ols_hac.summary

#Import the LM tests method
from pysal.spreg import LMtests
#Specify the files for the weights we want to try as a list
ws = [w, queen, dt_sb]
#Run the OLS
model = ps.spreg.OLS(df['emp_dens'].values[:, None], \
                     df[['pop_dens', 'dcbd']].values, \
                     spat_diag=False, nonspat_diag=False)
#Setup the loop over the weights files
for weights in ws:
    lms = LMtests(model, w)
    print '\tLM error: %.4f\t(%.4f)'%lms.lme
    print '\tLM lag:   %.4f\t(%.4f)'%lms.lml
    print '\tSARMA:    %.4f\t(%.4f)'%lms.sarma
    print '\tRobust LM error:   %.4f\t(%.4f)'%lms.rlme
    print '\tRobust LM lag:     %.4f\t(%.4f)'%lms.rlml
    print '----------------\n'