#############Packages###########
import pymc

import scipy.stats as st
import pickle
import sys
import operator
#import pykov 
import random
import scipy.integrate as spi
import numpy as np
import pylab as pl
import math as mt
import scipy.linalg as linear
from scipy.linalg import leslie
import decimal
from collections import defaultdict
from matplotlib import pyplot as plt
from matplotlib.pyplot import *
from xlrd import *
from operator import truediv


###########################################
######## Definitions##############
def is_empty(any_structure):
    if any_structure:
        return False
    else:
        return True
					
					


				
###################################################################
class NestedDict(dict):
   def __getitem__(self, key):
        if key in self: return self.get(key)
        return self.setdefault(key, NestedDict())
def extract(DictIn, Dictout):
	for key, value in DictIn.iteritems():
		if isinstance(value, dict): # If value itself is dictionary
			extract(value, Dictout)
		elif isinstance(value, list): # If value itself is list
			for i in value:
				extract(i, Dictout)
		else:
			Dictout[key] = value   

#############################################################
########################################################
scale='Year' 
################################################################
def null(A, eps=1e-15):
    u, s, vh = np.linalg.svd(A)
    null_space = np.compress(s <= eps, vh, axis=0)
   
    return s
################################################################
def scalechange(annualmort,scale):
	mu=0
	if annualmort>1:
		annualmort=1
	if scale=='Year':
		mu=annualmort
	elif scale=='Month':
		mu=1-pow((1-annualmort),float(1.0)/12)
	elif scale=='Week':
		mu=1-pow((1-annualmort),float(1.0)/52)
	else:
		mu=1-pow((1-annualmort),float(1.0)/365)
	return mu		
#################################################################
def fertility(annualfert,scale):
	b=0
	if scale=='Year':
		b=annualfert
	elif scale=='Month':
		b=-1+pow((1+annualfert),float(1.0)/12)
	elif scale=='Week':
		b=-1+pow((1+annualfert),float(1.0)/52)
	else:
		b=-1+pow((1+annualfert),float(1.0)/365)
	
	return b		
#################################################################	
def carriage(a,t,scale):
	carr=0
	expo=0.455
	rateuni=0.885
	rateexp=0.645
	if a==0:	
		carr+=0.1*rateuni+0.1*(float(mt.exp(-rateexp*pow(0.5,expo)))+float(mt.exp(-rateexp*pow(1.0,expo))))
		for m in range(1,5):
			carr+=0.2*(float(mt.exp(-rateexp*pow((m)+0.5,expo))))
	else:
			carr+=(float(mt.exp(-rateexp*pow(a*10+5,expo))))
	if carr<0:
		carr=0		
	return float(carr)
###########################################
def logPoisson(Alive,observed,simulation):
	log=0
	for i in range(0, len(observed)):
		log+=0.01*Alive[i][110]*observed[i]*mt.log(0.01*Alive[i][110]*simulation[i])-mt.log(0.01*Alive[i][110]*simulation[i])
	return log
