<jupyter_start><jupyter_text>conda install seaborn<jupyter_code>pop = pd.read_csv('Populacoes.csv')
#The state of the data that will be fitted
s='BA'

# Total Population
#N=9557071

# Extract population 

p = pop[pop.Sigla == s].reset_index()

N = p.Pop[0]
BA_data=[]
with open('AcompanhamentoCorona.txt', encoding="utf8") as data:                                                                                          
    data_reader = csv.reader(data, delimiter='\t')
    for data in data_reader:
        BA_data.append(data)


datas=[]
ts0=[]
infec =[]
leitos=[]
Uti= []
dth= []
        
        
for i in range(len(BA_data)):
    datas.append(BA_data[i][0])
    ts0.append(int(BA_data[i][1]))
    infec.append(int(BA_data[i][2]))
    leitos.append(int(BA_data[i][3]))
    Uti.append(int(BA_data[i][4]))
    dth.append(int(BA_data[i][5]))

    
ts0=np.array(ts0)
infec=np.array(infec)
leitos=np.array(leitos)
Uti=np.array(Uti)
dth=np.array(dth)
    
dth_cm=[sum(dth[:i]) for i in range(len(dth))]

dth_cm=np.array(dth_cm)

# Total Population
#N=14873064


#Parameters That Are Fixed

#delta = 1/2

kappa = 1/4

p = 0.2

gammaA = 1/3.5

gammaS = 1/4

#gammaH = 1/10

#gammaU = 1/8
 
muH = 0.15

muU = 0.4 

#h =  0.12

xi = 0.53 

omega_U = 0.29 
omega_H = 0.14



#Find when the death starts

for i in range(len(dth_cm)):
    if dth_cm[i] != 0:
        ti=i
        break
#ti = date of the start of the death curve


#Defining the Steap Function 
k=50
def H(t):
    h = 1.0/(1.0+ np.exp(-2.0*k*t))
    return h


#Defining the Beta(t)
def beta(t,t1,b,b1):
    beta = b*H(t1-t) + b1*H(t-t1) 
    return beta


#Defining the Model
def SEIRHU(f,t,parametros):
    
    #parameters
    b, b1,gammaH,gammaU, delta, h, t1 = parametros
    
    #variables
    S = f[0]
    E = f[1]
    IA = f[2]
    IS = f[3]
    H = f[4]
    U = f[5]
    R = f[6]
    D = f[7]
    Nw = f[8]
    
    
    #equations
    dS_dt = - beta(t,t1,b,b1)*S*(IS + delta*IA) 
    dE_dt = beta(t,t1,b,b1)*S*(IS + delta*IA) - kappa * E
    dIA_dt = (1-p)*kappa*E - gammaA*IA
    dIS_dt = p*kappa*E - gammaS*IS
    dH_dt = h*xi*gammaS*IS + (1- muU+ omega_U*muU)*gammaU*U - gammaH*H                 #(1-muU)*gammaU*U -gammaH*H
    dU_dt = h*(1-xi)*gammaS*IS + omega_H*gammaH*H -gammaU*U                 #+(1-(muH))*omega*gammaH*H -gammaU*U
    dR_dt = gammaA*IA + (1-(muH))*(1-omega_H)*gammaH*H + (1-h)*gammaS*IS
    dD_dt = (1-omega_H)*muH*gammaH*H + (1-omega_U)*muU*gammaU*U
    
    
    #epidemic curve
    dNw_dt = p*kappa*E
    
    #Returning the derivatives
    return [dS_dt,dE_dt,dIA_dt,dIS_dt,dH_dt,dU_dt,dR_dt,dD_dt,dNw_dt]


