import matplotlib
#matplotlib.use('Agg')
from pylab import *
from pymc import *
from operator import truediv 
from matplotlib import pyplot as plt
from pymc.Matplot import plot 
# load the model
import MCMC_model_vertical11 as mod #.pyforceinfectionColombia_postv 
reload(mod) # this reload streamlines interactive debugging
 


mc=pymc.MCMC(mod)#,db='pickle',dbname='../Colombia/Outputs/ResultsSimulation/NewForceColombiadb')
#mc.fit()
#mc.use_step_method(pymc.AdaptiveMetropolis, [mod.gammaA, mod.latency,mod.gammaC])
mc.sample(iter=20, burn=10, thin=1, verbose=1)
#
mc.write_csv("./SouthKorea/MCMCForceInfection/Vertical/ParametersKorea.csv" )

#


xKorea=[]



for i in range(0,len(mod.age)):
	if i!='__name__':

		xKorea.append(mod.age[i])		
############################################################################
# averagesero1970Korea=(mod.pobservationKorea_1970.stats()['mean'])[0:7]
# sero1970Korea25=(mod.pobservationKorea_1970.stats()['quantiles'][2.5])[0:7]	
# sero1970Korea97=(mod.pobservationKorea_1970.stats()['quantiles'][97.5])[0:7]

averagesero1980Korea_0=(mod.pobservationKorea_1980_0.stats()['mean'])
sero1980Korea25_0=(mod.pobservationKorea_1980_0.stats()['quantiles'][2.5])
sero1980Korea97_0=(mod.pobservationKorea_1980_0.stats()['quantiles'][97.5])
averagesero1980Korea_1=(mod.pobservationKorea_1980_1.stats()['mean'])
sero1980Korea25_1=(mod.pobservationKorea_1980_1.stats()['quantiles'][2.5])
sero1980Korea97_1=(mod.pobservationKorea_1980_1.stats()['quantiles'][97.5])
averagesero1980Korea_2=(mod.pobservationKorea_1980_2.stats()['mean'])
sero1980Korea25_2=(mod.pobservationKorea_1980_2.stats()['quantiles'][2.5])
sero1980Korea97_2=(mod.pobservationKorea_1980_2.stats()['quantiles'][97.5])
averagesero1980Korea_3=(mod.pobservationKorea_1980_3.stats()['mean'])
sero1980Korea25_3=(mod.pobservationKorea_1980_3.stats()['quantiles'][2.5])
sero1980Korea97_3=(mod.pobservationKorea_1980_3.stats()['quantiles'][97.5])
averagesero1980Korea_4=(mod.pobservationKorea_1980_4.stats()['mean'])
sero1980Korea25_4=(mod.pobservationKorea_1980_4.stats()['quantiles'][2.5])
sero1980Korea97_4=(mod.pobservationKorea_1980_4.stats()['quantiles'][97.5])
averagesero1980Korea_5=(mod.pobservationKorea_1980_5.stats()['mean'])
sero1980Korea25_5=(mod.pobservationKorea_1980_5.stats()['quantiles'][2.5])
sero1980Korea97_5=(mod.pobservationKorea_1980_5.stats()['quantiles'][97.5])
averagesero1980Korea_6=(mod.pobservationKorea_1980_6.stats()['mean'])
sero1980Korea25_6=(mod.pobservationKorea_1980_6.stats()['quantiles'][2.5])
sero1980Korea97_6=(mod.pobservationKorea_1980_6.stats()['quantiles'][97.5])
averagesero1980Korea=([averagesero1980Korea_0,averagesero1980Korea_1,averagesero1980Korea_2,averagesero1980Korea_3,averagesero1980Korea_4,averagesero1980Korea_5,averagesero1980Korea_6])
sero1980Korea25=([sero1980Korea25_0,sero1980Korea25_1,sero1980Korea25_2,sero1980Korea25_3,sero1980Korea25_4,sero1980Korea25_5,sero1980Korea25_6])
sero1980Korea97=([sero1980Korea97_0,sero1980Korea97_1,sero1980Korea97_2,sero1980Korea97_3,sero1980Korea97_4,sero1980Korea97_5,sero1980Korea97_6])