#########################################
# def PK(initialcondition,age,sero,infec,startyear,endyear,scale,fa,fc,popage,infant,forceaverage,latency,mu0,epsilon,alpha,T,b,trans,gammaA,gammaC,trans_decomp,trans_toHCC,trans_DecompCD,trans_HCCHBVD,trans_F,trans_FD,checktime,hbeag):
# 	Total=np.zeros(T)
# 	children=[]
# 	infantmod=[]
# 	m=9####number of groups
# 	nu=np.array([1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])###birth rate
# 	
# 	lambdaI=np.array([forceaverage[0],forceaverage[1],forceaverage[2],forceaverage[3],forceaverage[4],forceaverage[5],forceaverage[6],forceaverage[7],forceaverage[8]])
# 	p=np.array([carriage(0,0,scale),carriage(1,0,scale),carriage(2,0,scale),carriage(3,0,scale),carriage(4,0,scale),carriage(5,0,scale),carriage(6,0,scale),carriage(7,0,scale),carriage(8,0,scale)])
# 	fa=0.7111  # by acute mother 0.85# becoming carrier by perinatal infection
# 	fc=0.109
# 
# 	n=[]
# 	for l in range(0,m):
# 		n.append(popage[startyear][l])
# 	n=np.array(n)
# 	
# 	S0=[]
# 	E0=[]
# 	I0=[]
# 	C0=[]
# 	R0=[]
# 	Pop0=[]
# 	for l in range(0,m):
# 		S0.append(initialcondition[l][0]*n[l])
# 		E0.append(initialcondition[l][1]*n[l])
# 		I0.append(initialcondition[l][2]*n[l])
# 		C0.append(initialcondition[l][3]*n[l])
# 		R0.append(initialcondition[l][4]*n[l])
# 		Pop0.append(n[l])
# 	S0=np.array(S0)#####suscpetible
# 	E0=np.array(E0)
# 	I0=np.array(I0)
# 	C0= np.array(C0)
# 	R0= np.array(R0)
# 	Pop0= np.array(Pop0)
# 	ND=MaxTime=1
# 
# 	gammaA=gammaA
# 	gammaC=gammaC
# 	latency=latency
# 	year=startyear
# 	TS=1.0
# 
# 
# 	INPUT=np.hstack((S0,E0,I0,C0,R0,Pop0))
# 	mu=np.array([mu0[year][0],mu0[year][1],mu0[year][2],mu0[year][3],mu0[year][4],mu0[year][5],mu0[year][6],mu0[year][7],mu0[year][8]])####moratlity rate
# 	
# 	fert=np.array([0,b[1][0],b[2][0],b[3][0],b[4][0],b[5][0],b[6][0],b[7][0],b[8][0]])#([0,0,0,0,0,0,0,0,0])#b[1][0],b[2][0],b[3][0],b[4][0],b[5][0],b[6][0],b[7][0],b[8][0]])
# 	def diff_eqs(INP,t):  
# 		'''The main set of equations'''
# 		Y=np.zeros(54)
# 		V = INP
# 		Pop=[]
# 		for i in range(0,m):
# 			Pop.append(V[5*m+i])
# 		Pop=np.array(Pop)
# 		
# 		child=sum(fert[i]*Pop[i] for i in range(0,m))
# 		childA=sum(fert[i]*(fa*Y[(2*m+i)]+fc*Y[(3*m+i)]) for i in range(0,m))
# 		childS=sum(fert[i]*(Pop[i]-(fa*Y[(2*m+i)]+fc*Y[(3*m+i)])) for i in range(0,m))
# 		if child !=0:
# 			childA=infant[year+1]/child*childA
# 			childS=infant[year+1]/child*childS
# 
# 		for i in range(m):
# 			Inf=lambdaI[i]*V[i]
# 			mi=-mt.expm1(-mu[i])#1-mt.exp(-mu[i])#mortality
# 			Y[i] = nu[i]*childS  -Inf- mi * V[i]###Susceptible
# 			Y[(m+i)] =nu[i]*childA+ Inf - mi * V[(m+i)] - latency * V[(m+i)]####Latent
# 			Y[(2*m+i)] = latency * V[(m+i)] - gammaA * V[(2*m+i)] - mi * V[(2*m+i)]#####Acute
# 			Y[(3*m+i)] = gammaA *p[i]* V[(2*m+i)] - mi * V[(3*m+i)]-gammaC * V[(3*m+i)]#####Carrier
# 			Y[(4*m+i)] = gammaA *(1.0-p[i])* V[(2*m+i)]+gammaC * V[(3*m+i)] -mi * V[(4*m+i)]#####Recovered
# 			Y[(5*m+i)] = nu[i]*(childS+childA)- mi*(V[i]+ V[(m+i)]+ V[(2*m+i)]+ V[(3*m+i)]+ V[(4*m+i)])
# 		return Y   # For odeint
# 
# 	t_start = 0.0; t_end = ND; t_inc = TS
# 	t_range = np.arange(t_start, t_end+t_inc, t_inc)
# 
# 	RES2=np.zeros(54)
# 	k=1
# 	while k<=140:
# 		RES = spi.odeint(diff_eqs,INPUT,t_range)
# 		INPUT=RES[-1]
# 		year=startyear+k
# 		mu=np.array([mu0[year][0],mu0[year][1],mu0[year][2],mu0[year][3],mu0[year][4],mu0[year][5],mu0[year][6],mu0[year][7],mu0[year][8]])####moratlity rate
# 
# 		if k<len(b[1]):
# 			fert=np.array([0,b[1][k],b[2][k],b[3][k],b[4][k],b[5][k],b[6][k],b[7][k],b[8][k]])
# 		else:
# 			fert=fert	
# 		for l in range(0,6):
# 			INPUT[8+l*m]=INPUT[8+l*m]+INPUT[7+l*m]*0.1
# 			INPUT[7+l*m]=INPUT[7+l*m]+INPUT[6+l*m]*0.1-INPUT[7+l*m]*0.1
# 			INPUT[6+l*m]=INPUT[6+l*m]+INPUT[5+l*m]*0.1-INPUT[6+l*m]*0.1
# 			INPUT[5+l*m]=INPUT[5+l*m]+INPUT[4+l*m]*0.1-INPUT[5+l*m]*0.1
# 			INPUT[4+l*m]=INPUT[4+l*m]+INPUT[3+l*m]*0.1-INPUT[4+l*m]*0.1
# 			INPUT[3+l*m]=INPUT[3+l*m]+INPUT[2+l*m]*0.1-INPUT[3+l*m]*0.1
# 			INPUT[2+l*m]=INPUT[2+l*m]+INPUT[1+l*m]*0.1-INPUT[2+l*m]*0.1
# 			INPUT[1+l*m]=INPUT[1+l*m]+INPUT[0+l*m]*0.1-INPUT[1+l*m]*0.1
# 			INPUT[0+l*m]=INPUT[0+l*m]-INPUT[0+l*m]*0.1
# 		RES2=np.vstack((RES2,RES))
# 		k=k+1
# 
# 	RES=RES2[1:,]
# 	
# #################################################################	
# 	Susceptible=defaultdict(list)
# 	Latent=defaultdict(list)
# 	AcuteInfected=defaultdict(list)
# 	C0=defaultdict(list)
# 	Recovered=defaultdict(list)
# 	Alive=defaultdict(list)
# 
# 	for a in range(0,9):
# 		for t in range(0,141):
# 			Susceptible[a].append(RES[t,0+a])
# 			Latent[a].append(RES[t,1*m+a])
# 			AcuteInfected[a].append(RES[t,2*m+a])
# 			C0[a].append(RES[t,3*m+a])
# 			Recovered[a].append(RES[t,4*m+a])
# 			Alive[a].append(RES[t,5*m+a])
# 
# 
# 
# #################################################################################################
# 
# 	
# 
# 
# 	
# 
# 	return Susceptible,Latent,AcuteInfected,Recovered,C0,Alive

#######################################
color_dict = ( 'blue',  'green', 'purple', '#ff69b4', 'black',  'cyan', 'magenta', '#2e8b57',    '#da70d6',  '#ff7f50', '#cd853f', '#bc8f8f', '#5f9ea0', '#daa520')

startyear=1950
endyear=2090
checkyear=1984
Agerange=10	
TS=1.0
# if scale=='Year':
# 	ND=(endyear-startyear+1)
# elif scale=='Month':
# 	ND=(endyear-startyear+1)*12
# elif scale=='Week':
# 	ND=(endyear-startyear+1)*52
# else:
# 	ND=(endyear-startyear+1)*365



#############################################
ND=(endyear-startyear+1)
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)

if scale=='Year':
	checktime=(checkyear-startyear)
elif scale=='Month':
	checktime=12*(checkyear-startyear)	
elif scale=='Week':
	checktime=52*(checkyear-startyear)	
else:
	checktime=365*(checkyear-startyear)
############################################################################
##########################################################################################
filein='../SouthKorea/estimates_medium.txt'
IN=open(filein,'r')


popage=NestedDict()
infant=NestedDict()
Pop=NestedDict()
for y in range(startyear, endyear+1):
	Pop[int((y-startyear))]=0
	for a in range(0,100):
		if a<80:
			agegroup=int(a)/Agerange
			popage[y][a]=0
		else:
			popage[y][int(80)/Agerange]=0	
		


counter=0
for line in IN:
	if counter>0:
		elem0= line.split('\t')
		year=int(elem0[0])
		age=int(elem0[1])
		popf=int(elem0[2])
		popm=int(elem0[3])
		
		
		if year>=startyear and year<=endyear:
			Pop[int((year-startyear))]+=popf+popm
			if age<80:
				agegroup=int(age)/Agerange
				popage[year][agegroup]+=popf+popm	
				if age==0:
					infant[year]=popf+popm

			else:
				popage[year][int(80)/Agerange]+=popf+popm	

	counter+=1
###############################################################################################	


############################################################
hbeag=NestedDict()
filein='../SouthKorea/perinatal.txt'
IN=open(filein,'r')
for line in IN:	
	elem0= line.split('\t')
	age=int(elem0[0])/Agerange
	rate=float(elem0[1])
	hbeag[age]=rate