#Calculating the Error
def lq_SEIRHU(pars, ts0):
    
    #assaing the Least square given parameters to the numerical integration for the error calculation
    b, b1, gammaH,gammaU, delta, h, ia0, is0, e0, t1 = pars
    
    #initial conditions
    q0 = [1-ia0 - is0 -e0,e0,ia0,is0,0,0,0,0,0]
    
    #parameters to feed the E.D.Os
    parode = b, b1,gammaH,gammaU, delta, h, t1
    
    #Integrating
    qs = odeint(SEIRHU,q0, ts0,args=(parode,),mxstep=1000000)

    #sinf the epidemic curve
    #sdth the death curve

    sinf = qs[:,-1]
    sdth = qs[:,-2]
    slet = qs[:,4]
    suti = qs[:,5]
    
    #define the standardized residuals
    erri = (N*sinf - infec) / np.sqrt(N*sinf+1)
    errd = (N*sdth - dth_cm) / np.sqrt(N*sdth+1)
    errh = (N*slet - leitos) / np.sqrt(N*sinf+1)
    erru = (N*suti - Uti) / np.sqrt(N*sdth+1)
    
    return np.r_[erri,errd] 



#define your data


#assing an array of zeros to the days before the first death



#create an array for the time series starting in t=1
#we need it to start in 1 because we need to integrate from 0-1 to define the number of cases of the first day.
#This helps us because we can set the numer of new cases in t=0 equal 0, than there is no need to fit this parameter. 
ts0 = np.array(list(range(len(infec))))+1

 

#DEEFINE THE INTERVALS OF PARS b, b1, xi, delta, h, ia0, is0, e0, t1
intervals = np.array(   [ [0., 1.5],#b
                          [0., 1.],#b1
                          [1/14,1/5],#gammaH
                          [1/14, 1/5],#gammaU
                          [0., 1.],#delta
                          [0.05, 0.35],#h
                          [0.,10./N],#ia0
                          [0.,10./N],#is0
                          [0.,10./N],#e0
                          [0.,30]#t1
                          ])

#NUMBER OF FITS
n_tries = 10

#best error starts at inft 
best_err = np.inf
#aj the number of best fists 
j = 0

for i in range(n_tries):
    
    #create a set of ramdom parameters
    par0 = np.random.rand(len(intervals))
    #Limit those parameters to the interval defined
    par0 = intervals[:,0] + par0 * (intervals[:,1] - intervals[:,0])
    
    #try to fit

    res = least_squares(lambda pars: lq_SEIRHU(pars, ts0), par0, bounds=(intervals[:,0], \
                                                 intervals[:,1]))
    ier = res.status
    #if converges
    if ier >= 1 and ier <= 4:
        if res.cost < best_err:
            #best_err = erro
            best_err = res.cost
            #best_pop = parameters
            best_pop = res.x
            #+1 best fit
            j = j +1
            print('achou',j)



#put best pop in these parameters  
b, b1, gammaH,gammaU, delta, h, ia0, is0, e0, t1 = best_pop
parode =b, b1,gammaH,gammaU, delta, h, t1

q0 = [1-ia0 - is0 -e0,e0,ia0,is0,0,0,0,0,0]

#plot the data
qs = odeint(SEIRHU,q0, ts0,args=(parode,),mxstep=1000000)
pl.figure()
qq =qs[:,-1]
qd =qs[:,-2]
qh= qs[:,4]
qu= qs[:,5]

initial_date=date(int(2020), int(2), int(28))

dates = []
dates.append(initial_date.strftime('%d-%m-%Y'))  

for i in range(len(ts0)-1):
    d=initial_date + timedelta(days=i)
    dates.append(d.strftime('%d-%m-%Y'))   


    


# Create figure and plot space
fig, ((ax1, ax2), (ax3, ax4)) = pl.subplots(2, 2,figsize=(13, 10))



# Add x-axis and y-axis

ax1.plot(dates, N*qq,'b-', linewidth=2.5,label='Casos Previstos ')
ax2.plot(dates[30:], N*qd[30:],'r-', linewidth=2.5,label='Mortes Previstas')
ax3.plot(dates[30:], N*qh[30:],'g-', linewidth=2.5,label='Leitos de Enfermaria Previstos ')
ax4.plot(dates[30:], N*qu[30:],color='darkviolet', linewidth=2.5,label='Leitos de UTI Previstas ')

