def stan_cache(model_code, model_name=None, **kwargs):
    """Use just as you would `stan`"""
    code_hash = md5(model_code.encode('ascii')).hexdigest()
    if model_name is None:
        cache_fn = 'cached-model-{}.pkl'.format(code_hash)
    else:
        cache_fn = 'cached-{}-{}.pkl'.format(model_name, code_hash)
    try:
        sm = pickle.load(open(cache_fn, 'rb'))
    except:
        sm = pystan.StanModel(model_code=model_code)
        with open(cache_fn, 'wb') as f:
            pickle.dump(sm, f)
    else:
        print("Using cached StanModel")
    return sm.sampling(**kwargs)

def toenglish(s):
    spanish = ['ene', 'abr', 'ago', 'dic', 'de mayo de']
    english = ['jan', 'apr', 'aug', 'dec', 'may']
    for (j, month) in enumerate(spanish):
        s = s.replace(month, english[j])
    return s

def getPercentage(s):
    if (s[0] not in ['0','1','2','3','4','5','6','7','8','9']):
        return 0
    else:
        if (s.find('%') != -1):
            return float(s.split('%')[0].replace(',','.')) / 100.0
        else:
            return float(s.split('\n')[0].replace(',','.')) / 100.0

def getSigma(s):
    left = s.find('(')
    right = s.find(')')    
    if (s[left+1:right] in ['?', '-']):
        return 0.03
    else:
        return 1.0 / np.sqrt(float(s[left+1:right]))

def weeksDifference(d1, d2):
    monday1 = (d1 - datetime.timedelta(days=d1.weekday()))
    monday2 = (d2 - datetime.timedelta(days=d2.weekday()))

    return int((monday2 - monday1).days / 7)

"""
Explanation: Análisis de los sondeos electorales 2016
Vamos a plantearnos como ejercicio analizar los sondeos electorales. Para ello cogemos como datos los sondeos recopilados en la página web de la Wikipedia para las elecciones de 2015 y para las de 2016. Consideramos sólo los partidos PP, PSOE, IU, Podemos y Ciudadanos. El resto de partidos se consideran metidos dentro de un mismo saco, que vienen siendo alrededor del 15% en cada sondeo.
Para modelar la evolución temporal de los sondeos para cada partido, usamos un modelado con variables latentes. La idea es que los resultados de los sondeos son sólo una estimación ruidosa y sesgada de la intención de voto real de los votantes. Por tanto, considamos que $\theta_i(t_j)$ es la subyacente intención de voto para el partido $i$ en el tiempo $t_j$. Esta variable se observa de forma ruidosa para cada empresa $k$ que realiza el sondeo, y además incluye un cierto sesgo. Por tanto, nuestro modelo de observación es:
$$
y_{ik}(t_j) \sim N(\theta_i(t_j) + b_{ik}, \sigma_{k}(t_j)^2)
$$
donde $b_{ik}$ es el bias que introduce la empresa de sondeos $k$ sobre el partido $i$ y $\sigma_{ik}^2$ es la varianza del sondeo realizado en el tiempo $t_j$ por la empresa $k$. Para los resultados de las elecciones, elegimos una varianza muy pequeña y no incluimos los sesgos, porque es el único momento en el que tenemos acceso directo a las intenciones subyacentes:
$$
y_{ik}(\mathrm{elecciones}) \sim N(\theta_i(t_j), \sigma_{k}(\mathrm{elecciones})^2).
$$
Es obvio que la suma de porcentajes de cada partido debe sumar 1 en cada sondeo, por lo cual debemos asegurarnos de que
$$
\sum_i \theta_i(t_j) = 1.
$$
Además, es lógico asumir que, dadas suficientes casas de sondeo, los biases sumen 0 para cada partido:
$$
\sum_k b_{ik} = 0.
$$
Por último, asumimos un proceso autoregresivo Markoviano para modelar la variación temporal de la intención de voto real para cada partido. Es decir, en cada sondeo, la intención de voto estará dada por una variable Gaussiana con media igual a la intención de voto anterior, y una varianza que debemos estimar para cada partido:
$$
\theta_i(t_j) \sim N(\theta_i(t_{j-1}), \sigma_i^2)
$$
Como es un parámetro de escala, definimos un prior poco informativo sobre $\sigma_i$:
$$
\sigma_i \sim IG(0.01,0.01).
$$
Por último, debemos definir un prior para $\theta_i(t_0)$, para el cual escogemos una distribución de Dirichlet con $\alpha=1$:
$$
\theta_i(t_0) \sim \mathrm{Dir}(\alpha).
$$
Modelo
Definimos algunas funciones útiles para el parsing de la base de datos.
End of explanation
"""