# averagesero1990Korea=(mod.pobservationKorea_1990.stats()['mean'])[0:7]
# sero1990Korea25=(mod.pobservationKorea_1990.stats()['quantiles'][2.5])[0:7]	
# sero1990Korea97=(mod.pobservationKorea_1990.stats()['quantiles'][97.5])[0:7]
# 
# averagesero2000Korea=(mod.pobservationKorea_2000.stats()['mean'])[0:7]
# sero2000Korea25=(mod.pobservationKorea_2000.stats()['quantiles'][2.5])[0:7]	
# sero2000Korea97=(mod.pobservationKorea_2000.stats()['quantiles'][97.5])[0:7]
# 
# averagesero2020Korea=(mod.pobservationKorea_2020.stats()['mean'])[0:7]
# sero2020Korea25=(mod.pobservationKorea_2020.stats()['quantiles'][2.5])[0:7]	
# sero2020Korea97=(mod.pobservationKorea_2020.stats()['quantiles'][97.5])[0:7]
# 
# averagesero2040Korea=(mod.pobservationKorea_2040.stats()['mean'])[0:7]
# sero2040Korea25=(mod.pobservationKorea_2040.stats()['quantiles'][2.5])[0:7]	
# sero2040Korea97=(mod.pobservationKorea_2040.stats()['quantiles'][97.5])[0:7]

ageSero=[10*i+5 for i in range(0,7)]		
figure(2)

# plt.plot(xKorea,mod.sero1970, color='black', linestyle='None', marker='o')
plt.plot(xKorea,mod.sero1980, color='black', linestyle='None', marker='o')
# plt.plot(xKorea,mod.sero1990, color='black', linestyle='None', marker='o')
# plt.plot(xKorea,mod.sero2000, color='black', linestyle='None', marker='o')
# plt.plot(xKorea,mod.sero2020, color='black', linestyle='None', marker='o')
# plt.plot(xKorea,mod.sero2040, color='black', linestyle='None', marker='o')
# plt.plot(xKorea,mod.sero2060, color='black', linestyle='None', marker='o')

# plt.plot(ageSero,averagesero1970Korea, color='red',label='1970')
# 
# plt.fill_between(ageSero,sero1970Korea25,sero1970Korea97, facecolor='r',edgecolor='r',alpha=0.1)

plt.plot(ageSero,averagesero1980Korea, color='green',label='1980 ')

plt.fill_between(ageSero,sero1980Korea25,sero1980Korea97, facecolor='g',edgecolor='g',alpha=0.1)

# plt.plot(ageSero,averagesero1990Korea, color='blue',label='1990')
# 
# plt.fill_between(ageSero,sero1990Korea25,sero1990Korea97, facecolor='b',edgecolor='b',alpha=0.1)
# 
# plt.plot(ageSero,averagesero2000Korea, color='black',label='2000')
# 
# plt.fill_between(ageSero,sero2000Korea25,sero2000Korea97, facecolor='k',edgecolor='k',alpha=0.1)
# 
# plt.plot(ageSero,averagesero2020Korea, color='orange',label='2020')
# 
# plt.fill_between(ageSero,sero2020Korea25,sero2020Korea97, facecolor='orange',edgecolor='orange',alpha=0.1)
# 
# plt.plot(ageSero,averagesero2040Korea, color='magenta',label='2040')
# 
# plt.fill_between(ageSero,sero2040Korea25,sero2040Korea97, facecolor='magenta',edgecolor='magenta',alpha=0.1)



plt.legend(loc=4)
ylabel('Seroprevalence')
xlabel('Age Group')
savefig('./SouthKorea/MCMCForceInfection/Vertical/Seroprevalence.png')
plt.close()
print "Done SeroPrevalence"
#########################################################################