ax1.plot(dates, infec, 'ko',label='Casos Reportados ')
ax2.plot(dates[30:], dth_cm[30:], 'ko',label='Mortes Reportadas ')
ax3.plot(dates[30:], leitos[30:], 'ko',label='Leitos Enfermaria Ocupados ')
ax4.plot(dates[30:], Uti[30:], 'ko',label='Leitos UTI Ocupados ')

ax1.axvline(t1, 0, 600,c='r',linestyle='--',label='Intervenção')


    
# Set title and labels for axes
fig.suptitle("Epidemia Covid-19 " + '(' + s + ')',fontsize=20)


ax1.set_xlabel('Data',fontsize=17)
ax1.set_ylabel('Casos Confirmados Acumulados',fontsize=17)
ax1.legend(loc='upper left',fontsize=15)

ax2.set_xlabel('Data',fontsize=17)
ax2.set_ylabel('Mortes',fontsize=17)
ax2.legend(loc='upper left',fontsize=15)

ax3.set_xlabel('Data',fontsize=17)
ax3.set_ylabel('Leitos de Enfermaria',fontsize=17)
ax3.legend(loc='upper left',fontsize=15)

ax4.set_xlabel('Data',fontsize=17)
ax4.set_ylabel('Leitos de UTI',fontsize=17)
ax4.legend(loc='upper left',fontsize=15)


# Rotate tick marks on x-axis
ax1.xaxis.set_major_locator(pl.MaxNLocator(10))
pl.setp(ax1.get_xticklabels(), rotation=45)

ax2.xaxis.set_major_locator(pl.MaxNLocator(10))
pl.setp(ax2.get_xticklabels(), rotation=45)

ax3.xaxis.set_major_locator(pl.MaxNLocator(10))
pl.setp(ax3.get_xticklabels(), rotation=45)

ax4.xaxis.set_major_locator(pl.MaxNLocator(10))
pl.setp(ax4.get_xticklabels(), rotation=45)






pl.show()
best_err
#fit(self, x, y, d, pesoMorte=0.5, kappa=0.25, p=0.2, gammaA=0.2857142857142857, gammaS=0.25, muH=0.15, muU=0.4, xi=0.53, omegaU=0.29, omegaH=0.14, bound=[[0, 0.125, 0.08333333333333333, 0, 0.05], [2, 0.25, 0.3333333333333333, 0.7, 0.25]], stand_error=False, beta2=True, day_mudar=None, paramPSO={'options': {'c1': 0.1, 'c2': 0.3, 'w': 0.9, 'k': 5, 'p': 2}, 'n_particles': 300, 'iter': 1000})
modelo =md.SEIRHUD(N)
#bound beta1 beta2 t gammaH gammaU delta h ia0 is0 e0
bound = ([0,0,0,1/14,1/14,0,0.05,0,0,0],[1.5,1,30,1/5,1/5,1,0.35,10/N,10/N,10/N])
modelo.fit(range(1,len(infec)+1),infec,dth_cm,pesoMorte = 0.5,bound=bound,kappa=kappa,p=p,gammaA=gammaA,gammaS=gammaS,muH=muH,muU=muU,xi=xi,omegaU=omega_U,omegaH=omega_H,stand_error=True)
def rmse(ypre,yreal):
    yp = np.array(ypre)
    yr = np.array(yreal)
    
    return (((yp-yr)/ np.sqrt(yp+1))**2).mean()
modelo.predict(range(1,len(infec)+1))
rmse_I_pso = rmse(modelo.ypred,infec)
rmse_I = rmse(N*qq,infec)
res = {"rmse_infectados":(rmse_I,rmse_I_pso)}
rmse_D_pso = rmse(modelo.dpred,dth_cm)
rmse_D = rmse(N*qd,dth_cm)
res['rmse_mortes'] = (rmse_D,rmse_D_pso)
rmse_H_pso = rmse(modelo.H,leitos)
rmse_H = rmse(N*qh,leitos)
res['rmse_leitos'] = (rmse_H,rmse_H_pso)
rmse_U_pso = rmse(modelo.U,Uti)
rmse_U = rmse(N*qu,Uti)
res['rmse_UTI'] = (rmse_U,rmse_U_pso)