IN.close()


fertilityvec=defaultdict(list)
fertilityvec1984=defaultdict(list)
fert=NestedDict()

ratetot=NestedDict()
filein='../SouthKorea/fertility_10.txt'
IN=open(filein,'r')
counter=0
for line in IN:
	if counter>0:
		elem0= line.split('\t')
		year=int(elem0[0])
		agegroup=int(elem0[1])
		rate=float(elem0[2])
		if year>=startyear and year<endyear:
				if year in fert.keys():
					if agegroup in fert[year].keys():
						fert[year][agegroup]+=rate*infant[year+1]
					
					else:
						fert[year][agegroup]=rate*infant[year+1]
					ratetot[year]+=rate	
				else:
					
					fert[year][agegroup]=rate*infant[year+1]
					ratetot[year]=rate
					
			
	counter+=1
IN.close()
for y in range(startyear,endyear+1):
	

	if scale=='Year':
		for a in fert[y].keys():

			f=(fertility(float(fert[y][a]/(ratetot[y]*popage[y][a])),scale))
			
			fertilityvec[a].append(f)
	elif scale=='Month':
		for a in fert[y].keys():
			for i in range(0,12):
				f=(fertility(float(fert[y][a]/(ratetot[y]*popage[y][a])),scale))
				fertilityvec[a].append(f)
	elif scale=='Week':
		for a in fert[y].keys():
			for i in range(0,52):
				f=(fertility(float(fert[y][a]/(ratetot[y]*popage[y][a])),scale))
				#f=float(fert[y][a]/(52*ratetot[y]*popage[y][a]))
				fertilityvec[a].append(f)	
	else:
		for a in fert[y].keys():
			for i in range(0,365):
				f=(fertility(float(fert[y][a]/(ratetot[y]*popage[y][a])),scale))
				fertilityvec[a].append(f)
for y in range(startyear,endyear+1):
	if scale=='Year':
		for a in fert[y].keys():

			f=(fertility(float(fert[1984][a]/(ratetot[1984]*popage[1984][a])),scale))
			
			fertilityvec1984[a].append(f)
	elif scale=='Month':
		for a in fert[y].keys():
			for i in range(0,12):
				f=(fertility(float(fert[1984][a]/(ratetot[1984]*popage[1984][a])),scale))
				fertilityvec1984[a].append(f)
	elif scale=='Week':
		for a in fert[y].keys():
			for i in range(0,52):
				f=(fertility(float(fert[1984][a]/(ratetot[1984]*popage[1984][a])),scale))
				#f=float(fert[y][a]/(52*ratetot[y]*popage[y][a]))
				fertilityvec1984[a].append(f)	
	else:
		for a in fert[y].keys():
			for i in range(0,365):
				f=(fertility(float(fert[1984][a]/(ratetot[1984]*popage[1984][a])),scale))
				fertilityvec1984[a].append(f)				
######################################################################################
mortality=NestedDict()
mortality1984=NestedDict()
epsilon=NestedDict()
filein='../SouthKorea/mortality_11.csv'
IN=open(filein,'r')
counter=0
for line in IN:
	if counter>0:
		elem0= line.split(',')
		year=int(elem0[0])
		if year==1984:
			for y in range(startyear,endyear+1):
				for i in range(1,len(elem0)):
	
					mortality1984[y][i-1]=float(scalechange(float(elem0[i]),scale))
		for i in range(1,len(elem0)):
	
			mortality[year][i-1]=float(scalechange(float(elem0[i]),scale))
				

	counter+=1		




#########################		

############################################################################################################
forceaverage=np.zeros(9)
forceup=np.zeros(9)
forcedown=np.zeros(9)
filein='./SouthKorea/ForceInfection/forceinfection_vert.csv'
IN=open(filein,'r')
counter=0
for line in IN:
	if counter>0:
		elem0= line.split(',')
		if elem0[0].find('force') != -1:
			namevec=elem0[0].split('_')
			if int(namevec[1])==0:
				namevec[1]=int(namevec[1])
				
				forceaverage[int(10*namevec[1])/10]=float(elem0[1])
				forceup[int(10*namevec[1])/10]=float(elem0[-1])
				forcedown[int(10*namevec[1])/10]=float(elem0[-5])
				

			elif int(namevec[1])==6:
				for k in range(6,9):
					forceaverage[k]=float(elem0[1])
					forceup[k]=float(elem0[-1])
					forcedown[k]=float(elem0[-5])
			else:
				namevec[1]=int(namevec[1])
				
				forceaverage[int(10*namevec[1])/10]=float(elem0[1])
				forceup[int(10*namevec[1])/10]=float(elem0[-1])
				forcedown[int(10*namevec[1])/10]=float(elem0[-5])
						
	counter+=1



#forceaverage=([0.00434915116967,0.00390892223653,0.00319164506509,0.00262948231521, 0.00221458943872,0.00190431828876,0.00166634405814, 0.00147917629932,0.00132864300251,0.00120521736764,0.00110233053921,0.00101533625127,0.000940869217625,0.0008764390259,  0.000820166839682,0.000770610564981, 0.000770610564981 ])
age0=26
age=(5,15,25,35,45,55,65)
sero=(25.93194896,44.34151885,66.7939514,67.55553954,80.6335684,81.39473249,74.4361237)
infec=(6.20,7.06,8.71,8.77,8.23,5.82,4.37)

###################################################
initialcondition=defaultdict(list)	

for a in range(0,9):

	year0=startyear
	if a<6:	
		
		initialcondition[a].append(1.0-0.01*sero[a]-0.01*infec[a])
		initialcondition[a].append(0.01*infec[a])
		initialcondition[a].append(0.005*infec[a])
		initialcondition[a].append(0.005*infec[a])
		
		initialcondition[a].append(0.01*sero[a]-0.01*infec[a])
	else:
		initialcondition[a].append(1.0-0.01*sero[-1]-0.01*infec[-1])
		initialcondition[a].append(0.01*infec[-1])
		initialcondition[a].append(0.005*infec[-1])
		initialcondition[a].append(0.005*infec[-1])
		
		initialcondition[a].append(0.01*sero[-1]-0.01*infec[-1])
Total=0
pa=[]
for a in range(0,9):

	year0=startyear
	Total+=popage[year0][a]
#########################################################################
latency=pymc.Uniform("latency",0.05,0.75)
gammaA=pymc.Uniform("gammaA",0.05,0.75)
gammaC=pymc.Uniform("gammaC",0.005,0.05)
#gammaC=pymc.Uniform("gammaC",0.001,1.0)