wb = openpyxl.load_workbook("data/sondeos.xlsx")
ws = wb.active

empresas = ['GAD3', 'Encuestamos', 'GESOP', 'Metroscopia', 'Celeste-Tel','Demoscopia Servicios', 'Simple Lógica', 'CIS', 
            'TNS Demoscopia', 'Invymark', 'NC Report', 'El Espanol', 'DYM', 'Sondaxe', 'Sigma Dos', 'Resultados de las elecciones']

empresaSondeoAll = []
sondeosAll = []
dateAll = []
sigmaAll = []

wbnew = openpyxl.Workbook()

# grab the active worksheet
wsnew = wbnew.active
wsnew.append(['Encuestador','Fecha','Muestra','PP','PSOE','PODEMOS','Cs','IU'])

otrosPartidos = ['F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'R']

nSondeos = 302
loopSondeos = 0
for i in range(130):
    empresa = ws['A{0}'.format(i+2)].value
    for (loop, emp) in enumerate(empresas):
        if (empresa.find(emp) != -1):
            empresaSondeo = loop

    if (empresaSondeo == len(empresas)-1):
        sigma = 0.0001
    else:
        sigma = getSigma(empresa)
        
    PP = getPercentage(ws['C{0}'.format(i+2)].value)
    PSOE = getPercentage(ws['D{0}'.format(i+2)].value)
    IU = getPercentage(ws['E{0}'.format(i+2)].value)
    PODEMOS = getPercentage(ws['P{0}'.format(i+2)].value)
    CS = getPercentage(ws['Q{0}'.format(i+2)].value)

    total = PP + PSOE + IU + PODEMOS + CS
    otros = 1.0 - total    

    tmp = ws['B{0}'.format(i+2)].value
    if (isinstance(tmp, datetime.date)):
        date = tmp
    else:
        date = dateutil.parser.parse(toenglish(ws['B{0}'.format(i+2)].value.split('-')[-1].lower()))    

    tmp = date.year + (date.month-1.0) / 12.0

    sondeo = [PP, PSOE, IU+PODEMOS, CS, otros]
    
    if (0 not in sondeo):
        
        loopSondeos += 1

        sondeosAll.append(sondeo)
        sigmaAll.append(sigma)
    
        dateAll.append(date)
        empresaSondeoAll.append(empresaSondeo+1)

        print ("{0} - {1} {7} - s={8:4.2f} : PP={2:4.2f} - PSOE={3:4.2f} - UP={4:4.2f} - CS={5:4.2f} - Resto={6:4.2f}".format(i, 
            empresas[empresaSondeo], PP*100, PSOE*100, (IU+PODEMOS)*100, CS*100, otros*100, date, sigma*100))
        
        wsnew.append([empresas[empresaSondeo],'{0}/{1}/{2}'.format(date.day,date.month,date.year-2000),1.0/sigma**2,PP*100,PSOE*100,(IU+PODEMOS)*100,CS*100,0.0])       
        
wbnew.save('new.xlsx')
sondeosAll = np.array(sondeosAll)

nSondeos, nPartidos = sondeosAll.shape
nEmpresas = len(empresas)

print ("Número de sondeos incluidos : {0} de {1} posibles".format(loopSondeos,120))

# Compute week of every poll
weekAll = []
for i in range(nSondeos):
    weekAll.append(weeksDifference(dateAll[nSondeos-1], dateAll[i]) + 1)

nDates = max(weekAll)

# Reverse all lists
sondeosAll = sondeosAll[::-1]
empresaSondeoAll.reverse()
weekAll.reverse()
sigmaAll.reverse()

dictionary = {'NPartidos': nPartidos, 'NSondeos': nSondeos, 'NEmpresas': nEmpresas-1, 'NDates': nDates, 'empresa': empresaSondeoAll, 'sondeos': sondeosAll, 
'date': weekAll, 'sigmaSondeo': sigmaAll, 'alpha': np.ones(nPartidos)*1.0}

"""
Explanation: Leemos el fichero de datos obtenido de la Wikipedia y hacemos el parsing. Escogemos sólo aquellos sondeos en los que hay información para todos los partidos. Podríamos también escoger aquellos parciales porque el modelado Bayesiano los podría aceptar.
End of explanation
"""


date.year-2000

f = open('model.stan', 'r')
model = f.read()
f.close()

out = stan_cache(model, model_name='elecciones', data=dictionary, chains=1)

"""
Explanation: Leemos el modelo que hemos definido y hacemos el sampling.
End of explanation
"""


thetaChain = out.extract('theta')['theta']
theta = np.percentile(thetaChain, [50.0, 50.0-68/2.0, 50.0+68/2], axis=0)
houseChain = out.extract('house')['house']
house = np.percentile(houseChain, [50.0, 50.0-68/2.0, 50.0+68/2], axis=0)
sigmaChain = out.extract('sigma')['sigma']
sigma = np.percentile(sigmaChain, [50.0, 50.0-68/2.0, 50.0+68/2], axis=0)

"""
Explanation: Resultados
Extraemos algunas de las cantidades interesantes para hacer algunas representaciones.
End of explanation
"""


colors = ["blue", "red", "violet", "orange", "yellow"]
labelsPartidos = ['PP', 'PSOE', 'PODEMOS+IU', 'Cs', 'Otros']
labelsEmpresas = ['GAD3', 'Encuest', 'GESOP', 'Metroscopia', 'Celeste',' Demosc. S.', 'S. Lógica', 'CIS', 'TNS Demos.', 
                  'Invymark', 'NC Rep.', 'Netquest']
f, ax = pl.subplots(figsize=(13,8))
for i in range(5):
    ax.plot(weekAll, sondeosAll[:,i], '.', color=sns.xkcd_rgb[colors[i]], linewidth=2)
    ax.plot(np.arange(max(weekAll))+1, theta[0,:,i], color=sns.xkcd_rgb[colors[i]], linewidth=2)
    ax.fill_between(np.arange(max(weekAll))+1, theta[1,:,i], theta[2,:,i], color=sns.xkcd_rgb[colors[i]], alpha=0.3)
ax.set_xlabel('Semana')
ax.set_ylabel('Fracción de votos')
ax.axvline(49)
pl.savefig("sondeos.png")

"""
Explanation: Sondeos y variación subyacente
En esta representación vemos la variación de las estimaciones de voto en puntos para cada una de las encuestas que hemos considerado y la estimación subyacente que hemos inferido. Como se ve, la varianza de la estimación subyacente cambia con el tiempo, siendo muy pequeña en la fecha de las elecciones, marcada con una línea vertical.
End of explanation
"""


f, ax = pl.subplots()
sns.boxplot(data=thetaChain[:,-1,:], ax=ax)
ax.set_xticks(np.arange(5))
ax.set_xticklabels(labelsPartidos, rotation=90)
boxes = ax.artists
for (i, box) in enumerate(boxes):
    box.set_facecolor(sns.xkcd_rgb[colors[i]])
pl.savefig("estimacionActual.png")

"""
Explanation: Estimación de voto actual
A continuación pongo la estimación de voto actual y sus intervalos de confianza.
End of explanation
"""


labelsPartidos = ['PP', 'PSOE', 'PODEMOS+IU', 'Cs', 'Otros']
f, ax = pl.subplots(nrows=2, ncols=3, figsize=(10,8))
ax = ax.flatten()
for i in range(5):
    ax[i].plot(weekAll, sondeosAll[:,i], '.', color=sns.xkcd_rgb[colors[i]], linewidth=2)
    ax[i].set_ylim([0,0.35])
    ax[i].errorbar(weekAll, sondeosAll[:,i], fmt='none', yerr=sigmaAll, color=sns.xkcd_rgb[colors[i]], linewidth=2, ecolor=sns.xkcd_rgb[colors[i]])
    ax[i].plot(np.arange(max(weekAll))+1, theta[0,:,i], color=sns.xkcd_rgb[colors[i]], linewidth=2)
    ax[i].fill_between(np.arange(max(weekAll))+1, theta[1,:,i], theta[2,:,i], color=sns.xkcd_rgb[colors[i]], alpha=0.3)
    ax[i].set_xlabel('Semana')
    ax[i].set_xlabel('Porcentaje de votos')
    ax[i].set_title(labelsPartidos[i])
pl.tight_layout()
pl.savefig("porPartido.png")

"""
Explanation: A continuación represento la evolución temporal del voto por cada partido político, incluyendo las incertidumbres estimadas para cada sondeo.
End of explanation
"""


f, ax = pl.subplots(ncols=3, nrows=2, figsize=(12,10))
ax = ax.flatten()
for i in range(5):
    sns.boxplot(data=houseChain[:,:,i], ax=ax[i])
    ax[i].set_xticks(np.arange(12))
    ax[i].set_xticklabels(labelsEmpresas, rotation=90)
    ax[i].set_title(labelsPartidos[i])
pl.tight_layout()
pl.savefig("bias.png")

"""
Explanation: Sesgos por casa encuestadora
Represento a continuación las estimaciones de sesgo para cada una de las empresas encuestadoras, lo que muestra algún efecto curioso en Demoscopia Servicios que sería interesante analizar y entender por qué pasa. Nótese que un sesgo positivo implica que esa empresa subestima los votos y hace falta corregirlo por esa cantidad para obtener las observaciones.
End of explanation
"""


f, ax = pl.subplots(figsize=(10,6))
sns.boxplot(data=sigmaChain, ax=ax)
ax.set_xticks(np.arange(5))
ax.set_xticklabels(labelsPartidos, rotation=90)
boxes = ax.artists
for (i, box) in enumerate(boxes):
    box.set_facecolor(sns.xkcd_rgb[colors[i]])
pl.savefig("variabilidad.png", facecolor=f.get_facecolor(), edgecolor='none')

"""
Explanation: Variabilidad
Por último, aquí muestro la variabilidad intrínseca de la estimación subyacente de voto para cada partido. Es lógico que Podemos y Cs tengan una variabilidad mayor por su largo crecimiento en el último año.
End of explanation
"""


f = open('elecciones2015/index.html','r')
lines = f.readlines()
f.close()

fecha = datetime.datetime.now()
lines[29] = """<a id="welcome-to-github-pages" class="anchor" href="#welcome-to-github-pages" aria-hidden="true"><span aria-hidden="true" class="octicon octicon-link"></span></a>Resultado {0}/{1}/{2}</h3>\n""".format(fecha.day,fecha.month,fecha.year)

theta = np.percentile(thetaChain, [50.0, 50.0-95/2.0, 50.0+95/2], axis=0)

lines[33] = """<td style='height: 50px; background: #02B3F4; color: white; font-size: 25px; border-radius: 1px 10px 0px 0px;'>{0:4.1f}-{1:4.1f}-{2:4.1f}</td>\n""".format(theta[1,-1,0]*100,theta[0,-1,0]*100,theta[2,-1,0]*100)
lines[34] = """<td style='height: 50px; background: #8E1744; color: white; font-size: 25px; border-radius: 1px 10px 0px 0px;'>{0:4.1f}-{1:4.1f}-{2:4.1f}</td>\n""".format(theta[1,-1,2]*100,theta[0,-1,2]*100,theta[2,-1,2]*100)
lines[35] = """<td style='height: 50px; background: #FF0202; color: white; font-size: 25px; border-radius: 1px 10px 0px 0px;'>{0:4.1f}-{1:4.1f}-{2:4.1f}</td>\n""".format(theta[1,-1,1]*100,theta[0,-1,1]*100,theta[2,-1,1]*100)
lines[36] = """<td style='height: 50px; background: #FF800E; color: white; font-size: 25px; border-radius: 1px 10px 0px 0px;'>{0:4.1f}-{1:4.1f}-{2:4.1f}</td>\n""".format(theta[1,-1,3]*100,theta[0,-1,3]*100,theta[2,-1,3]*100)

f = open('elecciones2015/index.html','w')
f.writelines(lines)
f.close()

os.system('./update.sh')

"""
Explanation: Push de la página web
End of explanation
"""