ban = pd.DataFrame(res,index=['Trust_Region(algoritmo usado por Davi e Juliane)','PSO'])
ban
modelo.getCoef()
banco = pd.DataFrame({'data':dates,'casos':modelo.y,'casos_pred':modelo.ypred,'mortes':dth_cm,'mortes_pred':modelo.dpred,'leitos':leitos,'leitos_pred':modelo.H,'UTI':Uti,'UTI_pred':modelo.U})
banco
pl.figure()
fig, ((ax1, ax2), (ax3, ax4)) = pl.subplots(2, 2,figsize=(13, 10))
# Add x-axis and y-axis

ax1.plot(dates, modelo.ypred,'b-', linewidth=2.5,label='Casos Previstos ')
ax2.plot(dates[30:], modelo.dpred[30:],'r-', linewidth=2.5,label='Mortes Previstas')
ax3.plot(dates[30:], modelo.H[30:],'g-', linewidth=2.5,label='Leitos de Enfermaria Previstos ')
ax4.plot(dates[30:], modelo.U[30:],color='darkviolet', linewidth=2.5,label='Leitos de UTI Previstas ')

ax1.plot(dates, infec, 'ko',label='Casos Reportados ')
ax2.plot(dates[30:], dth_cm[30:], 'ko',label='Mortes Reportadas ')
ax3.plot(dates[30:], leitos[30:], 'ko',label='Leitos Enfermaria Ocupados ')
ax4.plot(dates[30:], Uti[30:], 'ko',label='Leitos UTI Ocupados ')

ax1.axvline(modelo.day_mudar, 0, 600,c='r',linestyle='--',label='Intervenção')


    
# Set title and labels for axes
fig.suptitle("Epidemia Covid-19 " + '(' + s + ')',fontsize=20)


ax1.set_xlabel('Data',fontsize=17)
ax1.set_ylabel('Casos Confirmados Acumulados',fontsize=17)
ax1.legend(loc='upper left',fontsize=15)

ax2.set_xlabel('Data',fontsize=17)
ax2.set_ylabel('Mortes',fontsize=17)
ax2.legend(loc='upper left',fontsize=15)

ax3.set_xlabel('Data',fontsize=17)
ax3.set_ylabel('Leitos de Enfermaria',fontsize=17)
ax3.legend(loc='upper left',fontsize=15)

ax4.set_xlabel('Data',fontsize=17)
ax4.set_ylabel('Leitos de UTI',fontsize=17)
ax4.legend(loc='upper left',fontsize=15)


# Rotate tick marks on x-axis
ax1.xaxis.set_major_locator(pl.MaxNLocator(10))
pl.setp(ax1.get_xticklabels(), rotation=45)

ax2.xaxis.set_major_locator(pl.MaxNLocator(10))
pl.setp(ax2.get_xticklabels(), rotation=45)

ax3.xaxis.set_major_locator(pl.MaxNLocator(10))
pl.setp(ax3.get_xticklabels(), rotation=45)

ax4.xaxis.set_major_locator(pl.MaxNLocator(10))
pl.setp(ax4.get_xticklabels(), rotation=45)






pl.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17676767676767677, 'fraction_numerical': 0.05844155844155844, 'mean_word_length': 3.569015245158632, 'pattern_counts': {'":': 1, '<': 7, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '437868', 'license': 'MIT', 'n_tokens_mistral': 5209, 'n_tokens_neox': 4877, 'n_words': 869, 'path': 'Reproducibility of published results/Evaluating the burden of COVID-19 in Bahia, Brazil: A modeling analysis of 14.8 million individuals/script/.ipynb_checkpoints/para_ontem-checkpoint.ipynb', 'repo_name': 'lhunlindeion/Mathematical-and-Statistical-Modeling-of-COVID19-in-Brazil', 'size': '11088'}