fa=0.7111  # by acute mother 0.85# becoming carrier by perinatal infection
fc=0.109 #by carrier mother

alpha=0.16

#latencyannual=float(lat)/12######annuallatency
transannual=0.0734547001096 ## transition between carriage states
#gammaCannual=0.0125#0.01 Chronic  recovery  rate### transition from carrier state to immune
#gammaAannual=float(g)/12 #Acute recovery rate in John paper is 4 gammaCannual
trans_decompannual=0.056### transition rate from cirrhosis to decompensated 
trans_HCCCarrierannual=0.016634## transition rate to carrier from HCC 
trans_toHCCannual=0.03##from cirrhosis  decompensate to HCC 
trans_DecompCarrierannual=0.043733## transition rate to carrier from Decompensated
trans_Fannual=0.0183####From Acute to Fulminnat Liver Disease

trans_DecompCDannual=0.25## from decompensated to death for HBV related cause
trans_HCCHBVDannual=0.56##from Cancer to Death for HBV related cause
trans_FDannual=0.0046#0.25136612#########From Fulminat to Death

trans=float(scalechange(transannual,scale))
#gammaA=scalechange(gammaAannual,scale)
#gammaC=float(scalechange(gammaCannual,scale))
#latency=float(scalechange(latencyannual,scale))
trans_decomp=float(scalechange(trans_decompannual,scale))
trans_HCCCarrier=float(scalechange(trans_HCCCarrierannual,scale))
trans_toHCC=float(scalechange(trans_toHCCannual,scale))
trans_DecompCarrier=float(scalechange(trans_DecompCarrierannual,scale))
trans_DecompCD=float(scalechange(trans_DecompCDannual,scale))
trans_HCCHBVD=float(scalechange(trans_HCCHBVDannual,scale))
trans_F=float(scalechange(trans_Fannual,scale))
trans_FD=float(scalechange(trans_FDannual,scale))
####################################################
@pymc.deterministic
def PK2(initialcondition=initialcondition,age=age,sero=sero,infec=infec,startyear=startyear,endyear=endyear,scale=scale,fa=fa,fc=fc,popage=popage,infant=infant,forceaverage=forceaverage,latency=latency,mu0=mortality1984,T=ND,b=fertilityvec1984,gammaA=gammaA,gammaC=gammaC):	
	Total=np.zeros(T)
	children=[]
	infantmod=[]
	m=9####number of groups
	nu=np.array([1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])###birth rate
	
	lambdaI=np.array([forceaverage[0],forceaverage[1],forceaverage[2],forceaverage[3],forceaverage[4],forceaverage[5],forceaverage[6],forceaverage[7],forceaverage[8]])
	p=np.array([carriage(0,0,scale),carriage(1,0,scale),carriage(2,0,scale),carriage(3,0,scale),carriage(4,0,scale),carriage(5,0,scale),carriage(6,0,scale),carriage(7,0,scale),carriage(8,0,scale)])
	fa=0.7111  # by acute mother 0.85# becoming carrier by perinatal infection
	fc=0.109

	n=[]
	for l in range(0,m):
		n.append(popage[startyear][l])
	n=np.array(n)
	
	S0=[]
	E0=[]
	I0=[]
	C0=[]
	R0=[]
	Pop0=[]
	for l in range(0,m):
		S0.append(initialcondition[l][0]*n[l])
		E0.append(initialcondition[l][1]*n[l])
		I0.append(initialcondition[l][2]*n[l])
		C0.append(initialcondition[l][3]*n[l])
		R0.append(initialcondition[l][4]*n[l])
		Pop0.append(n[l])
	S0=np.array(S0)#####suscpetible
	E0=np.array(E0)
	I0=np.array(I0)
	C0= np.array(C0)
	R0= np.array(R0)
	Pop0= np.array(Pop0)
	MaxTime=1

	gammaA=gammaA
	gammaC=gammaC
	latency=latency
	year=startyear
	TS=1.0


	INPUT=np.hstack((S0,E0,I0,C0,R0,Pop0))
	mu=np.array([mu0[year][0],mu0[year][1],mu0[year][2],mu0[year][3],mu0[year][4],mu0[year][5],mu0[year][6],mu0[year][7],mu0[year][8]])####moratlity rate
	
	fert=np.array([0,b[1][0],b[2][0],b[3][0],b[4][0],b[5][0],b[6][0],b[7][0],b[8][0]])#([0,0,0,0,0,0,0,0,0])#b[1][0],b[2][0],b[3][0],b[4][0],b[5][0],b[6][0],b[7][0],b[8][0]])
	def diff_eqs(INP,t):  
		'''The main set of equations'''
		Y=np.zeros(54)
		V = INP
		Pop=[]
		for i in range(0,m):
			Pop.append(V[5*m+i])
		Pop=np.array(Pop)
		
		child=sum(fert[i]*Pop[i] for i in range(0,m))
		childA=sum(fert[i]*(fa*Y[(2*m+i)]+fc*Y[(3*m+i)]) for i in range(0,m))
		childS=sum(fert[i]*(Pop[i]-(fa*Y[(2*m+i)]+fc*Y[(3*m+i)])) for i in range(0,m))
		if child !=0:
			childA=infant[year+1]/child*childA
			childS=infant[year+1]/child*childS

		for i in range(m):
			Inf=lambdaI[i]*V[i]
			mi=-mt.expm1(-mu[i])#1-mt.exp(-mu[i])#mortality
			Y[i] = nu[i]*childS  -Inf- mi * V[i]###Susceptible
			Y[(m+i)] =nu[i]*childA+ Inf - mi * V[(m+i)] - latency * V[(m+i)]####Latent
			Y[(2*m+i)] = latency * V[(m+i)] - gammaA * V[(2*m+i)] - mi * V[(2*m+i)]#####Acute
			Y[(3*m+i)] = gammaA *p[i]* V[(2*m+i)] - mi * V[(3*m+i)]-gammaC * V[(3*m+i)]#####Carrier
			Y[(4*m+i)] = gammaA *(1.0-p[i])* V[(2*m+i)]+gammaC * V[(3*m+i)] -mi * V[(4*m+i)]#####Recovered
			Y[(5*m+i)] = nu[i]*(childS+childA)- mi*(V[i]+ V[(m+i)]+ V[(2*m+i)]+ V[(3*m+i)]+ V[(4*m+i)])
		return Y   # For odeint

	t_start = 0.0; t_end = MaxTime; t_inc = TS
	t_range = np.arange(t_start, t_end+t_inc, t_inc)

	RES2=np.zeros(54)
	k=1
	while k<=140:
		RES = spi.odeint(diff_eqs,INPUT,t_range)
		INPUT=RES[-1]
		year=startyear+k
		mu=np.array([mu0[year][0],mu0[year][1],mu0[year][2],mu0[year][3],mu0[year][4],mu0[year][5],mu0[year][6],mu0[year][7],mu0[year][8]])####moratlity rate

		if k<len(b[1]):
			fert=np.array([0,b[1][k],b[2][k],b[3][k],b[4][k],b[5][k],b[6][k],b[7][k],b[8][k]])
		else:
			fert=fert	
		for l in range(0,6):
			INPUT[8+l*m]=INPUT[8+l*m]+INPUT[7+l*m]*0.1
			INPUT[7+l*m]=INPUT[7+l*m]+INPUT[6+l*m]*0.1-INPUT[7+l*m]*0.1
			INPUT[6+l*m]=INPUT[6+l*m]+INPUT[5+l*m]*0.1-INPUT[6+l*m]*0.1
			INPUT[5+l*m]=INPUT[5+l*m]+INPUT[4+l*m]*0.1-INPUT[5+l*m]*0.1
			INPUT[4+l*m]=INPUT[4+l*m]+INPUT[3+l*m]*0.1-INPUT[4+l*m]*0.1
			INPUT[3+l*m]=INPUT[3+l*m]+INPUT[2+l*m]*0.1-INPUT[3+l*m]*0.1
			INPUT[2+l*m]=INPUT[2+l*m]+INPUT[1+l*m]*0.1-INPUT[2+l*m]*0.1
			INPUT[1+l*m]=INPUT[1+l*m]+INPUT[0+l*m]*0.1-INPUT[1+l*m]*0.1
			INPUT[0+l*m]=INPUT[0+l*m]-INPUT[0+l*m]*0.1
		RES2=np.vstack((RES2,RES))
		k=k+1

	RES=RES2[1:,]
	
#################################################################	
	Susceptible=defaultdict(list)
	Latent=defaultdict(list)
	AcuteInfected=defaultdict(list)
	C0=defaultdict(list)
	Recovered=defaultdict(list)
	Alive=defaultdict(list)

	for a in range(0,9):
		for t in range(0,141):
			Susceptible[a].append(RES[t,0+a])
			Latent[a].append(RES[t,1*m+a])
			AcuteInfected[a].append(RES[t,2*m+a])
			C0[a].append(RES[t,3*m+a])
			Recovered[a].append(RES[t,4*m+a])
			Alive[a].append(RES[t,5*m+a])



#################################################################################################

###################################################################################################################################################
	initialcondition2=defaultdict(list)	

	for a in range(0,9):

		year0=startyear
		initialcondition2[a].append(Susceptible[a][138]/Alive[a][138])
		initialcondition2[a].append(Latent[a][138]/Alive[a][138])
		initialcondition2[a].append(AcuteInfected[a][138]/Alive[a][138])
		initialcondition2[a].append(C0[a][138]/Alive[a][138])

		initialcondition2[a].append(Recovered[a][138]/Alive[a][138])
	def mu(t,mu0):
			year=1950+int(t)
			mu=np.array([mu0[year][0],mu0[year][1],mu0[year][2],mu0[year][3],mu0[year][4],mu0[year][5],mu0[year][6],mu0[year][7],mu0[year][8]])####moratlity rate
			return mu

	t_range = np.linspace(0,138,138)	
	initialcondition=initialcondition2
	nu=np.array([1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])###birth rate
	
	lambdaI=np.array([forceaverage[0],forceaverage[1],forceaverage[2],forceaverage[3],forceaverage[4],forceaverage[5],forceaverage[6],forceaverage[7],forceaverage[8]])
	p=np.array([carriage(0,0,scale),carriage(1,0,scale),carriage(2,0,scale),carriage(3,0,scale),carriage(4,0,scale),carriage(5,0,scale),carriage(6,0,scale),carriage(7,0,scale),carriage(8,0,scale)])
	fa=0.7111  # by acute mother 0.85# becoming carrier by perinatal infection
	fc=0.109
	m=9
	n=[]
	for l in range(0,m):
		n.append(popage[startyear][l])
	n=np.array(n)

	S0=[]
	E0=[]
	I0=[]
	C0=[]
	R0=[]
	Pop0=[]
	for l in range(0,m):
		S0.append(initialcondition[l][0]*n[l])
		E0.append(initialcondition[l][1]*n[l])
		I0.append(initialcondition[l][2]*n[l])
		C0.append(initialcondition[l][3]*n[l])
		R0.append(initialcondition[l][4]*n[l])
		Pop0.append(n[l])
	S0=np.array(S0)#####suscpetible
	E0=np.array(E0)
	I0=np.array(I0)
	C0= np.array(C0)
	R0= np.array(R0)
	Pop0= np.array(Pop0)
	print "Done equilibrium"
	epsilonup=np.array([0,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1])	
	epsilond=np.array([0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0])	


	INPUT=np.hstack((S0,E0,I0,C0,R0,Pop0))
	m=9####number of groups
	mort=np.array([0.00037125846825,0.00010538894875,0.0002457375915,0.000503574828,0.00135061587525,0.00305347856725,0.0080537843695,0.027210962198,0.179291737335])

	def diff_eqs(INP,t):  
		
		'''The main set of equations'''

		Y=np.zeros(54)
		V = INP
		Pop=[]
		for i in range(0,m):
			Pop.append(V[5*m+i])
		Pop=np.array(Pop)
		year=startyear+int(t)
		
		if t<len(fertilityvec[1]):
			fert=np.array([0,fertilityvec[1][int(t)],fertilityvec[2][int(t)],fertilityvec[3][int(t)],fertilityvec[4][int(t)],fertilityvec[5][int(t)],fertilityvec[6][int(t)],fertilityvec[7][int(t)],fertilityvec[8][int(t)]])
		else:
			fert=np.array([0,fertilityvec[1][-1],fertilityvec[2][-1],fertilityvec[3][-1],fertilityvec[4][-1],fertilityvec[5][-1],fertilityvec[6][-1],fertilityvec[7][-1],fertilityvec[8][-1]])
		child=sum(fert[i]*Pop[i] for i in range(0,m))
		childA=sum(fert[i]*(fa*Y[(2*m+i)]+fc*Y[(3*m+i)]) for i in range(0,m))
		childS=sum(fert[i]*(Pop[i]-(fa*Y[(2*m+i)]+fc*Y[(3*m+i)])) for i in range(0,m))
		
		if child !=0:
			if year+1 not in infant.keys():
				print year+1
				
			childA=infant[year+1]/child*childA
			childS=infant[year+1]/child*childS
	

		for i in range(m):
			Inf=lambdaI[i]*V[i]
			mi=-mt.expm1(-mort[i])#1-mt.exp(-mu[i])#mortality
			Y[i] = nu[i]*childS  -Inf- mi * V[i]###Susceptible
			Y[(m+i)] =nu[i]*childA+ Inf - mi * V[(m+i)] - latency * V[(m+i)]####Latent
			Y[(2*m+i)] = latency * V[(m+i)] - gammaA * V[(2*m+i)] - mi * V[(2*m+i)]#####Acute
			Y[(3*m+i)] = gammaA *p[i]* V[(2*m+i)] - mi * V[(3*m+i)]-gammaC * V[(3*m+i)]#####Carrier
			Y[(4*m+i)] = gammaA *(1.0-p[i])* V[(2*m+i)]+gammaC * V[(3*m+i)] -mi * V[(4*m+i)]#####Recovered
			Y[(5*m+i)] = nu[i]*(childS+childA)- mi*(V[i]+ V[(m+i)]+ V[(2*m+i)]+ V[(3*m+i)]+ V[(4*m+i)])
		return Y   # For odeint
	Maxtime=1
	t_start = 0.0; t_end = Maxtime; t_inc = TS
	t_range = np.arange(t_start, t_end+t_inc, t_inc)

	RES2=np.zeros(54)
	k=1
	while k<=138:
		RES = spi.odeint(diff_eqs,INPUT,t_range)
		INPUT=RES[-1]
		year=startyear+k
	
		for l in range(0,6):
			INPUT[8+l*m]=INPUT[8+l*m]+INPUT[7+l*m]*0.1
			INPUT[7+l*m]=INPUT[7+l*m]+INPUT[6+l*m]*0.1-INPUT[7+l*m]*0.1
			INPUT[6+l*m]=INPUT[6+l*m]+INPUT[5+l*m]*0.1-INPUT[6+l*m]*0.1
			INPUT[5+l*m]=INPUT[5+l*m]+INPUT[4+l*m]*0.1-INPUT[5+l*m]*0.1
			INPUT[4+l*m]=INPUT[4+l*m]+INPUT[3+l*m]*0.1-INPUT[4+l*m]*0.1
			INPUT[3+l*m]=INPUT[3+l*m]+INPUT[2+l*m]*0.1-INPUT[3+l*m]*0.1
			INPUT[2+l*m]=INPUT[2+l*m]+INPUT[1+l*m]*0.1-INPUT[2+l*m]*0.1
			INPUT[1+l*m]=INPUT[1+l*m]+INPUT[0+l*m]*0.1-INPUT[1+l*m]*0.1
			INPUT[0+l*m]=INPUT[0+l*m]-INPUT[0+l*m]*0.1
		RES2=np.vstack((RES2,RES))
		k=k+1

	RES=RES2[1:,]
	
#################################################################	
	Susceptible=defaultdict(list)
	Latent=defaultdict(list)
	AcuteInfected=defaultdict(list)
	C0=defaultdict(list)
	Recovered=defaultdict(list)
	Alive=defaultdict(list)

	for a in range(0,9):
		for t in range(0,139):
			Susceptible[a].append(RES[t,0+a])
			Latent[a].append(RES[t,1*m+a])
			AcuteInfected[a].append(RES[t,2*m+a])
			C0[a].append(RES[t,3*m+a])
			Recovered[a].append(RES[t,4*m+a])
			Alive[a].append(RES[t,5*m+a])

	# for a in range(0,9):
# 		for t in range(0,141):	
# 			Alive[a].append(Susceptible[a][t]+Latent[a][t]+AcuteInfected[a][t]+C0[a][t]+Recovered[a][t])
	T=140
	Total=np.zeros(T)
	for a in sorted(Alive.keys()):
			if a !='__name__':
				for t in range(0,len(Alive[a])):
					if t !='__name__' and t<=T:
				
						Total[t]+=Alive[a][t]
	susc=np.zeros(T)
	carriers2=np.zeros(T)
	acute=np.zeros(T)
	susckeys=sorted(Susceptible.keys())
	for s in susckeys:
		if s !='__name__':
			for t in range(0,len(Susceptible[s])):
				if t<T:
					if t !='__name__': 
					
						susc[t]+=float(Susceptible[s][t])/Total[t]
						carriers2[t]+=float(C0[s][t])/Total[t] 
						acute[t]+=float(AcuteInfected[s][t])/Total[t]#
	rec=np.zeros(T)

	reckeys=sorted(Recovered.keys())
	for s in reckeys:
		if s !='__name__':
			for t in range(0,len(Recovered[s])):
				if t<T:
					if t !='__name__': 
						rec[t]+=float(Recovered[s][t])/Total[t]
	Serop1970=np.zeros(138)
	Serop1980=np.zeros(138)
	Serop1990=np.zeros(138)
	Serop2000=np.zeros(138)
	Serop2020=np.zeros(138)
	Serop2040=np.zeros(138)
	Serop2060=np.zeros(138)
	Infec1970=np.zeros(138)
	Infec1980=np.zeros(138)
	Infec1990=np.zeros(138)
	Infec2000=np.zeros(138)
	Infec2020=np.zeros(138)
	Infec2040=np.zeros(138)
	Infec2060=np.zeros(138)
	for s in range(0,7):
		if s !='__name__':
			t=20
			carriers=C0[s][t]
			Serop1970[s]=(int(AcuteInfected[s][t]+Recovered[s][t]+carriers))
			Infec1970[s]=(int(AcuteInfected[s][t]+carriers))
			t=30
			carriers=C0[s][t]
			Serop1980[s]=(int(AcuteInfected[s][t]+Recovered[s][t]+carriers))
			Infec1980[s]=(int(AcuteInfected[s][t]+carriers))
			t=40
			carriers=C0[s][t]
			Serop1990[s]=(int(AcuteInfected[s][t]+Recovered[s][t]+carriers))
			Infec1990[s]=(int(AcuteInfected[s][t]+carriers))
			t=50
			carriers=C0[s][t]
			Serop2000[s]=(int(AcuteInfected[s][t]+Recovered[s][t]+carriers))	
			Infec2000[s]=(int(AcuteInfected[s][t]+carriers))
			t=70
			carriers=C0[s][t]
			Serop2020[s]=(int(AcuteInfected[s][t]+Recovered[s][t]+carriers))
			Infec2020[s]=(int(AcuteInfected[s][t]+carriers))
			t=90
			carriers=C0[s][t]
			Serop2040[s]=(int(AcuteInfected[s][t]+Recovered[s][t]+carriers))
			Infec2040[s]=(int(AcuteInfected[s][t]+carriers))
			t=110
			carriers=C0[s][t]
			Serop2060[s]=(int(AcuteInfected[s][t]+Recovered[s][t]+carriers))								
			Infec2060[s]=(int(AcuteInfected[s][t]+carriers))	
	


	
	return (Serop1980,Infec1980)

	#return np.array(Serop1970),np.array(Serop1980),np.array(Serop1990),np.array(Serop2000),np.array(Serop2020),np.array(Serop2040),np.array(Serop2060),np.array(Infec1970),np.array(Infec1980),np.array(Infec1990),np.array(Infec2000),np.array(Infec2020),np.array(Infec2040),np.array(Infec2060)
#Serop1970=pymc.Lambda('Serop1970', lambda PK2=PK2: PK2[0])
Serop1980=pymc.Lambda('Serop1980', lambda PK2=PK2: PK2[0])
# Serop1990=pymc.Lambda('Serop1990', lambda PK2=PK2: PK2[2])
# Serop2000=pymc.Lambda('Serop2000', lambda PK2=PK2: PK2[3])
# Serop2020=pymc.Lambda('Serop2020', lambda PK2=PK2: PK2[4])
# Serop2040=pymc.Lambda('Serop2040', lambda PK2=PK2: PK2[5])
# Serop2060=pymc.Lambda('Serop2060', lambda PK2=PK2: PK2[6])
# Infec1970=pymc.Lambda('Infec1970', lambda PK2=PK2: PK2[7])
Infec1980=pymc.Lambda('Infec1980', lambda PK2=PK2: PK2[1])
# Infec1990=pymc.Lambda('Infec1990', lambda PK2=PK2: PK2[9])
# Infec2000=pymc.Lambda('Infec2000', lambda PK2=PK2: PK2[10])
# Infec2020=pymc.Lambda('Infec2020', lambda PK2=PK2: PK2[11])
# Infec2040=pymc.Lambda('Infec2040', lambda PK2=PK2: PK2[12])
# Infec2060=pymc.Lambda('Infec2060', lambda PK2=PK2: PK2[13])

print "Done Model"

# sero1970=[]
sero1980=[]
# sero1990=[]
# sero2000=[]
# sero2020=[]
# sero2040=[]
# sero2060=[]
# 
# infec1970=[]
infec1980=[]
# infec1990=[]
# infec2000=[]
# infec2020=[]
# infec2040=[]
# infec2060=[]
for a in range(0,7):
# 	sero1970.append(int(0.01*sero[a]*popage[1970][a]))
	sero1980.append(int(0.01*sero[a]*popage[1980][a]))
# 	sero1990.append(int(0.01*sero[a]*popage[1990][a]))
# 	sero2000.append(int(0.01*sero[a]*popage[2000][a]))
# 	sero2020.append(int(0.01*sero[a]*popage[2020][a]))
# 	sero2040.append(int(0.01*sero[a]*popage[2040][a]))
# 	sero2060.append(int(0.01*sero[a]*popage[2060][a]))
# 	infec1970.append(int(0.01*infec[a]*popage[1970][a]))
	infec1980.append(int(0.01*infec[a]*popage[1980][a]))
# 	infec1990.append(int(0.01*infec[a]*popage[1990][a]))
# 	infec2000.append(int(0.01*infec[a]*popage[2000][a]))
# 	infec2020.append(int(0.01*infec[a]*popage[2020][a]))
# 	infec2040.append(int(0.01*infec[a]*popage[2040][a]))
# 	infec2060.append(int(0.01*infec[a]*popage[2060][a]))



#########################################################################################################################################	
	
#################################################################################################

@pymc.deterministic
def predict(Serop1980=Serop1980):
	predict=[]
	sero=Serop1980
	
	for i in range(0, 7):
		predict.append(sero[i])
	
	predict=np.array(predict)
	return predict

print predict.value
#serobservationKorea_1970 = pymc.Poisson('serobservationKorea_1970', mu=Serop1970av,  value=sero1970,observed=True) 
serobservationKorea_1980_0 = pymc.Poisson('serobservationKorea_1980_0', mu=predict[0],value=sero1980[0],observed=True) 
serobservationKorea_1980_1 = pymc.Poisson('serobservationKorea_1980_1', mu=predict[1],value=sero1980[1],observed=True) 
serobservationKorea_1980_2 = pymc.Poisson('serobservationKorea_1980_2', mu=predict[2],value=sero1980[2],observed=True) 
serobservationKorea_1980_3 = pymc.Poisson('serobservationKorea_1980_3', mu=predict[3],value=sero1980[3],observed=True) 
serobservationKorea_1980_4 = pymc.Poisson('serobservationKorea_1980_4', mu=predict[4],value=sero1980[4],observed=True) 
serobservationKorea_1980_5 = pymc.Poisson('serobservationKorea_1980_5', mu=predict[5],value=sero1980[5],observed=True) 
serobservationKorea_1980_6 = pymc.Poisson('serobservationKorea_1980_6', mu=predict[6],value=sero1980[6],observed=True) 
#serobservationKorea_1990 = pymc.Poisson('serobservationKorea_1990', mu=Serop1990av,   value=sero1990,observed=True) 
# serobservationKorea_2000 = pymc.Poisson('serobservationKorea_2000', mu=Serop2000av,   value=sero2000,observed=True) 
# serobservationKorea_2020 = pymc.Poisson('serobservationKorea_2020', mu=Serop2020av,   value=sero2020,observed=True) 
# serobservationKorea_2040 = pymc.Poisson('serobservationKorea_2040', mu=Serop2040av,   value=sero2040,observed=True) 
# serobservationKorea_2060 = pymc.Poisson('serobservationKorea_2060', mu=Serop2060av,   value=sero2060,observed=True) 

# surfobservationKorea_1970 = pymc.Poisson('surfobservationKorea_1970', mu=Infec1970av,  value=infec1970,observed=True) 
# surfobservationKorea_1980 = pymc.Poisson('surfobservationKorea_1980', mu=Infec1980,value=infec1980,observed=True) 
#surfobservationKorea_1990 = pymc.Poisson('surfobservationKorea_1990', mu=Infec1990av,   value=infec1990,observed=True) 
# surfobservationKorea_2000 = pymc.Poisson('surfobservationKorea_2000', mu=Infec2000av,   value=infec2000,observed=True) 
# surfobservationKorea_2020 = pymc.Poisson('surfobservationKorea_2020', mu=Infec2020av,   value=infec2020,observed=True) 
# surfobservationKorea_2040 = pymc.Poisson('surfobservationKorea_2040', mu=Infec2040av,   value=infec2040,observed=True) 
# surfobservationKorea_2060 = pymc.Poisson('surfobservationKorea_2060', mu=Infec2060av,   value=infec2060,observed=True) 



 
# pobservationKorea_1970 = pymc.Poisson('pobservationKorea_1970', mu=Serop1970) 
pobservationKorea_1980_0 = pymc.Poisson('pobservationKorea_1980_0', mu=predict[0]) 
pobservationKorea_1980_1 = pymc.Poisson('pobservationKorea_1980_1', mu=predict[1]) 
pobservationKorea_1980_2 = pymc.Poisson('pobservationKorea_1980_2', mu=predict[2]) 
pobservationKorea_1980_3 = pymc.Poisson('pobservationKorea_1980_3', mu=predict[3]) 
pobservationKorea_1980_4 = pymc.Poisson('pobservationKorea_1980_4', mu=predict[4]) 
pobservationKorea_1980_5 = pymc.Poisson('pobservationKorea_1980_5', mu=predict[5]) 
pobservationKorea_1980_6 = pymc.Poisson('pobservationKorea_1980_6', mu=predict[6]) 
# pobservationKorea_1990 = pymc.Poisson('pobservationKorea_1990', mu=Serop1990) 
# pobservationKorea_2000 = pymc.Poisson('pobservationKorea_2000', mu=Serop2000) 
# pobservationKorea_2020 = pymc.Poisson('pobservationKorea_2020', mu=Serop2020) 
# pobservationKorea_2040 = pymc.Poisson('pobservationKorea_2040', mu=Serop2040) 
# pobservationKorea_2060 = pymc.Poisson('pobservationKorea_2060', mu=Serop2060) 

##################################################################################################################



# 
# infecobservationKorea_1970 = pymc.Poisson('infecobservationKorea_1970', mu=Infec1970) 
#infecobservationKorea_1980 = pymc.Poisson('infecobservationKorea_1980', mu=Infec1980) 
# infecobservationKorea_1990 = pymc.Poisson('infecobservationKorea_1990', mu=Infec1990) 
# infecobservationKorea_2000 = pymc.Poisson('infecobservationKorea_2000', mu=Infec2000) 
# infecobservationKorea_2020 = pymc.Poisson('infecobservationKorea_2020', mu=Infec2020) 
# infecobservationKorea_2040 = pymc.Poisson('infecobservationKorea_2040', mu=Infec2040) 
# infecobservationKorea_2060 = pymc.Poisson('infecobservationKorea_2060', mu=Infec2060) 
# 	
# 		for i in range(0,m):
# 			Inf=lambdaI[i]*V[i]
# 			mi=-mt.expm1(-mort[i])#1-mt.exp(-mu[i])#mortality
# 			if i<1:
# 				Y[i] = nu[i]*childS  -Inf -mi * V[i]-epsilond[i]*V[i]###Susceptible
# 				Y[(m+i)] =nu[i]*childA+ Inf - mi * V[(m+i)] - latency * V[(m+i)]-epsilond[i]*V[m+i]####Latent
# 				Y[(2*m+i)] = latency * V[(m+i)] - gammaA * V[(2*m+i)] - mi * V[(2*m+i)]-epsilond[i]*V[2*m+i]#####Acute
# 				Y[(3*m+i)] = gammaA *p[i]* V[(2*m+i)] - mi * V[(3*m+i)]-gammaC * V[(3*m+i)]-epsilond[i]*V[3*m+i]#####Carrier
# 				Y[(4*m+i)] = gammaA *(1.0-p[i])* V[(2*m+i)]+gammaC * V[(3*m+i)] -mi * V[(4*m+i)]-epsilond[i]*V[4*m+i]#####Recovered
# 				Y[(5*m+i)] = nu[i]*(childS+childA)- mi*(V[(5*m+i)])-epsilond[i]*V[5*m+i]
# 			elif i>=1 and i<8:	
# 				Y[i] = nu[i]*childS  -Inf -mi * V[i]-epsilond[i]*V[i]+epsilonup[i]*V[i-1]###Susceptible
# 				Y[(m+i)] =nu[i]*childA+ Inf - mi * V[(m+i)] - latency * V[(m+i)]-epsilond[i]*V[m+i]+epsilonup[i]*V[m+i-1]####Latent
# 				Y[(2*m+i)] = latency * V[(m+i)] - gammaA * V[(2*m+i)] - mi * V[(2*m+i)]-epsilond[i]*V[2*m+i]+epsilonup[i]*V[2*m+i-1]#####Acute
# 				Y[(3*m+i)] = gammaA *p[i]* V[(2*m+i)] - mi * V[(3*m+i)]-gammaC * V[(3*m+i)]-epsilond[i]*V[3*m+i]+epsilonup[i]*V[3*m+i-1]#####Carrier
# 				Y[(4*m+i)] = gammaA *(1.0-p[i])* V[(2*m+i)]+gammaC * V[(3*m+i)] -mi * V[(4*m+i)]-epsilond[i]*V[4*m+i]+epsilonup[i]*V[4*m+i-1]#####Recovered
# 				Y[(5*m+i)] = nu[i]*(childS+childA)- mi*(V[(5*m+i)])-epsilond[i]*V[5*m+i]+epsilonup[i]*V[5*m+i-1]
# 			else:
# 				Y[i] = nu[i]*childS  -Inf -mi * V[i]+epsilonup[i]*V[i-1]###Susceptible
# 				Y[(m+i)] =nu[i]*childA+ Inf - mi * V[(m+i)] - latency * V[(m+i)]+epsilonup[i]*V[m+i-1]####Latent
# 				Y[(2*m+i)] = latency * V[(m+i)] - gammaA * V[(2*m+i)] - mi * V[(2*m+i)]+epsilonup[i]*V[2*m+i-1]#####Acute
# 				Y[(3*m+i)] = gammaA *p[i]* V[(2*m+i)] - mi * V[(3*m+i)]-gammaC * V[(3*m+i)]+epsilonup[i]*V[3*m+i-1]#####Carrier
# 				Y[(4*m+i)] = gammaA *(1.0-p[i])* V[(2*m+i)]+gammaC * V[(3*m+i)] -mi * V[(4*m+i)]+epsilonup[i]*V[4*m+i-1]#####Recovered
# 				Y[(5*m+i)] = nu[i]*(childS+childA)- mi*(V[(5*m+i)])+epsilonup[i]*V[5*m+i-1]
# 					
# 		return Y   # For odeint
# 
# 	RES = spi.odeint(diff_eqs,INPUT,t_range)