#nbi:hide_in
def RMPlot(Country,Type,name2,time_fit):
name = Country
# New form of data acquisition. It pulls updated data from the github repository used by JHU
url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv'
df_recovered = pandas.read_csv(url, error_bad_lines=False)
url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv'
df_deaths = pandas.read_csv(url, error_bad_lines=False)
url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv'
df_confirmed = pandas.read_csv(url, error_bad_lines=False)
####### Italy data #######
df_deaths = df_deaths.groupby(['Country/Region']).aggregate(np.sum).reset_index()
df_confirmed = df_confirmed.groupby(['Country/Region']).aggregate(np.sum).reset_index()
#print(df_deaths['Country/Region'].values.tolist())
#print(df_confirmed['Country/Region'].values.tolist())
deaths = df_deaths[(df_deaths['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float')[0][np.where(df_deaths[(df_deaths['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float') > 0)[1][0]:len(df_deaths.columns.values)]
recovered = df_recovered[(df_recovered['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float')[0][np.where(df_deaths[(df_deaths['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float') > 0)[1][0]:len(df_deaths.columns.values)]
confirmed = df_confirmed[(df_confirmed['Country/Region'] == name)][df_confirmed.columns.values[3:len(df_confirmed.columns.values)]].to_numpy(dtype='float')[0][np.where(df_confirmed[(df_confirmed['Country/Region'] == name)][df_confirmed.columns.values[3:len(df_confirmed.columns.values)]].to_numpy(dtype='float') > 0)[1][0]:len(df_confirmed.columns.values)]
if Type == 'Deaths':
kind1 = deaths
if Type == 'Cases':
kind1 = confirmed
kind1 = kind1[0:time_fit]
t_forecast = 1 + 0/100
t = np.linspace(0, int(t_forecast*len(kind1)), 2000)
t0 = np.array(range(len(kind1)))
#print(globals()[str(name)+'_deaths'])
# fcn2min é uma das funções mais importantes do pacote lmfit. Basicamente ela define o modelo que o minimizador vai utilizar para fitar
# os pontos dos dados
def expfit(x,C0,r):
return C0*np.exp(r*x)
t_exp = np.array(range(len(kind1[0:20])))
ydata = kind1[0:20]
popt, pcov = curve_fit(expfit, t_exp, ydata)
C0 = popt[0]
#print(C0)
def fcn2min(params, t, data):
"""Generalized Richard's Model"""
K = params['K']
alpha = params['alpha']
r = params['r']
#tc = params['tc']
tc = np.log((-1+(K/C0)**alpha)/alpha)/(alpha*r)
model = K/(1+alpha*np.exp(-alpha*r*(t-tc)))**(1/alpha)
return model - data
# A próxima etaa é definir os chutes iniciais. Eu sempre tento fixando os limites ou deixando sem limite
##### Parâmetros livres quarantine after day 30
params = Parameters()
params.add('K', value=np.amax(kind1), min=0)
params.add('alpha', value=.1, min=0, max=1)
params.add('r', value=0.2, min=0, max=1)
#params.add('tc', value=10, min=0)
start =timer()
# create Minimizer
minner = Minimizer(fcn2min, params, fcn_args=(t0,kind1))
# first solve with Nelder-Mead algorithm
out1 = minner.minimize(method= 'least_squares')
end = timer()
#print('time spent first minimization:', end - start)
#lmfit.report_fit(out1)
start = timer()
out2 = minner.minimize(method='leastsq', params=out1.params)
end = timer()
#print('time spent second minimization:', end - start)
#lmfit.report_fit(out2)
R2 = 1 - out2.residual.var() / np.var(kind1)
global R2_generalc
R2_generalc = R2
params = out2.params
K = params['K']
alpha = params['alpha']
r = params['r']
#tc = params['tc']
K = np.float(params['K'].value)
alpha = np.float(params['alpha'].value)
r = np.float(params['r'].value)
tc = np.log((-1+(K/C0)**alpha)/alpha)/(alpha*r)
tc = np.float(tc)
tjp = tc - np.log(1.5 + 0.5*alpha - 0.5*np.sqrt(5 + 6*alpha + alpha**2))/(r*alpha)
tjm = tc - np.log(1.5 + 0.5*alpha + 0.5*np.sqrt(5 + 6*alpha + alpha**2))/(r*alpha)
ts = np.real(( alpha )**( -1 ) * ( r )**( -1 ) * np.log( ( 1/3 * ( 7 + 4 * alpha ) * ( np.e )**( alpha * r * tc ) + ( -1/3 * ( 2 )**( 1/3 ) * ( -31 * ( np.e )**( 2 * alpha * r * tc ) + ( -44 * alpha * ( np.e )**( 2 * alpha * r * tc ) + -13 * ( alpha )**( 2 ) * ( np.e )**( 2 * alpha * r * tc ) ) ) * ( ( 335 * ( np.e )**( 3 * alpha * r * tc ) + ( 708 * alpha * ( np.e )**( 3 * alpha * r * tc ) + ( 465 * ( alpha )**( 2 ) * ( np.e )**( 3 * alpha * r * tc ) + ( 92 * ( alpha )**( 3 ) * ( np.e )**( 3 * alpha * r * tc ) + 3 * ( 3 )**( 1/2 ) * ( ( -257 * ( np.e )**( 6 * alpha * r * tc ) + ( -1224 * alpha * ( np.e )**( 6 * alpha * r * tc ) + ( -2122 * ( alpha )**( 2 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -1712 * ( alpha )**( 3 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -681 * ( alpha )**( 4 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -136 * ( alpha )**( 5 ) * ( np.e )**( 6 * alpha * r * tc ) + -12 * ( alpha )**( 6 ) * ( np.e )**( 6 * alpha * r * tc ) ) ) ) ) ) ) )**( 1/2 ) ) ) ) ) )**( -1/3 ) + 1/3 * ( 2 )**( -1/3 ) * ( ( 335 * ( np.e )**( 3 * alpha * r * tc ) + ( 708 * alpha * ( np.e )**( 3 * alpha * r * tc ) + ( 465 * ( alpha )**( 2 ) * ( np.e )**( 3 * alpha * r * tc ) + ( 92 * ( alpha )**( 3 ) * ( np.e )**( 3 * alpha * r * tc ) + 3 * ( 3 )**( 1/2 ) * ( ( -257 * ( np.e )**( 6 * alpha * r * tc ) + ( -1224 * alpha * ( np.e )**( 6 * alpha * r * tc ) + ( -2122 * ( alpha )**( 2 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -1712 * ( alpha )**( 3 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -681 * ( alpha )**( 4 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -136 * ( alpha )**( 5 ) * ( np.e )**( 6 * alpha * r * tc ) + -12 * ( alpha )**( 6 ) * ( np.e )**( 6 * alpha * r * tc ) ) ) ) ) ) ) )**( 1/2 ) ) ) ) ) )**( 1/3 ) ) ) ))
K = params['K']
alpha = params['alpha']
r = params['r']
global stage
if t[-1]>tjm and t[-1]<=tc:
stage = 'early intermediate'
if t[-1]>tc and t[-1]<=tjp:
stage = 'late intermediate'
if t[-1]>tjp and t[-1]<=ts:
stage = 'early saturation'
if t[-1]>ts:
stage = 'late saturation'
def func(t,r,alpha,K):
tc = np.log((-1+(K/C0)**alpha)/alpha)/(alpha*r)
return K/(1+alpha*np.exp(-alpha*r*(t-tc)))**(1/alpha)
return Country, Type, name2, t0, kind1, K, r, alpha, C0
def RMPlotPlot(Country, Type, name2, t0, kind1, K, r, alpha, C0, Forecast, interv_time, interv_param, togg1, togg2, x_axis, y_axis, toggdetails1):
name = Country
#t_forecast = 1 + Forecast/100
t = np.arange(0, len(kind1)+Forecast, 1)
def func(t,r,alpha,K):
tc = np.log((-1+(K/C0)**alpha)/alpha)/(alpha*r)
return K/(1+alpha*np.exp(-alpha*r*(t-tc)))**(1/alpha)
alpha1 = alpha
K1 = K
r1 = r
K = np.float(K.value)
alpha = np.float(alpha.value)
r = np.float(r.value)
tc = np.log((-1+(K/C0)**alpha)/alpha)/(alpha*r)
tc = np.float(tc)
tjp = tc - np.log(1.5 + 0.5*alpha - 0.5*np.sqrt(5 + 6*alpha + alpha**2))/(r*alpha)
tjm = tc - np.log(1.5 + 0.5*alpha + 0.5*np.sqrt(5 + 6*alpha + alpha**2))/(r*alpha)
ts = np.real(( alpha )**( -1 ) * ( r )**( -1 ) * np.log( ( 1/3 * ( 7 + 4 * alpha ) * ( np.e )**( alpha * r * tc ) + ( -1/3 * ( 2 )**( 1/3 ) * ( -31 * ( np.e )**( 2 * alpha * r * tc ) + ( -44 * alpha * ( np.e )**( 2 * alpha * r * tc ) + -13 * ( alpha )**( 2 ) * ( np.e )**( 2 * alpha * r * tc ) ) ) * ( ( 335 * ( np.e )**( 3 * alpha * r * tc ) + ( 708 * alpha * ( np.e )**( 3 * alpha * r * tc ) + ( 465 * ( alpha )**( 2 ) * ( np.e )**( 3 * alpha * r * tc ) + ( 92 * ( alpha )**( 3 ) * ( np.e )**( 3 * alpha * r * tc ) + 3 * ( 3 )**( 1/2 ) * ( ( -257 * ( np.e )**( 6 * alpha * r * tc ) + ( -1224 * alpha * ( np.e )**( 6 * alpha * r * tc ) + ( -2122 * ( alpha )**( 2 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -1712 * ( alpha )**( 3 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -681 * ( alpha )**( 4 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -136 * ( alpha )**( 5 ) * ( np.e )**( 6 * alpha * r * tc ) + -12 * ( alpha )**( 6 ) * ( np.e )**( 6 * alpha * r * tc ) ) ) ) ) ) ) )**( 1/2 ) ) ) ) ) )**( -1/3 ) + 1/3 * ( 2 )**( -1/3 ) * ( ( 335 * ( np.e )**( 3 * alpha * r * tc ) + ( 708 * alpha * ( np.e )**( 3 * alpha * r * tc ) + ( 465 * ( alpha )**( 2 ) * ( np.e )**( 3 * alpha * r * tc ) + ( 92 * ( alpha )**( 3 ) * ( np.e )**( 3 * alpha * r * tc ) + 3 * ( 3 )**( 1/2 ) * ( ( -257 * ( np.e )**( 6 * alpha * r * tc ) + ( -1224 * alpha * ( np.e )**( 6 * alpha * r * tc ) + ( -2122 * ( alpha )**( 2 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -1712 * ( alpha )**( 3 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -681 * ( alpha )**( 4 ) * ( np.e )**( 6 * alpha * r * tc ) + ( -136 * ( alpha )**( 5 ) * ( np.e )**( 6 * alpha * r * tc ) + -12 * ( alpha )**( 6 ) * ( np.e )**( 6 * alpha * r * tc ) ) ) ) ) ) ) )**( 1/2 ) ) ) ) ) )**( 1/3 ) ) ) ))
global date_first
if Type == 'Deaths':
date_first_peak1 = date_first + timedelta(days=int(tc))
if Type == 'Cases':
date_first_peak1 = date_first + timedelta(days=int(tc))
lang = html_widget_lang.value
if lang == 'Por':
date_first1 = date_first.strftime('%d/%m/%y')
date_first_peak1 = date_first_peak1.strftime('%d/%m/%y')
if lang == 'Eng':
date_first1 = date_first.strftime('%m/%d/%y')
date_first_peak1 = date_first_peak1.strftime('%m/%d/%y')
alpha = alpha1
K = K1
r = r1
def funcprime(t,r,alpha,K):
tc = np.log((-1+(K/C0)**alpha)/alpha)/(alpha*r)
return (K*r*alpha*np.exp(-alpha*r*(t-tc)))/(1+alpha*np.exp(-alpha*r*(t-tc)))**(1 + 1/alpha)
if Type=='Deaths':
tipo = 'Mortes'
tipo1 = 'morte'
tipo2 = 'a primeira morte'
if Type=='Cases':
tipo = 'Casos'
tipo1 = 'caso'
tipo2 = 'o primeiro caso'
curve_forecast = func(t,r,alpha,K)
curve = func(np.arange(0, len(kind1), 1),r,alpha,K)
increment = curve_forecast[-1] - curve[-1]
fig, ax = plt.subplots(figsize=(6,4))
ax.plot(t0,kind1,'ro', alpha=0.5, label=str(subs[lang][Type])+' '+str(name))
ax.plot(t,func(t,r,alpha,K), '-k', label=methods_name1[lang][1])
if togg1 == True:
if interv_param>=0:
ax.plot(t1,func1(t1,r1,alpha1,K1), '-b', label='Intervention')
if interv_param<0:
ax.plot(t1,func1(t1,r1,alpha1,K1), '-r', label='Relaxation')
ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=r'$\gamma$= '+str('%.2f' % eta))
ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=first_death[lang][Type]+str(date_first1))
#ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=first_plateau[lang]+str(int(K.value)))
if Forecast>0:
ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=forecast_message[lang][str(Forecast)]+'+'+str(increment)+' '+subs[lang][Type])
ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=forecast_message[lang][str(Forecast)]+str(curve_forecast[-1])+' '+subs[lang][Type])
#ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=r'$T_c= %.1f$' %(tc))
ax.plot(t[0],kind1[0], 'ko', alpha=0, label=estagio[lang])
global stage
if stage == 'early intermediate':
ax.plot(np.ones(100)*tjm,np.linspace(0,func(tjm,r,alpha,K),100), '-', alpha=1, color='darkorange')
ax.plot(t[0],kind1[0], 'o', alpha=0, label=fase[lang][1], color='darkorange')
if stage == 'late intermediate':
ax.plot(np.ones(100)*tjm,np.linspace(0,func(tjm,r,alpha,K),100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,func(tc,r,alpha,K),100), '-', alpha=1, color='yellow')
ax.plot(t[0],kind1[0], 'o', alpha=0, label=fase[lang][2], color='yellow')
if stage == 'early saturation':
ax.plot(np.ones(100)*tjm,np.linspace(0,func(tjm,r,alpha,K),100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,func(tc,r,alpha,K),100), '-', alpha=1, color='yellow')
ax.plot(np.ones(100)*tjp,np.linspace(0,func(tjp,r,alpha,K),100), '-', alpha=1, color='green')
ax.plot(t[0],kind1[0], 'o', alpha=0, label=fase[lang][3], color='green')
if stage == 'late saturation':
ax.plot(np.ones(100)*tjm,np.linspace(0,func(tjm,r,alpha,K),100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,func(tc,r,alpha,K),100), '-', alpha=1, color='yellow')
ax.plot(np.ones(100)*tjp,np.linspace(0,func(tjp,r,alpha,K),100), '-', alpha=1, color='green')
ax.plot(np.ones(100)*ts,np.linspace(0,func(ts,r,alpha,K),100), '-', alpha=1, color='royalblue')
ax.plot(t[0],kind1[0], 'o', alpha=0, label=fase[lang][4], color='royalblue')
ax.plot(t[0],kind1[0], 'o', alpha=0, label=r'$\copyright$2021 ModInterv')
ax.set_xlabel(label_x[lang][Type])
ax.set_ylabel(label_y[lang][Type])
ax.set_title(name+' '+methods_name[lang][1],y=1.05)
plt.figtext(0.51,0.9,time_dic[lang]+' UTC-3', fontsize=8,ha='center')
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True)
plt.xscale(dic[x_axis])
plt.yscale(dic[y_axis])
if x_axis=='Linear' and y_axis=='Linear':
legend = ax.legend(loc='best', fontsize='small')
line = legend.get_lines()
text = legend.get_texts()
text[len(text)-2].set_color(line[len(text)-2].get_color())
if x_axis=='Logarithmic':
ax.set_xlim(left=0.8)
if y_axis=='Logarithmic':
ax.set_ylim(bottom=0.8)
for spine in ('top', 'right', 'bottom', 'left'):
ax.spines[spine].set_visible(False)
plt.savefig('fig.png', format='png', dpi=200)
ax.set_rasterized(True)
plt.savefig('fig.eps', format='eps')
plt.show()
if toggdetails1 == True:
fig, ax = plt.subplots(figsize=(6,4))
ax.plot(t0,kind1,'ro', alpha=0.5, label=str(subs[lang][Type])+' '+str(name))
ax.plot(t,func(t,r,alpha,K), '-k', label=methods_name1[lang][1])
if togg1 == True:
if interv_param>=0:
ax.plot(t1,func1(t1,r1,alpha1,K1), '-b', label='Intervention')
if interv_param<0:
ax.plot(t1,func1(t1,r1,alpha1,K1), '-r', label='Relaxation')
ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=r'$\gamma$= '+str('%.2f' % eta))
ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=r'$K= %.0f \pm %.0f$' %(K.value, K.stderr))
ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=r'$r= %.3f \pm %.3f$' %(r.value, r.stderr))
ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=r'$\alpha= %.3f \pm %.3f$' %(alpha.value, alpha.stderr))
ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=first_death[lang][Type]+str(date_first1))
ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=first_plateau[lang]+str(int(K.value)))
#ax.plot(t[0],func(t[0],r,alpha,K), 'wo', alpha=0, label=r'$T_c= %.1f$' %(tc))
ax.plot(t[0],kind1[0], 'ko', alpha=0, label=estagio[lang])
#global stage
#print(stage)
if stage == 'early intermediate':
ax.plot(np.ones(100)*tjm,np.linspace(0,func(tjm,r,alpha,K),100), '-', alpha=1, color='darkorange')
ax.plot(t[0],kind1[0], 'o', alpha=0, label=fase[lang][1], color='darkorange')
if stage == 'late intermediate':
ax.plot(np.ones(100)*tjm,np.linspace(0,func(tjm,r,alpha,K),100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,func(tc,r,alpha,K),100), '-', alpha=1, color='yellow')
ax.plot(t[0],kind1[0], 'o', alpha=0, label=fase[lang][2], color='yellow')
if stage == 'early saturation':
ax.plot(np.ones(100)*tjm,np.linspace(0,func(tjm,r,alpha,K),100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,func(tc,r,alpha,K),100), '-', alpha=1, color='yellow')
ax.plot(np.ones(100)*tjp,np.linspace(0,func(tjp,r,alpha,K),100), '-', alpha=1, color='green')
ax.plot(t[0],kind1[0], 'o', alpha=0, label=fase[lang][3], color='green')
if stage == 'late saturation':
ax.plot(np.ones(100)*tjm,np.linspace(0,func(tjm,r,alpha,K),100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,func(tc,r,alpha,K),100), '-', alpha=1, color='yellow')
ax.plot(np.ones(100)*tjp,np.linspace(0,func(tjp,r,alpha,K),100), '-', alpha=1, color='green')
ax.plot(np.ones(100)*ts,np.linspace(0,func(ts,r,alpha,K),100), '-', alpha=1, color='royalblue')
ax.plot(t[0],kind1[0], 'o', alpha=0, label=fase[lang][4], color='royalblue')
ax.plot(t[0],kind1[0], 'o', alpha=0, label=r'$\copyright$2021 ModInterv')
ax.set_xlabel(label_x[lang][Type])
ax.set_ylabel(label_y[lang][Type])
ax.set_title(name+' '+methods_name[lang][1],y=1.05)
plt.figtext(0.51,0.9,time_dic[lang]+' UTC-3', fontsize=8,ha='center')
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True)
plt.xscale(dic[x_axis])
plt.yscale(dic[y_axis])
if x_axis=='Linear' and y_axis=='Linear':
legend = ax.legend(loc='best', fontsize='small')
line = legend.get_lines()
text = legend.get_texts()
text[len(text)-2].set_color(line[len(text)-2].get_color())
if x_axis=='Logarithmic':
ax.set_xlim(left=0.8)
if y_axis=='Logarithmic':
ax.set_ylim(bottom=0.8)
for spine in ('top', 'right', 'bottom', 'left'):
ax.spines[spine].set_visible(False)
plt.savefig('fig.png', format='png', dpi=200)
ax.set_rasterized(True)
plt.savefig('fig.eps', format='eps')
plt.show()
ratio = tc-max(t0)
abs_ratio = abs(tc-max(t0))
text_string = 'For the curve of '+str(Type)+' of '+str(name)+', the Richards model predicts an inflection point at Tc = '+str('%.1f' % tc)+' days after the first '+str(name2)+'. The total number of '+str(Type)+' predicted at the end of the epidemic is K = '+str('%.0f' % K.value)+' '+str(Type)+', assuming that the current trend continues.'
if ratio>7:
text_string1 = 'Warning: Tc is '+str('%.0f' % ratio)+' days beyond the last available datapoint. This is an indication that the epidemic curve is in its early growth regime. The user is therefore advised to use instead the q-exponential model.'
if ratio<=7 and ratio>=0:
text_string1 = 'Warning: Tc appears '+str('%.0f' % ratio)+' days after the last available datapoint. This indicates that the inflection is expected to occur soon. However, because Tc is too close to the last point, the current estimate for the appearance of the inflection point should be viewed as tentative only.'
if ratio>=-7 and ratio<0:
text_string1 = 'Warning: Tc appears '+str('%.0f' % abs_ratio)+' days before the last available datapoint. This indicates that the inflection has just occurred. However, because Tc is too close to the last point, the current estimate for the appearance of the inflection point should be viewed as tentative only.'
text_string_br = 'Para a curva de '+str(tipo)+' para '+str(name)+', o modelo Richards prevê um ponto de inflexão em Tc = '+str('%.1f' % tc)+' dias após '+str(tipo2)+'. O número total de '+str(tipo)+' previsto no final da epidemia é de K = '+str('%.0f' % K.value)+' '+str(tipo)+', supondo que a tendência atual se mantenha.'
if ratio>7:
text_string_br1 = 'Advertência: O ponto de inflexão, Tc, dá-se '+str('%.0f' % ratio)+' dias após a data do último dado mostrado no gráfico. Essa é uma indicação de que a curva epidêmica encontra-se no regime de crescimento inicial. Aconselha-se, portanto, a utilizar para esse caso o modelo q-exponencial.'
if ratio<=7 and ratio>=0:
text_string_br1 = 'Advertência: O ponto de inflexão, Tc, acontece '+str('%.0f' % abs_ratio)+' dias depois da data do último dado mostrado no gráfico. Isto indica que a inflexão acabou de ocorrer. No entanto, como o Tc está muito próximo do último ponto, a presente estimativa para o ponto de inflexão deve ser vista com cautela.'
if ratio>=-7 and ratio<0:
text_string_br1 = 'Advertência: O ponto de inflexão, Tc, acontece '+str('%.0f' % abs_ratio)+' dias antes da data do último dado mostrado no gráfico. Isto indica que a inflexão está prevista para acontecer em breve. No entanto, como o Tc está muito próximo do último ponto, a presente estimativa para o ponto de inflexão deve ser vista com cautela.'
if togg2==True:
daily=[]
for i in range(0,len(kind1)-1):
daily.append(kind1[i+1]-kind1[i])
fig, ax = plt.subplots(figsize=(6,4))
ax.plot(t0[0:len(t0)-1],daily,'ro-', markersize=4, alpha=0.5, label=label_y_daily[lang][Type]+' '+str(name))
ax.plot(t,funcprime(t,r,alpha,K), '-k', label=methods_name1[lang][1])
#ax.plot(tc*np.ones(100),np.linspace(0,funcprime(tc,r,alpha,K),100),'-',color='yellow',alpha=1,linewidth=2,label=r'$t_c$')
ax.plot(tc,funcprime(tc,r,alpha,K),'ok', markersize=8,alpha=1)
ax.plot(tc,funcprime(tc,r,alpha,K),'ok',alpha=0,label=first_peak[lang]+str(date_first_peak1))
if x_axis=='Linear' and y_axis=='Linear':
legend = ax.legend(loc='best', fontsize='small')
if x_axis=='Logarithmic':
ax.set_xlim(left=0.8)
if y_axis=='Logarithmic':
ax.set_ylim(bottom=0.8)
ax.set_xlabel(label_x[lang][Type])
ax.set_ylabel(label_y_daily[lang][Type])
ax.set_title(name+' '+methods_name[lang][1],y=1.05)
plt.figtext(0.51,0.9,time_dic[lang]+' UTC-3', fontsize=8,ha='center')
plt.xscale(dic[x_axis])
plt.yscale(dic[y_axis])
plt.savefig('fig1.png', format='png', dpi=200)
ax.set_rasterized(True)
plt.savefig('fig1.eps', format='eps')
plt.show()
image_local = 5
pdf = FPDF('P','cm','A4')
pdf.add_page()
pdf.set_font('Arial', 'B', 14)
pdf.cell(19, 1, name+' '+Type+' of COVID-19',0,1,'C')
pdf.set_font('Arial', '', 12)
pdf.multi_cell(19, 1, text_string)
if ratio>=-7:
pdf.multi_cell(19, 1, text_string1)
image_local = 10
pdf.image('fig.png', 1, image_local, 10)
if togg2==True:
if ratio>=-7:
pdf.image('fig1.png', 1, 17, 10)
else:
pdf.image('fig1.png', 1, 12, 10)
pdf.add_page()
pdf.set_font('Arial', 'B', 14)
pdf.cell(19, 1, tipo+' do COVID-19 para '+name ,0,1,'C')
pdf.set_font('Arial', '', 12)
pdf.multi_cell(19, 1, text_string_br)
if ratio>=-7:
pdf.multi_cell(19, 1, text_string_br1)
image_local = 10
pdf.image('fig.png', 1, image_local, 10)
if togg2==True:
if ratio>=-7:
pdf.image('fig1.png', 1, 17, 10)
else:
pdf.image('fig1.png', 1, 12, 10)
pdf.output('report.pdf', 'F')
def GRMPlot(Country,Type,name2,time_fit):
name = Country
# New form of data acquisition. It pulls updated data from the github repository used by JHU
url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv'
df_recovered = pd.read_csv(url, error_bad_lines=False)
url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv'
df_deaths = pd.read_csv(url, error_bad_lines=False)
url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv'
df_confirmed = pd.read_csv(url, error_bad_lines=False)
url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/UID_ISO_FIPS_LookUp_Table.csv'
df_population = pd.read_csv(url, error_bad_lines=False)
####### Data #######
df_deaths = df_deaths.groupby(['Country/Region']).aggregate(np.sum).reset_index()
df_confirmed = df_confirmed.groupby(['Country/Region']).aggregate(np.sum).reset_index()
df_recovered = df_recovered.groupby(['Country/Region']).aggregate(np.sum).reset_index()
#print(df_population[(df_population['Combined_Key'] == name)]['Population'].to_numpy(dtype='float'))
deaths = df_deaths[(df_deaths['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float')[0][np.where(df_deaths[(df_deaths['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float') > 0)[1][0]:len(df_deaths.columns.values)]
recovered = df_recovered[(df_recovered['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float')[0][np.where(df_deaths[(df_deaths['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float') > 0)[1][0]:len(df_deaths.columns.values)]
confirmed = df_confirmed[(df_confirmed['Country/Region'] == name)][df_confirmed.columns.values[3:len(df_confirmed.columns.values)]].to_numpy(dtype='float')[0][np.where(df_confirmed[(df_confirmed['Country/Region'] == name)][df_confirmed.columns.values[3:len(df_confirmed.columns.values)]].to_numpy(dtype='float') > 0)[1][0]:len(df_confirmed.columns.values)]
if Type == 'Deaths':
kind1 = deaths
if Type == 'Cases':
kind1 = confirmed
x = np.arange(0,len(kind1),1)
data=kind1
#plt.plot(x[:],data[:])
def f(y,q,alpha):
return hyp2f1(1,(1-q)/alpha,1+(1-q)/alpha,y)
# The implicity solution for the Generalized Richards Model is given by:
def f_imp(C,C0,q,alpha,r,K):
return (C**(1-q)*f((C/K)**alpha,q,alpha) - (C0)**(1-q)*f((C0/K)**alpha,q,alpha))/(r*(1-q))
def fcn2min(params, x, data):
"""q_exp Model"""
A = params['A']
r = params['r']
q = params['q']
#tc = params['tc']
model = (A+(1-q)*r*x)**(1/(1-q))
return model - data
# Initial values for parameters
params = Parameters()
params.add('A', value=1.3, min=0.)
params.add('q', value=.5, min=0., max=1.0)
params.add('r', value=.5, min=0.)
R2 = [0]
r2 =0
for i in range(0,10):
for j in range(25,50):
try:
minner = Minimizer(fcn2min, params, fcn_args=(x[i:j],data[i:j]))
result = minner.minimize(method='nelder')
result = minner.minimize(method='leastsq', params = result.params)
r2 = 1 - result.residual.var() / np.var(data)
if r2==1.0:
r2=0
if r2>max(R2):
A1 = result.params.get('A').value
q1 = result.params.get('q').value
r1 = result.params.get('r').value
C01= (A1)**(1/(1- q1))
R2.append(r2)
except:
pass
try:
def fcn2min(params, x, data):
"""Generalized Richard's Model."""
K = params['K']
alpha = params['alpha']
q = params['q']
r = params['r']
C0= params['C0']
#C0=1.
def f_ric(t,y,*args):
return r*(y)**q*(1-(y/K)**alpha)
sol = solve_ivp(f_ric,[0,len(x)],[C0],method='LSODA',t_eval=[i for i in x],rtol=1e-3,atol=1e-6)
model = sol.y[0]
return model - data
params = Parameters()
params.add('K', value=max(data))
params.add('alpha', value=0.5, min=0, vary = True)
params.add('q', value=q1, min=0, max=1.,vary=True)
params.add('r', value=r1, min=0,vary=True)
params.add('C0', value=C01, min=0., vary=True)
minner = Minimizer(fcn2min, params, fcn_args=(x[:],data[:]))
start =timer()
# first solve with Nelder-Mead algorithm
result = minner.minimize(method= 'least_squares')
end = timer()
#print('time spent first minimization:', end - start)
#print('first min: R2=',1 - result.residual.var() / np.var(data))
#lmfit.report_fit(result.params, min_correl=0.5)
start =timer()
result = minner.minimize(method= 'leastsq', params=result.params)
#print('second min: R2=',1 - result.residual.var() / np.var(data))
end = timer()
#print('time spent second minimization:', end - start)
#lmfit.report_fit(result.params, min_correl=0.5)
alpha2=result.params.get('alpha').value
q2=result.params.get('q').value
r2=result.params.get('r').value
K2=result.params.get('K').value
C02=result.params.get('C0').value
args2 = (alpha2, q2, r2, K2)
def f_ric1(t,y,*args2):
return r2*(y)**q2*(1-(y/K2)**alpha2)
t_forecast = 1 + 0/100
x1 = np.arange(0,t_forecast*max(x),1)
sol1 = solve_ivp(f_ric1,[0,len(x1)],[C02],method='LSODA',t_eval=[i for i in x1])
time1=sol1.t
sol1=sol1.y[0]
f2 = InterpolatedUnivariateSpline(time1, sol1)
der = []
for i in range(len(time1)):
h = 1e-4
der.append( ( f2(time1[i]+h)-f2(time1[i]-h) )/(2*h) )
f3 = InterpolatedUnivariateSpline(time1, der)
der2 = []
for i in range(len(time1)):
h = 1e-4
der2.append( ( f3(time1[i]+h)-f3(time1[i]-h) )/(2*h) )
der_max = max(der)
t_pos = der.index(der_max)
time_max = time1[t_pos]
except ValueError:
pass
return Country, Type, name2, x, kind1, alpha2, q2, r2, K2, C02
def GRMPlotPlot(Country, Type, name2, x, kind1, alpha2, q2, r2, beta2, K2, C02, Forecast, interv_time, interv_param,togg1, togg2, x_axis, y_axis, toggdetails1):
with warnings.catch_warnings():
warnings.simplefilter("ignore")
name = Country
data = kind1
t0 = x
args2 = (alpha2, q2, r2,beta2, K2)
def f_ric1(t,y,*args2):
return r2*(y)**q2*((1-(y/K2)**alpha2))**beta2
def f(y,beta,q,alpha):
return hyp2f1(beta,(1-q)/alpha,1+(1-q)/alpha,y)
def f_imp(C,C0,beta,q,alpha,r,K):
return (C**(1-q)*f((C/K)**alpha,beta,q,alpha) - (C0)**(1-q)*f((C0/K)**alpha,beta,q,alpha))/(r*(1-q))
#find the K that maps to the time of the last data point
def timefit(C_curve):
return (f_imp(C_curve*K2,C02,beta2,q2,alpha2,r2,K2) - x[-1])**2
def timefit1(C_curve):
return (f_imp(C_curve*K2,C02,beta2,q2,alpha2,r2,K2) - (x[-1]+Forecast))**2
res = minimize_scalar(timefit, bounds=(0, 1), method='bounded')
C_plot = res.x
res1 = minimize_scalar(timefit1, bounds=(0, 1), method='bounded')
C_plot1 = res1.x
const = 0.9999 - C_plot
t_forecast1 = const*Forecast/100
c_forecast = np.arange(C02,C_plot*K2+1+Forecast,0.1)
alpha3 = alpha2
q3 = q2
beta3 = beta2
alpha2 = np.float(alpha2.value)
q2 = np.float(q2.value)
beta2 = np.float(beta2.value)
c_tc = ( ( q2 )**( -1 ) * ( alpha2 * beta2 + q2 ) )**( -1 * ( alpha2 )**( -1 ) )
C_tc = K2*c_tc
c_ts = ( 6 )**( -1 * ( alpha2 )**( -1 ) ) * ( ( ( 6 * ( alpha2 )**( 3 ) * ( beta2 )**( 3 ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) * ( -7 + 18 * q2 ) + ( q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) )**( -1 ) * ( 2 * ( ( alpha2 )**( 3 ) * beta2 * ( -1 + 7 * beta2 ) + ( 3 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( 4 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) + ( alpha2 )**( 2 ) * beta2 * ( -3 + ( 7 * q2 + beta2 * ( -7 + 18 * q2 ) ) ) ) ) ) + ( 2 * ( 2 )**( 1/3 ) * ( alpha2 )**( 2 ) * beta2 * ( -3 * ( 3 + ( 3 * alpha2 + -7 * q2 ) ) * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 3 ) * ( 13 + ( 31 * ( alpha2 )**( 2 ) + 2 * alpha2 * ( -22 + 63 * q2 ) ) ) + ( alpha2 * ( beta2 )**( 2 ) * ( -14 + ( -14 * ( alpha2 )**( 3 ) + ( 26 * q2 + ( ( alpha2 )**( 2 ) * ( -7 + 8 * q2 ) + alpha2 * ( 35 + ( -289 * q2 + 378 * ( q2 )**( 2 ) ) ) ) ) ) ) + beta2 * ( 4 + ( ( alpha2 )**( 4 ) + ( ( alpha2 )**( 3 ) * ( 6 + -14 * q2 ) + ( -14 * q2 + ( 13 * ( q2 )**( 2 ) + ( ( alpha2 )**( 2 ) * ( -5 + ( 56 * q2 + -77 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * ( -3 + ( 70 * q2 + ( -223 * ( q2 )**( 2 ) + 189 * ( q2 )**( 3 ) ) ) ) ) ) ) ) ) ) ) ) ) * ( ( 2 * ( alpha2 )**( 9 ) * ( beta2 )**( 3 ) * ( -1 + ( 21 * beta2 + ( -120 * ( beta2 )**( 2 ) + 154 * ( beta2 )**( 3 ) ) ) ) + ( 3 * ( alpha2 )**( 8 ) * ( beta2 )**( 3 ) * ( -6 + ( beta2 * ( 49 + -106 * q2 ) + ( 14 * q2 + ( -49 * ( beta2 )**( 2 ) * ( -1 + 4 * q2 ) + 2 * ( beta2 )**( 3 ) * ( -91 + 279 * q2 ) ) ) ) ) + ( 3 * ( alpha2 )**( 7 ) * ( beta2 )**( 3 ) * ( -4 + ( -14 * q2 + ( 28 * ( q2 )**( 2 ) + ( ( beta2 )**( 3 ) * ( 56 + 90 * q2 ) + ( 3 * ( beta2 )**( 2 ) * ( 32 + ( -399 * q2 + 528 * ( q2 )**( 2 ) ) ) + -2 * beta2 * ( 56 + ( -451 * q2 + 602 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( 70 * ( beta2 )**( 4 ) + ( 27 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( 3 * ( beta2 )**( 3 ) * ( -19 + ( -189 * q2 + 360 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 90 + ( -1596 * q2 + ( 4971 * ( q2 )**( 2 ) + -4228 * ( q2 )**( 3 ) ) ) ) + 3 * ( beta2 )**( 2 ) * ( 49 + ( 216 * q2 + ( -987 * ( q2 )**( 2 ) + 756 * ( q2 )**( 3 ) ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 5 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 3 ) * ( -46 + 70 * q2 ) + ( 12 * q2 * ( 6 + ( -35 * q2 + ( 67 * ( q2 )**( 2 ) + -42 * ( q2 )**( 3 ) ) ) ) + ( ( beta2 )**( 2 ) * ( -28 + ( 394 * q2 + ( -903 * ( q2 )**( 2 ) + 540 * ( q2 )**( 3 ) ) ) ) + beta2 * ( -26 + ( 238 * q2 + ( -1170 * ( q2 )**( 2 ) + ( 2401 * ( q2 )**( 3 ) + -1656 * ( q2 )**( 4 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 3 ) * beta2 * ( -54 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 2 * ( beta2 )**( 2 ) * ( -8 + ( 42 * q2 + ( -69 * ( q2 )**( 2 ) + 35 * ( q2 )**( 3 ) ) ) ) + 9 * beta2 * q2 * ( 12 + ( -70 * q2 + ( 144 * ( q2 )**( 2 ) + ( -119 * ( q2 )**( 3 ) + 30 * ( q2 )**( 4 ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 4 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 2 ) * ( 28 + ( -92 * q2 + 70 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 12 + ( -238 * q2 + ( 845 * ( q2 )**( 2 ) + ( -1015 * ( q2 )**( 3 ) + 360 * ( q2 )**( 4 ) ) ) ) ) + -3 * q2 * ( 18 + ( -231 * q2 + ( 868 * ( q2 )**( 2 ) + ( -1295 * ( q2 )**( 3 ) + 678 * ( q2 )**( 4 ) ) ) ) ) ) ) + 3 * ( 3 )**( 1/2 ) * ( -1 * ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( ( 6 * ( alpha2 )**( 3 ) * ( beta2 )**( 3 ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) * ( -7 + 18 * q2 ) + ( q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) )**( 2 ) * ( -108 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 4 * beta2 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) * ( 27 + ( 63 * q2 + ( -396 * ( q2 )**( 2 ) + ( 343 * ( q2 )**( 3 ) + ( -9 * ( alpha2 )**( 2 ) * ( -3 + 7 * q2 ) + -54 * alpha2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( beta2 )**( 4 ) * ( 4 + ( 25 * ( alpha2 )**( 6 ) + ( 4 * alpha2 * ( -29 + 56 * q2 ) + ( 2 * ( alpha2 )**( 5 ) * ( -88 + 217 * q2 ) + ( ( alpha2 )**( 2 ) * ( 425 + ( -2114 * q2 + 2473 * ( q2 )**( 2 ) ) ) + ( ( alpha2 )**( 4 ) * ( 482 + ( -2534 * q2 + 3193 * ( q2 )**( 2 ) ) ) + ( alpha2 )**( 3 ) * ( -644 + ( 3990 * q2 + ( -9554 * ( q2 )**( 2 ) + 8232 * ( q2 )**( 3 ) ) ) ) ) ) ) ) ) ) + ( -2 * ( beta2 )**( 3 ) * ( 7 * ( alpha2 )**( 6 ) + ( ( alpha2 )**( 5 ) * ( -42 + 103 * q2 ) + ( -4 * ( 7 + ( -29 * q2 + 28 * ( q2 )**( 2 ) ) ) + ( 7 * ( alpha2 )**( 4 ) * ( 12 + ( -61 * q2 + 77 * ( q2 )**( 2 ) ) ) + ( alpha2 * ( 84 + ( -1142 * q2 + ( 3171 * ( q2 )**( 2 ) + -2473 * ( q2 )**( 3 ) ) ) ) + ( -1 * ( alpha2 )**( 3 ) * ( 42 + ( 365 * q2 + ( -1533 * ( q2 )**( 2 ) + 1249 * ( q2 )**( 3 ) ) ) ) + ( alpha2 )**( 2 ) * ( -63 + ( 1715 * q2 + ( -9667 * ( q2 )**( 2 ) + ( 19108 * ( q2 )**( 3 ) + -12348 * ( q2 )**( 4 ) ) ) ) ) ) ) ) ) ) ) + ( beta2 )**( 2 ) * ( -28 + ( ( alpha2 )**( 6 ) + ( -336 * q2 + ( 2284 * ( q2 )**( 2 ) + ( -4228 * ( q2 )**( 3 ) + ( 2473 * ( q2 )**( 4 ) + ( -2 * ( alpha2 )**( 4 ) * ( 15 + ( -84 * q2 + 103 * ( q2 )**( 2 ) ) ) + ( -18 * ( alpha2 )**( 3 ) * ( -6 + ( 70 * q2 + ( -201 * ( q2 )**( 2 ) + 168 * ( q2 )**( 3 ) ) ) ) + ( ( alpha2 )**( 2 ) * ( -159 + ( 1680 * q2 + ( -7154 * ( q2 )**( 2 ) + ( 13076 * ( q2 )**( 3 ) + -8471 * ( q2 )**( 4 ) ) ) ) ) + 2 * alpha2 * ( 54 + ( -126 * q2 + ( -3213 * ( q2 )**( 2 ) + ( 15344 * ( q2 )**( 3 ) + ( -23885 * ( q2 )**( 4 ) + 12348 * ( q2 )**( 5 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) )**( 1/2 ) ) ) ) ) ) ) ) )**( -1/3 ) + ( 2 )**( 2/3 ) * ( ( 2 * ( alpha2 )**( 9 ) * ( beta2 )**( 3 ) * ( -1 + ( 21 * beta2 + ( -120 * ( beta2 )**( 2 ) + 154 * ( beta2 )**( 3 ) ) ) ) + ( 3 * ( alpha2 )**( 8 ) * ( beta2 )**( 3 ) * ( -6 + ( beta2 * ( 49 + -106 * q2 ) + ( 14 * q2 + ( -49 * ( beta2 )**( 2 ) * ( -1 + 4 * q2 ) + 2 * ( beta2 )**( 3 ) * ( -91 + 279 * q2 ) ) ) ) ) + ( 3 * ( alpha2 )**( 7 ) * ( beta2 )**( 3 ) * ( -4 + ( -14 * q2 + ( 28 * ( q2 )**( 2 ) + ( ( beta2 )**( 3 ) * ( 56 + 90 * q2 ) + ( 3 * ( beta2 )**( 2 ) * ( 32 + ( -399 * q2 + 528 * ( q2 )**( 2 ) ) ) + -2 * beta2 * ( 56 + ( -451 * q2 + 602 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( 70 * ( beta2 )**( 4 ) + ( 27 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( 3 * ( beta2 )**( 3 ) * ( -19 + ( -189 * q2 + 360 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 90 + ( -1596 * q2 + ( 4971 * ( q2 )**( 2 ) + -4228 * ( q2 )**( 3 ) ) ) ) + 3 * ( beta2 )**( 2 ) * ( 49 + ( 216 * q2 + ( -987 * ( q2 )**( 2 ) + 756 * ( q2 )**( 3 ) ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 5 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 3 ) * ( -46 + 70 * q2 ) + ( 12 * q2 * ( 6 + ( -35 * q2 + ( 67 * ( q2 )**( 2 ) + -42 * ( q2 )**( 3 ) ) ) ) + ( ( beta2 )**( 2 ) * ( -28 + ( 394 * q2 + ( -903 * ( q2 )**( 2 ) + 540 * ( q2 )**( 3 ) ) ) ) + beta2 * ( -26 + ( 238 * q2 + ( -1170 * ( q2 )**( 2 ) + ( 2401 * ( q2 )**( 3 ) + -1656 * ( q2 )**( 4 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 3 ) * beta2 * ( -54 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 2 * ( beta2 )**( 2 ) * ( -8 + ( 42 * q2 + ( -69 * ( q2 )**( 2 ) + 35 * ( q2 )**( 3 ) ) ) ) + 9 * beta2 * q2 * ( 12 + ( -70 * q2 + ( 144 * ( q2 )**( 2 ) + ( -119 * ( q2 )**( 3 ) + 30 * ( q2 )**( 4 ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 4 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 2 ) * ( 28 + ( -92 * q2 + 70 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 12 + ( -238 * q2 + ( 845 * ( q2 )**( 2 ) + ( -1015 * ( q2 )**( 3 ) + 360 * ( q2 )**( 4 ) ) ) ) ) + -3 * q2 * ( 18 + ( -231 * q2 + ( 868 * ( q2 )**( 2 ) + ( -1295 * ( q2 )**( 3 ) + 678 * ( q2 )**( 4 ) ) ) ) ) ) ) + 3 * ( 3 )**( 1/2 ) * ( -1 * ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( ( 6 * ( alpha2 )**( 3 ) * ( beta2 )**( 3 ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) * ( -7 + 18 * q2 ) + ( q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) )**( 2 ) * ( -108 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 4 * beta2 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) * ( 27 + ( 63 * q2 + ( -396 * ( q2 )**( 2 ) + ( 343 * ( q2 )**( 3 ) + ( -9 * ( alpha2 )**( 2 ) * ( -3 + 7 * q2 ) + -54 * alpha2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( beta2 )**( 4 ) * ( 4 + ( 25 * ( alpha2 )**( 6 ) + ( 4 * alpha2 * ( -29 + 56 * q2 ) + ( 2 * ( alpha2 )**( 5 ) * ( -88 + 217 * q2 ) + ( ( alpha2 )**( 2 ) * ( 425 + ( -2114 * q2 + 2473 * ( q2 )**( 2 ) ) ) + ( ( alpha2 )**( 4 ) * ( 482 + ( -2534 * q2 + 3193 * ( q2 )**( 2 ) ) ) + ( alpha2 )**( 3 ) * ( -644 + ( 3990 * q2 + ( -9554 * ( q2 )**( 2 ) + 8232 * ( q2 )**( 3 ) ) ) ) ) ) ) ) ) ) + ( -2 * ( beta2 )**( 3 ) * ( 7 * ( alpha2 )**( 6 ) + ( ( alpha2 )**( 5 ) * ( -42 + 103 * q2 ) + ( -4 * ( 7 + ( -29 * q2 + 28 * ( q2 )**( 2 ) ) ) + ( 7 * ( alpha2 )**( 4 ) * ( 12 + ( -61 * q2 + 77 * ( q2 )**( 2 ) ) ) + ( alpha2 * ( 84 + ( -1142 * q2 + ( 3171 * ( q2 )**( 2 ) + -2473 * ( q2 )**( 3 ) ) ) ) + ( -1 * ( alpha2 )**( 3 ) * ( 42 + ( 365 * q2 + ( -1533 * ( q2 )**( 2 ) + 1249 * ( q2 )**( 3 ) ) ) ) + ( alpha2 )**( 2 ) * ( -63 + ( 1715 * q2 + ( -9667 * ( q2 )**( 2 ) + ( 19108 * ( q2 )**( 3 ) + -12348 * ( q2 )**( 4 ) ) ) ) ) ) ) ) ) ) ) + ( beta2 )**( 2 ) * ( -28 + ( ( alpha2 )**( 6 ) + ( -336 * q2 + ( 2284 * ( q2 )**( 2 ) + ( -4228 * ( q2 )**( 3 ) + ( 2473 * ( q2 )**( 4 ) + ( -2 * ( alpha2 )**( 4 ) * ( 15 + ( -84 * q2 + 103 * ( q2 )**( 2 ) ) ) + ( -18 * ( alpha2 )**( 3 ) * ( -6 + ( 70 * q2 + ( -201 * ( q2 )**( 2 ) + 168 * ( q2 )**( 3 ) ) ) ) + ( ( alpha2 )**( 2 ) * ( -159 + ( 1680 * q2 + ( -7154 * ( q2 )**( 2 ) + ( 13076 * ( q2 )**( 3 ) + -8471 * ( q2 )**( 4 ) ) ) ) ) + 2 * alpha2 * ( 54 + ( -126 * q2 + ( -3213 * ( q2 )**( 2 ) + ( 15344 * ( q2 )**( 3 ) + ( -23885 * ( q2 )**( 4 ) + 12348 * ( q2 )**( 5 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) )**( 1/2 ) ) ) ) ) ) ) ) )**( 1/3 ) ) ) )**( ( alpha2 )**( -1 ) )
C_ts = K2*np.real(c_ts)
c_tjm = ( ( ( 4 * ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) + ( 2 * q2 * ( -1 + 2 * q2 ) + 2 * alpha2 * beta2 * ( -1 + 4 * q2 ) ) ) )**( -1 ) * ( ( alpha2 )**( 2 ) * beta2 + ( 2 * q2 * ( -1 + 2 * q2 ) + -1 * alpha2 * ( beta2 + ( -4 * beta2 * q2 + ( beta2 )**( 1/2 ) * ( ( 4 * q2 * ( -1 + 2 * q2 ) + beta2 * ( 1 + ( -2 * alpha2 + ( ( alpha2 )**( 2 ) + 8 * alpha2 * q2 ) ) ) ) )**( 1/2 ) ) ) ) ) )**( ( alpha2 )**( -1 ) )
C_tjm = K2*c_tjm
c_tjp = ( ( ( 4 * ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) + ( 2 * q2 * ( -1 + 2 * q2 ) + 2 * alpha2 * beta2 * ( -1 + 4 * q2 ) ) ) )**( -1 ) * ( ( alpha2 )**( 2 ) * beta2 + ( 2 * q2 * ( -1 + 2 * q2 ) + alpha2 * ( -1 * beta2 + ( 4 * beta2 * q2 + ( beta2 )**( 1/2 ) * ( ( 4 * q2 * ( -1 + 2 * q2 ) + beta2 * ( 1 + ( -2 * alpha2 + ( ( alpha2 )**( 2 ) + 8 * alpha2 * q2 ) ) ) ) )**( 1/2 ) ) ) ) ) )**( ( alpha2 )**( -1 ) )
C_tjp = K2*c_tjp
alpha2 = alpha3
q2 = q3
beta2 = beta3
tc = f_imp(C_tc,C02,beta2,q2,alpha2,r2,K2)
tjm = f_imp(C_tjm,C02,beta2,q2,alpha2,r2,K2)
tjp = f_imp(C_tjp,C02,beta2,q2,alpha2,r2,K2)
ts = f_imp(C_ts,C02,beta2,q2,alpha2,r2,K2)
global date_first
if Type == 'Deaths':
date_first_peak1 = date_first + timedelta(days=int(tc))
if Type == 'Cases':
date_first_peak1 = date_first + timedelta(days=int(tc))
lang = html_widget_lang.value
if lang == 'Por':
date_first1 = date_first.strftime('%d/%m/%y')
date_first_peak1 = date_first_peak1.strftime('%d/%m/%y')
if lang == 'Eng':
date_first1 = date_first.strftime('%m/%d/%y')
date_first_peak1 = date_first_peak1.strftime('%m/%d/%y')
if Type=='Deaths':
tipo = 'Mortes'
tipo1 = 'morte'
tipo2 = 'a primeira morte'
if Type=='Cases':
tipo = 'Casos'
tipo1 = 'caso'
tipo2 = 'o primeiro caso'
fig, ax = plt.subplots(figsize=(6,4))
ax.plot(x,data, 'ro',alpha=0.5, label = str(subs[lang][Type])+' '+str(name))
ax.plot(f_imp(c_forecast,C02,beta2,q2,alpha2,r2,K2),c_forecast, '-k', label = methods_name1[lang][2])
if togg1 == True:
if interv_param>=0:
ax.plot(time3[0:len(time3)-1],sol3[0:len(sol3)-1], '-b', label = 'Intervention')
if interv_param<0:
ax.plot(time3[0:len(time3)-1],sol3[0:len(sol3)-1], '-r', label = 'Relaxation')
ax.plot(x[0],data[0], 'wo', alpha=0, label=r'$\gamma$ = %.2f' %(eta))
ax.plot(x[0],data[0], 'wo', alpha=0, label=first_death[lang][Type]+str(date_first1))
#ax.plot(x[0],data[0], 'wo', alpha=0, label=first_plateau[lang]+str(int(K2.value)))
if Forecast>0:
ax.plot(x[0],data[0], 'wo', alpha=0, label=forecast_message[lang][str(Forecast)]+'+'+str(int(C_plot1*K2-C_plot*K2))+' '+subs[lang][Type])
ax.plot(x[0],data[0], 'wo', alpha=0, label=forecast_message[lang][str(Forecast)]+str(int(C_plot1*K2))+' '+subs[lang][Type])
#ax.plot(x[0],data[0], 'wo', alpha=0, label=r'$T_c$ = %.1f' %(tc))
ax.plot(x[0],data[0], 'ko', alpha=0, label=estagio[lang])
global stage
if stage == 'early intermediate':
ax.plot(np.ones(100)*tjm,np.linspace(0,C_tjm,100), '-', alpha=1, color='darkorange')
ax.plot(x[0],data[0], 'o', alpha=0, label=fase[lang][1], color='darkorange')
if stage == 'late intermediate':
ax.plot(np.ones(100)*tjm,np.linspace(0,C_tjm,100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,C_tc,100), '-', alpha=1, color='yellow')
ax.plot(x[0],data[0], 'o', alpha=0, label=fase[lang][2], color='yellow')
if stage == 'early saturation':
ax.plot(np.ones(100)*tjm,np.linspace(0,C_tjm,100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,C_tc,100), '-', alpha=1, color='yellow')
ax.plot(np.ones(100)*tjp,np.linspace(0,C_tjp,100), '-', alpha=1, color='green')
ax.plot(x[0],data[0], 'o', alpha=0, label=fase[lang][3], color='green')
if stage == 'late saturation':
ax.plot(np.ones(100)*tjm,np.linspace(0,C_tjm,100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,C_tc,100), '-', alpha=1, color='yellow')
ax.plot(np.ones(100)*tjp,np.linspace(0,C_tjp,100), '-', alpha=1, color='green')
ax.plot(np.ones(100)*ts,np.linspace(0,C_ts,100), '-', alpha=1, color='royalblue')
ax.plot(x[0],data[0], 'o', alpha=0, label=fase[lang][4], color='royalblue')
ax.plot(x[0],data[0], 'o', alpha=0, label=r'$\copyright$2020 ModInterv')
ax.set_title(name+' '+methods_name[lang][2],y=1.05)
plt.figtext(0.51,0.9,time_dic[lang]+' UTC-3', fontsize=8,ha='center')
ax.set_xlabel(label_x[lang][Type])
ax.set_ylabel(label_y[lang][Type])
plt.xscale(dic[x_axis])
plt.yscale(dic[y_axis])
if x_axis=='Linear' and y_axis=='Linear':
legend = ax.legend(loc='best', fontsize='small')
line = legend.get_lines()
text = legend.get_texts()
text[len(text)-2].set_color(line[len(text)-2].get_color())
if x_axis=='Logarithmic':
ax.set_xlim(left=0.8)
if y_axis=='Logarithmic':
ax.set_ylim(bottom=0.8)
plt.savefig('fig.png', format='png', dpi=200)
ax.set_rasterized(True)
plt.savefig('fig.eps', format='eps')
if toggdetails1 == True:
fig, ax = plt.subplots(figsize=(6,4))
ax.plot(x,data, 'ro',alpha=0.5, label = str(subs[lang][Type])+' '+str(name))
ax.plot(f_imp(c_forecast,C02,beta2,q2,alpha2,r2,K2),c_forecast, '-k', label = methods_name1[lang][2])
if togg1 == True:
if interv_param>=0:
ax.plot(time3[0:len(time3)-1],sol3[0:len(sol3)-1], '-b', label = 'Intervention')
if interv_param<0:
ax.plot(time3[0:len(time3)-1],sol3[0:len(sol3)-1], '-r', label = 'Relaxation')
ax.plot(x[0],data[0], 'wo', alpha=0, label=r'$\gamma$ = %.2f' %(eta))
ax.plot(x[0],data[0], 'wo', alpha=0, label=r'$K = %.0f \pm %.0f$' %(K2.value, K2.stderr))
ax.plot(x[0],data[0], 'wo', alpha=0, label=r'$\alpha = %.3f \pm %.3f$' %(alpha2.value, alpha2.stderr))
#ax.plot(x[0],data[0], 'wo', alpha=0, label=r'$p = %.3f \pm %.3f$' %(beta2.value, beta2.stderr))
ax.plot(x[0],data[0], 'wo', alpha=0, label=r'$q = %.3f \pm %.3f$' %(q2.value, q2.stderr))
ax.plot(x[0],data[0], 'wo', alpha=0, label=r'$r = %.3f \pm %.3f$' %(r2.value, r2.stderr))
ax.plot(x[0],data[0], 'wo', alpha=0, label=first_death[lang][Type]+str(date_first1))
ax.plot(x[0],data[0], 'wo', alpha=0, label=first_plateau[lang]+str(int(K2.value)))
#ax.plot(x[0],data[0], 'wo', alpha=0, label=r'$T_c$ = %.1f' %(tc))
ax.plot(x[0],data[0], 'ko', alpha=0, label=estagio[lang])
#global stage
if stage == 'early intermediate':
ax.plot(np.ones(100)*tjm,np.linspace(0,C_tjm,100), '-', alpha=1, color='darkorange')
ax.plot(x[0],data[0], 'o', alpha=0, label=fase[lang][1], color='darkorange')
if stage == 'late intermediate':
ax.plot(np.ones(100)*tjm,np.linspace(0,C_tjm,100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,C_tc,100), '-', alpha=1, color='yellow')
ax.plot(x[0],data[0], 'o', alpha=0, label=fase[lang][2], color='yellow')
if stage == 'early saturation':
ax.plot(np.ones(100)*tjm,np.linspace(0,C_tjm,100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,C_tc,100), '-', alpha=1, color='yellow')
ax.plot(np.ones(100)*tjp,np.linspace(0,C_tjp,100), '-', alpha=1, color='green')
ax.plot(x[0],data[0], 'o', alpha=0, label=fase[lang][3], color='green')
if stage == 'late saturation':
ax.plot(np.ones(100)*tjm,np.linspace(0,C_tjm,100), '-', alpha=1, color='darkorange')
ax.plot(np.ones(100)*tc,np.linspace(0,C_tc,100), '-', alpha=1, color='yellow')
ax.plot(np.ones(100)*tjp,np.linspace(0,C_tjp,100), '-', alpha=1, color='green')
ax.plot(np.ones(100)*ts,np.linspace(0,C_ts,100), '-', alpha=1, color='royalblue')
ax.plot(x[0],data[0], 'o', alpha=0, label=fase[lang][4], color='royalblue')
ax.plot(x[0],data[0], 'o', alpha=0, label=r'$\copyright$2020 ModInterv')
ax.set_title(name+' '+methods_name[lang][2],y=1.05)
plt.figtext(0.51,0.9,time_dic[lang]+' UTC-3', fontsize=8,ha='center')
ax.set_xlabel(label_x[lang][Type])
ax.set_ylabel(label_y[lang][Type])
plt.xscale(dic[x_axis])
plt.yscale(dic[y_axis])
if x_axis=='Linear' and y_axis=='Linear':
legend = ax.legend(loc='best', fontsize='small')
line = legend.get_lines()
text = legend.get_texts()
text[len(text)-2].set_color(line[len(text)-2].get_color())
if x_axis=='Logarithmic':
ax.set_xlim(left=0.8)
if y_axis=='Logarithmic':
ax.set_ylim(bottom=0.8)
plt.savefig('fig.png', format='png', dpi=200)
ax.set_rasterized(True)
plt.savefig('fig.eps', format='eps')
ratio = tc-max(t0)
abs_ratio = abs(tc-max(t0))
text_string = 'For the curve of '+str(Type)+' of '+str(name)+', the Beta Logistic model predicts an inflection point at Tc = '+str('%.1f' % tc)+' days after the first '+str(name2)+'. The total number of '+str(Type)+' predicted at the end of the epidemic is K = '+str('%.0f' % K2)+' '+str(Type)+', assuming that the current trend continues.'
if ratio>7:
text_string1 = 'Warning: Tc is '+str('%.0f' % ratio)+' days beyond the last available datapoint. This is an indication that the epidemic curve is in its early growth regime. The user is therefore advised to use instead the q-exponential model.'
if ratio<=7 and ratio>=0:
text_string1 = 'Warning: Tc appears '+str('%.0f' % ratio)+' days after the last available datapoint. This indicates that the inflection is expected to occur soon. However, because Tc is too close to the last point, the current estimate for the appearance of the inflection point should be viewed as tentative only.'
if ratio>=-7 and ratio<0:
text_string1 = 'Warning: Tc appears '+str('%.0f' % abs_ratio)+' days before the last available datapoint. This indicates that the inflection has just occurred. However, because Tc is too close to the last point, the current estimate for the appearance of the inflection point should be viewed as tentative only.'
text_string_br = 'Para a curva de '+str(tipo)+' para '+str(name)+', o modelo Beta LogÃstico prevê um ponto de inflexão em Tc = '+str('%.1f' % tc)+' dias após '+str(tipo2)+'. O número total de '+str(tipo)+' previsto no final da epidemia é de K = '+str('%.0f' % K2)+' '+str(tipo)+', supondo que a tendência atual se mantenha.'
if ratio>7:
text_string_br1 = 'Advertência: O ponto de inflexão, Tc, dá-se '+str('%.0f' % ratio)+' dias após a data do último dado mostrado no gráfico. Essa é uma indicação de que a curva epidêmica encontra-se no regime de crescimento inicial. Aconselha-se, portanto, a utilizar para esse caso o modelo q-exponencial.'
if ratio<=7 and ratio>=0:
text_string_br1 = 'Advertência: O ponto de inflexão, Tc, acontece '+str('%.0f' % abs_ratio)+' dias depois da data do último dado mostrado no gráfico. Isto indica que a inflexão acabou de ocorrer. No entanto, como o Tc está muito próximo do último ponto, a presente estimativa para o ponto de inflexão deve ser vista com cautela.'
if ratio>=-7 and ratio<0:
text_string_br1 = 'Advertência: O ponto de inflexão, Tc, acontece '+str('%.0f' % abs_ratio)+' dias antes da data do último dado mostrado no gráfico. Isto indica que a inflexão está prevista para acontecer em breve. No entanto, como o Tc está muito próximo do último ponto, a presente estimativa para o ponto de inflexão deve ser vista com cautela.'
if togg2==True:
daily=[]
daily.append(data[0])
for i in range(0,len(data)-1):
daily.append(data[i+1]-data[i])
fig, ax = plt.subplots(figsize=(6,4))
ax.plot(x[0:len(x)],daily,'ro-', markersize=4, alpha=0.5, label=label_y_daily[lang][Type]+' '+str(name))
ax.plot(f_imp(c_forecast,C02,beta2,q2,alpha2,r2,K2),f_ric1(0,c_forecast,*args2), '-k', label=methods_name1[lang][2])
#ax.plot(tc*np.ones(100),np.linspace(0,f_ric1(0,C_tc,*args2),100),'-',color='yellow',alpha=1,linewidth=2,label=r'$t_c$')
ax.plot(tc,f_ric1(0,C_tc,*args2),'ok', markersize=8,alpha=1)
ax.plot(tc,f_ric1(0,C_tc,*args2),'ok',alpha=0,label=first_peak[lang]+str(date_first_peak1))
if x_axis=='Linear' and y_axis=='Linear':
legend = ax.legend(loc='best', fontsize='small')
if x_axis=='Logarithmic':
ax.set_xlim(left=0.8)
if y_axis=='Logarithmic':
ax.set_ylim(bottom=0.8)
ax.set_xlabel(label_x[lang][Type])
ax.set_ylabel(label_y_daily[lang][Type])
ax.set_title(name+' '+methods_name[lang][2],y=1.05)
plt.figtext(0.51,0.9,time_dic[lang]+' UTC-3', fontsize=8,ha='center')
plt.xscale(dic[x_axis])
plt.yscale(dic[y_axis])
plt.savefig('fig1.png', format='png', dpi=200)
ax.set_rasterized(True)
plt.savefig('fig1.eps', format='eps')
plt.show()
image_local = 5
pdf = FPDF('P','cm','A4')
pdf.add_page()
pdf.set_font('Arial', 'B', 14)
pdf.cell(19, 1, name+' '+Type+' of COVID-19',0,1,'C')
pdf.set_font('Arial', '', 12)
pdf.multi_cell(19, 1, text_string)
if ratio>=-7:
pdf.multi_cell(19, 1, text_string1)
image_local = 10
pdf.image('fig.png', 1, image_local, 10)
if togg2==True:
if ratio>=-7:
pdf.image('fig1.png', 1, 17, 10)
else:
pdf.image('fig1.png', 1, 12, 10)
pdf.add_page()
pdf.set_font('Arial', 'B', 14)
pdf.cell(19, 1, tipo+' do COVID-19 para '+name ,0,1,'C')
pdf.set_font('Arial', '', 12)
pdf.multi_cell(19, 1, text_string_br)
if ratio>=-7:
pdf.multi_cell(19, 1, text_string_br1)
image_local = 10
pdf.image('fig.png', 1, image_local, 10)
if togg2==True:
if ratio>=-7:
pdf.image('fig1.png', 1, 17, 10)
else:
pdf.image('fig1.png', 1, 12, 10)
pdf.output('report.pdf', 'F')
plt.show()
def BLM_fit(Country,Type,name2,time_fit):
with warnings.catch_warnings():
warnings.simplefilter("ignore")
name = Country
# New form of data acquisition. It pulls updated data from the github repository used by JHU
url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_recovered_global.csv'
df_recovered = pandas.read_csv(url, error_bad_lines=False)
url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv'
df_deaths = pandas.read_csv(url, error_bad_lines=False)
url = 'https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv'
df_confirmed = pandas.read_csv(url, error_bad_lines=False)
####### Italy data #######
df_deaths = df_deaths.groupby(['Country/Region']).aggregate(np.sum).reset_index()
df_confirmed = df_confirmed.groupby(['Country/Region']).aggregate(np.sum).reset_index()
#print(df_deaths['Country/Region'].values.tolist())
#print(df_confirmed['Country/Region'].values.tolist())
deaths = df_deaths[(df_deaths['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float')[0][np.where(df_deaths[(df_deaths['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float') > 0)[1][0]:len(df_deaths.columns.values)]
recovered = df_recovered[(df_recovered['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float')[0][np.where(df_deaths[(df_deaths['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float') > 0)[1][0]:len(df_deaths.columns.values)]
confirmed = df_confirmed[(df_confirmed['Country/Region'] == name)][df_confirmed.columns.values[3:len(df_confirmed.columns.values)]].to_numpy(dtype='float')[0][np.where(df_confirmed[(df_confirmed['Country/Region'] == name)][df_confirmed.columns.values[3:len(df_confirmed.columns.values)]].to_numpy(dtype='float') > 0)[1][0]:len(df_confirmed.columns.values)]
#active = confirmed - recovered - deaths
global date_first
if Type == 'Deaths':
date_first = df_deaths.columns[3+np.where(df_deaths[(df_deaths['Country/Region'] == name)][df_deaths.columns.values[3:len(df_deaths.columns.values)]].to_numpy(dtype='float') > 0)[1][0]]
date_first = datetime.strptime(date_first, '%m/%d/%y')
if Type == 'Cases':
date_first = df_confirmed.columns[3+np.where(df_confirmed[(df_confirmed['Country/Region'] == name)][df_confirmed.columns.values[3:len(df_confirmed.columns.values)]].to_numpy(dtype='float') > 0)[1][0]]
date_first = datetime.strptime(date_first, '%m/%d/%y')
if Type == 'Deaths':
kind1 = deaths
if Type == 'Cases':
kind1 = confirmed
kind1 = kind1[0:time_fit]
x = np.arange(0,len(kind1),1)
daily = []
daily.append(kind1[0])
for j in range(0,len(kind1)-1):
delta = kind1[j+1] - kind1[j]
daily.append(delta)
def f(y,beta,q,alpha):
return hyp2f1(beta,(1-q)/alpha,1+(1-q)/alpha,y)
def f_imp(C,C0,beta,q,alpha,r,K):
return (C**(1-q)*f((C/K)**alpha,beta,q,alpha) - (C0)**(1-q)*f((C0/K)**alpha,beta,q,alpha))/(r*(1-q))
def fcn2min(params, x, data):
"""hyper-generalized logistic Model (HGLM)."""
K = params['K']
alpha = params['alpha']
q = params['q']
r = params['r']
#C0= params['C0']
C0=kind1[0]
beta=params['beta']
model = f_imp(x,C0,beta,q,alpha,r,K)
return model - data
params = Parameters()
params.add('K', value=1.1*max(kind1),min=max(kind1))
params.add('alpha', value=0.6, min=0, vary = True)
params.add('q', value=0.9, min=0.,max=.99999, vary=True)
params.add('r', value=0.5, min=0.,vary=True)
try:
params.add('beta', value=1, min=1.,vary=True)
minner = Minimizer(fcn2min, params,
fcn_args=(kind1, x))
result = minner.minimize(method= 'least_squares')
#result = minner.minimize(method= 'least_squares', params = result.params)
result = minner.minimize(method= 'leastsq', params = result.params)
R2 = 1 - result.residual.var() / np.var(kind1)
global R2_generalc
R2_generalc = R2
params = result.params
alpha2=params['alpha']
q2=params['q']
r2=params['r']
beta2=params['beta']
K2=params['K']
C02=data[0]
# Calculating the inflection point t_c:
alpha2=np.float(params['alpha'].value)
q2=np.float(params['q'].value)
beta2=np.float(params['beta'].value)
c_tc = ( ( q2 )**( -1 ) * ( alpha2 * beta2 + q2 ) )**( -1 * ( alpha2 )**( -1 ) )
C_tc = K2*c_tc
c_ts = np.real( ( 6 )**( -1 * ( alpha2 )**( -1 ) ) * ( ( ( 6 * ( alpha2 )**( 3 ) * ( beta2 )**( 3 ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) * ( -7 + 18 * q2 ) + ( q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) )**( -1 ) * ( 2 * ( ( alpha2 )**( 3 ) * beta2 * ( -1 + 7 * beta2 ) + ( 3 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( 4 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) + ( alpha2 )**( 2 ) * beta2 * ( -3 + ( 7 * q2 + beta2 * ( -7 + 18 * q2 ) ) ) ) ) ) + ( 2 * ( 2 )**( 1/3 ) * ( alpha2 )**( 2 ) * beta2 * ( -3 * ( 3 + ( 3 * alpha2 + -7 * q2 ) ) * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 3 ) * ( 13 + ( 31 * ( alpha2 )**( 2 ) + 2 * alpha2 * ( -22 + 63 * q2 ) ) ) + ( alpha2 * ( beta2 )**( 2 ) * ( -14 + ( -14 * ( alpha2 )**( 3 ) + ( 26 * q2 + ( ( alpha2 )**( 2 ) * ( -7 + 8 * q2 ) + alpha2 * ( 35 + ( -289 * q2 + 378 * ( q2 )**( 2 ) ) ) ) ) ) ) + beta2 * ( 4 + ( ( alpha2 )**( 4 ) + ( ( alpha2 )**( 3 ) * ( 6 + -14 * q2 ) + ( -14 * q2 + ( 13 * ( q2 )**( 2 ) + ( ( alpha2 )**( 2 ) * ( -5 + ( 56 * q2 + -77 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * ( -3 + ( 70 * q2 + ( -223 * ( q2 )**( 2 ) + 189 * ( q2 )**( 3 ) ) ) ) ) ) ) ) ) ) ) ) ) * ( ( 2 * ( alpha2 )**( 9 ) * ( beta2 )**( 3 ) * ( -1 + ( 21 * beta2 + ( -120 * ( beta2 )**( 2 ) + 154 * ( beta2 )**( 3 ) ) ) ) + ( 3 * ( alpha2 )**( 8 ) * ( beta2 )**( 3 ) * ( -6 + ( beta2 * ( 49 + -106 * q2 ) + ( 14 * q2 + ( -49 * ( beta2 )**( 2 ) * ( -1 + 4 * q2 ) + 2 * ( beta2 )**( 3 ) * ( -91 + 279 * q2 ) ) ) ) ) + ( 3 * ( alpha2 )**( 7 ) * ( beta2 )**( 3 ) * ( -4 + ( -14 * q2 + ( 28 * ( q2 )**( 2 ) + ( ( beta2 )**( 3 ) * ( 56 + 90 * q2 ) + ( 3 * ( beta2 )**( 2 ) * ( 32 + ( -399 * q2 + 528 * ( q2 )**( 2 ) ) ) + -2 * beta2 * ( 56 + ( -451 * q2 + 602 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( 70 * ( beta2 )**( 4 ) + ( 27 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( 3 * ( beta2 )**( 3 ) * ( -19 + ( -189 * q2 + 360 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 90 + ( -1596 * q2 + ( 4971 * ( q2 )**( 2 ) + -4228 * ( q2 )**( 3 ) ) ) ) + 3 * ( beta2 )**( 2 ) * ( 49 + ( 216 * q2 + ( -987 * ( q2 )**( 2 ) + 756 * ( q2 )**( 3 ) ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 5 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 3 ) * ( -46 + 70 * q2 ) + ( 12 * q2 * ( 6 + ( -35 * q2 + ( 67 * ( q2 )**( 2 ) + -42 * ( q2 )**( 3 ) ) ) ) + ( ( beta2 )**( 2 ) * ( -28 + ( 394 * q2 + ( -903 * ( q2 )**( 2 ) + 540 * ( q2 )**( 3 ) ) ) ) + beta2 * ( -26 + ( 238 * q2 + ( -1170 * ( q2 )**( 2 ) + ( 2401 * ( q2 )**( 3 ) + -1656 * ( q2 )**( 4 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 3 ) * beta2 * ( -54 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 2 * ( beta2 )**( 2 ) * ( -8 + ( 42 * q2 + ( -69 * ( q2 )**( 2 ) + 35 * ( q2 )**( 3 ) ) ) ) + 9 * beta2 * q2 * ( 12 + ( -70 * q2 + ( 144 * ( q2 )**( 2 ) + ( -119 * ( q2 )**( 3 ) + 30 * ( q2 )**( 4 ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 4 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 2 ) * ( 28 + ( -92 * q2 + 70 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 12 + ( -238 * q2 + ( 845 * ( q2 )**( 2 ) + ( -1015 * ( q2 )**( 3 ) + 360 * ( q2 )**( 4 ) ) ) ) ) + -3 * q2 * ( 18 + ( -231 * q2 + ( 868 * ( q2 )**( 2 ) + ( -1295 * ( q2 )**( 3 ) + 678 * ( q2 )**( 4 ) ) ) ) ) ) ) + 3 * ( 3 )**( 1/2 ) * ( -1 * ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( ( 6 * ( alpha2 )**( 3 ) * ( beta2 )**( 3 ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) * ( -7 + 18 * q2 ) + ( q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) )**( 2 ) * ( -108 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 4 * beta2 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) * ( 27 + ( 63 * q2 + ( -396 * ( q2 )**( 2 ) + ( 343 * ( q2 )**( 3 ) + ( -9 * ( alpha2 )**( 2 ) * ( -3 + 7 * q2 ) + -54 * alpha2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( beta2 )**( 4 ) * ( 4 + ( 25 * ( alpha2 )**( 6 ) + ( 4 * alpha2 * ( -29 + 56 * q2 ) + ( 2 * ( alpha2 )**( 5 ) * ( -88 + 217 * q2 ) + ( ( alpha2 )**( 2 ) * ( 425 + ( -2114 * q2 + 2473 * ( q2 )**( 2 ) ) ) + ( ( alpha2 )**( 4 ) * ( 482 + ( -2534 * q2 + 3193 * ( q2 )**( 2 ) ) ) + ( alpha2 )**( 3 ) * ( -644 + ( 3990 * q2 + ( -9554 * ( q2 )**( 2 ) + 8232 * ( q2 )**( 3 ) ) ) ) ) ) ) ) ) ) + ( -2 * ( beta2 )**( 3 ) * ( 7 * ( alpha2 )**( 6 ) + ( ( alpha2 )**( 5 ) * ( -42 + 103 * q2 ) + ( -4 * ( 7 + ( -29 * q2 + 28 * ( q2 )**( 2 ) ) ) + ( 7 * ( alpha2 )**( 4 ) * ( 12 + ( -61 * q2 + 77 * ( q2 )**( 2 ) ) ) + ( alpha2 * ( 84 + ( -1142 * q2 + ( 3171 * ( q2 )**( 2 ) + -2473 * ( q2 )**( 3 ) ) ) ) + ( -1 * ( alpha2 )**( 3 ) * ( 42 + ( 365 * q2 + ( -1533 * ( q2 )**( 2 ) + 1249 * ( q2 )**( 3 ) ) ) ) + ( alpha2 )**( 2 ) * ( -63 + ( 1715 * q2 + ( -9667 * ( q2 )**( 2 ) + ( 19108 * ( q2 )**( 3 ) + -12348 * ( q2 )**( 4 ) ) ) ) ) ) ) ) ) ) ) + ( beta2 )**( 2 ) * ( -28 + ( ( alpha2 )**( 6 ) + ( -336 * q2 + ( 2284 * ( q2 )**( 2 ) + ( -4228 * ( q2 )**( 3 ) + ( 2473 * ( q2 )**( 4 ) + ( -2 * ( alpha2 )**( 4 ) * ( 15 + ( -84 * q2 + 103 * ( q2 )**( 2 ) ) ) + ( -18 * ( alpha2 )**( 3 ) * ( -6 + ( 70 * q2 + ( -201 * ( q2 )**( 2 ) + 168 * ( q2 )**( 3 ) ) ) ) + ( ( alpha2 )**( 2 ) * ( -159 + ( 1680 * q2 + ( -7154 * ( q2 )**( 2 ) + ( 13076 * ( q2 )**( 3 ) + -8471 * ( q2 )**( 4 ) ) ) ) ) + 2 * alpha2 * ( 54 + ( -126 * q2 + ( -3213 * ( q2 )**( 2 ) + ( 15344 * ( q2 )**( 3 ) + ( -23885 * ( q2 )**( 4 ) + 12348 * ( q2 )**( 5 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) )**( 1/2 ) ) ) ) ) ) ) ) )**( -1/3 ) + ( 2 )**( 2/3 ) * ( ( 2 * ( alpha2 )**( 9 ) * ( beta2 )**( 3 ) * ( -1 + ( 21 * beta2 + ( -120 * ( beta2 )**( 2 ) + 154 * ( beta2 )**( 3 ) ) ) ) + ( 3 * ( alpha2 )**( 8 ) * ( beta2 )**( 3 ) * ( -6 + ( beta2 * ( 49 + -106 * q2 ) + ( 14 * q2 + ( -49 * ( beta2 )**( 2 ) * ( -1 + 4 * q2 ) + 2 * ( beta2 )**( 3 ) * ( -91 + 279 * q2 ) ) ) ) ) + ( 3 * ( alpha2 )**( 7 ) * ( beta2 )**( 3 ) * ( -4 + ( -14 * q2 + ( 28 * ( q2 )**( 2 ) + ( ( beta2 )**( 3 ) * ( 56 + 90 * q2 ) + ( 3 * ( beta2 )**( 2 ) * ( 32 + ( -399 * q2 + 528 * ( q2 )**( 2 ) ) ) + -2 * beta2 * ( 56 + ( -451 * q2 + 602 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( 70 * ( beta2 )**( 4 ) + ( 27 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( 3 * ( beta2 )**( 3 ) * ( -19 + ( -189 * q2 + 360 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 90 + ( -1596 * q2 + ( 4971 * ( q2 )**( 2 ) + -4228 * ( q2 )**( 3 ) ) ) ) + 3 * ( beta2 )**( 2 ) * ( 49 + ( 216 * q2 + ( -987 * ( q2 )**( 2 ) + 756 * ( q2 )**( 3 ) ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 5 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 3 ) * ( -46 + 70 * q2 ) + ( 12 * q2 * ( 6 + ( -35 * q2 + ( 67 * ( q2 )**( 2 ) + -42 * ( q2 )**( 3 ) ) ) ) + ( ( beta2 )**( 2 ) * ( -28 + ( 394 * q2 + ( -903 * ( q2 )**( 2 ) + 540 * ( q2 )**( 3 ) ) ) ) + beta2 * ( -26 + ( 238 * q2 + ( -1170 * ( q2 )**( 2 ) + ( 2401 * ( q2 )**( 3 ) + -1656 * ( q2 )**( 4 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 3 ) * beta2 * ( -54 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 2 * ( beta2 )**( 2 ) * ( -8 + ( 42 * q2 + ( -69 * ( q2 )**( 2 ) + 35 * ( q2 )**( 3 ) ) ) ) + 9 * beta2 * q2 * ( 12 + ( -70 * q2 + ( 144 * ( q2 )**( 2 ) + ( -119 * ( q2 )**( 3 ) + 30 * ( q2 )**( 4 ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 4 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 2 ) * ( 28 + ( -92 * q2 + 70 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 12 + ( -238 * q2 + ( 845 * ( q2 )**( 2 ) + ( -1015 * ( q2 )**( 3 ) + 360 * ( q2 )**( 4 ) ) ) ) ) + -3 * q2 * ( 18 + ( -231 * q2 + ( 868 * ( q2 )**( 2 ) + ( -1295 * ( q2 )**( 3 ) + 678 * ( q2 )**( 4 ) ) ) ) ) ) ) + 3 * ( 3 )**( 1/2 ) * ( -1 * ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( ( 6 * ( alpha2 )**( 3 ) * ( beta2 )**( 3 ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) * ( -7 + 18 * q2 ) + ( q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) )**( 2 ) * ( -108 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 4 * beta2 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) * ( 27 + ( 63 * q2 + ( -396 * ( q2 )**( 2 ) + ( 343 * ( q2 )**( 3 ) + ( -9 * ( alpha2 )**( 2 ) * ( -3 + 7 * q2 ) + -54 * alpha2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( beta2 )**( 4 ) * ( 4 + ( 25 * ( alpha2 )**( 6 ) + ( 4 * alpha2 * ( -29 + 56 * q2 ) + ( 2 * ( alpha2 )**( 5 ) * ( -88 + 217 * q2 ) + ( ( alpha2 )**( 2 ) * ( 425 + ( -2114 * q2 + 2473 * ( q2 )**( 2 ) ) ) + ( ( alpha2 )**( 4 ) * ( 482 + ( -2534 * q2 + 3193 * ( q2 )**( 2 ) ) ) + ( alpha2 )**( 3 ) * ( -644 + ( 3990 * q2 + ( -9554 * ( q2 )**( 2 ) + 8232 * ( q2 )**( 3 ) ) ) ) ) ) ) ) ) ) + ( -2 * ( beta2 )**( 3 ) * ( 7 * ( alpha2 )**( 6 ) + ( ( alpha2 )**( 5 ) * ( -42 + 103 * q2 ) + ( -4 * ( 7 + ( -29 * q2 + 28 * ( q2 )**( 2 ) ) ) + ( 7 * ( alpha2 )**( 4 ) * ( 12 + ( -61 * q2 + 77 * ( q2 )**( 2 ) ) ) + ( alpha2 * ( 84 + ( -1142 * q2 + ( 3171 * ( q2 )**( 2 ) + -2473 * ( q2 )**( 3 ) ) ) ) + ( -1 * ( alpha2 )**( 3 ) * ( 42 + ( 365 * q2 + ( -1533 * ( q2 )**( 2 ) + 1249 * ( q2 )**( 3 ) ) ) ) + ( alpha2 )**( 2 ) * ( -63 + ( 1715 * q2 + ( -9667 * ( q2 )**( 2 ) + ( 19108 * ( q2 )**( 3 ) + -12348 * ( q2 )**( 4 ) ) ) ) ) ) ) ) ) ) ) + ( beta2 )**( 2 ) * ( -28 + ( ( alpha2 )**( 6 ) + ( -336 * q2 + ( 2284 * ( q2 )**( 2 ) + ( -4228 * ( q2 )**( 3 ) + ( 2473 * ( q2 )**( 4 ) + ( -2 * ( alpha2 )**( 4 ) * ( 15 + ( -84 * q2 + 103 * ( q2 )**( 2 ) ) ) + ( -18 * ( alpha2 )**( 3 ) * ( -6 + ( 70 * q2 + ( -201 * ( q2 )**( 2 ) + 168 * ( q2 )**( 3 ) ) ) ) + ( ( alpha2 )**( 2 ) * ( -159 + ( 1680 * q2 + ( -7154 * ( q2 )**( 2 ) + ( 13076 * ( q2 )**( 3 ) + -8471 * ( q2 )**( 4 ) ) ) ) ) + 2 * alpha2 * ( 54 + ( -126 * q2 + ( -3213 * ( q2 )**( 2 ) + ( 15344 * ( q2 )**( 3 ) + ( -23885 * ( q2 )**( 4 ) + 12348 * ( q2 )**( 5 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) )**( 1/2 ) ) ) ) ) ) ) ) )**( 1/3 ) ) ) )**( ( alpha2 )**( -1 ) ) )
C_ts = K2*c_ts
c_tjm = ( ( ( 4 * ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) + ( 2 * q2 * ( -1 + 2 * q2 ) + 2 * alpha2 * beta2 * ( -1 + 4 * q2 ) ) ) )**( -1 ) * ( ( alpha2 )**( 2 ) * beta2 + ( 2 * q2 * ( -1 + 2 * q2 ) + -1 * alpha2 * ( beta2 + ( -4 * beta2 * q2 + ( beta2 )**( 1/2 ) * ( ( 4 * q2 * ( -1 + 2 * q2 ) + beta2 * ( 1 + ( -2 * alpha2 + ( ( alpha2 )**( 2 ) + 8 * alpha2 * q2 ) ) ) ) )**( 1/2 ) ) ) ) ) )**( ( alpha2 )**( -1 ) )
C_tjm = K2*c_tjm
c_tjp = ( ( ( 4 * ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) + ( 2 * q2 * ( -1 + 2 * q2 ) + 2 * alpha2 * beta2 * ( -1 + 4 * q2 ) ) ) )**( -1 ) * ( ( alpha2 )**( 2 ) * beta2 + ( 2 * q2 * ( -1 + 2 * q2 ) + alpha2 * ( -1 * beta2 + ( 4 * beta2 * q2 + ( beta2 )**( 1/2 ) * ( ( 4 * q2 * ( -1 + 2 * q2 ) + beta2 * ( 1 + ( -2 * alpha2 + ( ( alpha2 )**( 2 ) + 8 * alpha2 * q2 ) ) ) ) )**( 1/2 ) ) ) ) ) )**( ( alpha2 )**( -1 ) )
C_tjp = K2*c_tjp
alpha2=params['alpha']
q2=params['q']
beta2=params['beta']
tc = f_imp(C_tc,C02,beta2,q2,alpha2,r2,K2)
tjm = f_imp(C_tjm,C02,beta2,q2,alpha2,r2,K2)
tjp = f_imp(C_tjp,C02,beta2,q2,alpha2,r2,K2)
ts = f_imp(C_ts,C02,beta2,q2,alpha2,r2,K2)
global stage
if x[-1]>tjm and x[-1]<=tc:
stage = 'early intermediate'
if x[-1]>tc and x[-1]<=tjp:
stage = 'late intermediate'
if x[-1]>tjp and x[-1]<=ts:
stage = 'early saturation'
if x[-1]>ts:
stage = 'late saturation'
test_vec = result.covar
test_vec1 = [alpha2.stderr/alpha2.value, q2.stderr/q2.value, r2.stderr/r2.value, beta2.stderr/beta2.value, K2.stderr/K2.value]
for i in range(0, len(test_vec1)):
if test_vec1[i]>0.2:
raise ValueError
Method1 = 'BLM'
return Method1, Country, Type, name2, x, kind1, alpha2, q2, r2, beta2, K2, C02
except:
try:
params.add('beta', value=1, min=1.,vary=False)
minner = Minimizer(fcn2min, params,
fcn_args=(kind1, x))
result = minner.minimize(method= 'least_squares')
#result = minner.minimize(method= 'least_squares', params = result.params)
result = minner.minimize(method= 'leastsq', params = result.params)
R2 = 1 - result.residual.var() / np.var(kind1)
#global R2_generalc
R2_generalc = R2
params = result.params
alpha2=params['alpha']
q2=params['q']
r2=params['r']
beta2=params['beta']
K2=params['K']
C02=data[0]
# Calculating the inflection point t_c:
alpha2=np.float(params['alpha'].value)
q2=np.float(params['q'].value)
beta2=np.float(params['beta'].value)
c_tc = ( ( q2 )**( -1 ) * ( alpha2 * beta2 + q2 ) )**( -1 * ( alpha2 )**( -1 ) )
C_tc = K2*c_tc
c_ts = np.real( ( 6 )**( -1 * ( alpha2 )**( -1 ) ) * ( ( ( 6 * ( alpha2 )**( 3 ) * ( beta2 )**( 3 ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) * ( -7 + 18 * q2 ) + ( q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) )**( -1 ) * ( 2 * ( ( alpha2 )**( 3 ) * beta2 * ( -1 + 7 * beta2 ) + ( 3 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( 4 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) + ( alpha2 )**( 2 ) * beta2 * ( -3 + ( 7 * q2 + beta2 * ( -7 + 18 * q2 ) ) ) ) ) ) + ( 2 * ( 2 )**( 1/3 ) * ( alpha2 )**( 2 ) * beta2 * ( -3 * ( 3 + ( 3 * alpha2 + -7 * q2 ) ) * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 3 ) * ( 13 + ( 31 * ( alpha2 )**( 2 ) + 2 * alpha2 * ( -22 + 63 * q2 ) ) ) + ( alpha2 * ( beta2 )**( 2 ) * ( -14 + ( -14 * ( alpha2 )**( 3 ) + ( 26 * q2 + ( ( alpha2 )**( 2 ) * ( -7 + 8 * q2 ) + alpha2 * ( 35 + ( -289 * q2 + 378 * ( q2 )**( 2 ) ) ) ) ) ) ) + beta2 * ( 4 + ( ( alpha2 )**( 4 ) + ( ( alpha2 )**( 3 ) * ( 6 + -14 * q2 ) + ( -14 * q2 + ( 13 * ( q2 )**( 2 ) + ( ( alpha2 )**( 2 ) * ( -5 + ( 56 * q2 + -77 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * ( -3 + ( 70 * q2 + ( -223 * ( q2 )**( 2 ) + 189 * ( q2 )**( 3 ) ) ) ) ) ) ) ) ) ) ) ) ) * ( ( 2 * ( alpha2 )**( 9 ) * ( beta2 )**( 3 ) * ( -1 + ( 21 * beta2 + ( -120 * ( beta2 )**( 2 ) + 154 * ( beta2 )**( 3 ) ) ) ) + ( 3 * ( alpha2 )**( 8 ) * ( beta2 )**( 3 ) * ( -6 + ( beta2 * ( 49 + -106 * q2 ) + ( 14 * q2 + ( -49 * ( beta2 )**( 2 ) * ( -1 + 4 * q2 ) + 2 * ( beta2 )**( 3 ) * ( -91 + 279 * q2 ) ) ) ) ) + ( 3 * ( alpha2 )**( 7 ) * ( beta2 )**( 3 ) * ( -4 + ( -14 * q2 + ( 28 * ( q2 )**( 2 ) + ( ( beta2 )**( 3 ) * ( 56 + 90 * q2 ) + ( 3 * ( beta2 )**( 2 ) * ( 32 + ( -399 * q2 + 528 * ( q2 )**( 2 ) ) ) + -2 * beta2 * ( 56 + ( -451 * q2 + 602 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( 70 * ( beta2 )**( 4 ) + ( 27 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( 3 * ( beta2 )**( 3 ) * ( -19 + ( -189 * q2 + 360 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 90 + ( -1596 * q2 + ( 4971 * ( q2 )**( 2 ) + -4228 * ( q2 )**( 3 ) ) ) ) + 3 * ( beta2 )**( 2 ) * ( 49 + ( 216 * q2 + ( -987 * ( q2 )**( 2 ) + 756 * ( q2 )**( 3 ) ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 5 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 3 ) * ( -46 + 70 * q2 ) + ( 12 * q2 * ( 6 + ( -35 * q2 + ( 67 * ( q2 )**( 2 ) + -42 * ( q2 )**( 3 ) ) ) ) + ( ( beta2 )**( 2 ) * ( -28 + ( 394 * q2 + ( -903 * ( q2 )**( 2 ) + 540 * ( q2 )**( 3 ) ) ) ) + beta2 * ( -26 + ( 238 * q2 + ( -1170 * ( q2 )**( 2 ) + ( 2401 * ( q2 )**( 3 ) + -1656 * ( q2 )**( 4 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 3 ) * beta2 * ( -54 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 2 * ( beta2 )**( 2 ) * ( -8 + ( 42 * q2 + ( -69 * ( q2 )**( 2 ) + 35 * ( q2 )**( 3 ) ) ) ) + 9 * beta2 * q2 * ( 12 + ( -70 * q2 + ( 144 * ( q2 )**( 2 ) + ( -119 * ( q2 )**( 3 ) + 30 * ( q2 )**( 4 ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 4 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 2 ) * ( 28 + ( -92 * q2 + 70 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 12 + ( -238 * q2 + ( 845 * ( q2 )**( 2 ) + ( -1015 * ( q2 )**( 3 ) + 360 * ( q2 )**( 4 ) ) ) ) ) + -3 * q2 * ( 18 + ( -231 * q2 + ( 868 * ( q2 )**( 2 ) + ( -1295 * ( q2 )**( 3 ) + 678 * ( q2 )**( 4 ) ) ) ) ) ) ) + 3 * ( 3 )**( 1/2 ) * ( -1 * ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( ( 6 * ( alpha2 )**( 3 ) * ( beta2 )**( 3 ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) * ( -7 + 18 * q2 ) + ( q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) )**( 2 ) * ( -108 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 4 * beta2 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) * ( 27 + ( 63 * q2 + ( -396 * ( q2 )**( 2 ) + ( 343 * ( q2 )**( 3 ) + ( -9 * ( alpha2 )**( 2 ) * ( -3 + 7 * q2 ) + -54 * alpha2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( beta2 )**( 4 ) * ( 4 + ( 25 * ( alpha2 )**( 6 ) + ( 4 * alpha2 * ( -29 + 56 * q2 ) + ( 2 * ( alpha2 )**( 5 ) * ( -88 + 217 * q2 ) + ( ( alpha2 )**( 2 ) * ( 425 + ( -2114 * q2 + 2473 * ( q2 )**( 2 ) ) ) + ( ( alpha2 )**( 4 ) * ( 482 + ( -2534 * q2 + 3193 * ( q2 )**( 2 ) ) ) + ( alpha2 )**( 3 ) * ( -644 + ( 3990 * q2 + ( -9554 * ( q2 )**( 2 ) + 8232 * ( q2 )**( 3 ) ) ) ) ) ) ) ) ) ) + ( -2 * ( beta2 )**( 3 ) * ( 7 * ( alpha2 )**( 6 ) + ( ( alpha2 )**( 5 ) * ( -42 + 103 * q2 ) + ( -4 * ( 7 + ( -29 * q2 + 28 * ( q2 )**( 2 ) ) ) + ( 7 * ( alpha2 )**( 4 ) * ( 12 + ( -61 * q2 + 77 * ( q2 )**( 2 ) ) ) + ( alpha2 * ( 84 + ( -1142 * q2 + ( 3171 * ( q2 )**( 2 ) + -2473 * ( q2 )**( 3 ) ) ) ) + ( -1 * ( alpha2 )**( 3 ) * ( 42 + ( 365 * q2 + ( -1533 * ( q2 )**( 2 ) + 1249 * ( q2 )**( 3 ) ) ) ) + ( alpha2 )**( 2 ) * ( -63 + ( 1715 * q2 + ( -9667 * ( q2 )**( 2 ) + ( 19108 * ( q2 )**( 3 ) + -12348 * ( q2 )**( 4 ) ) ) ) ) ) ) ) ) ) ) + ( beta2 )**( 2 ) * ( -28 + ( ( alpha2 )**( 6 ) + ( -336 * q2 + ( 2284 * ( q2 )**( 2 ) + ( -4228 * ( q2 )**( 3 ) + ( 2473 * ( q2 )**( 4 ) + ( -2 * ( alpha2 )**( 4 ) * ( 15 + ( -84 * q2 + 103 * ( q2 )**( 2 ) ) ) + ( -18 * ( alpha2 )**( 3 ) * ( -6 + ( 70 * q2 + ( -201 * ( q2 )**( 2 ) + 168 * ( q2 )**( 3 ) ) ) ) + ( ( alpha2 )**( 2 ) * ( -159 + ( 1680 * q2 + ( -7154 * ( q2 )**( 2 ) + ( 13076 * ( q2 )**( 3 ) + -8471 * ( q2 )**( 4 ) ) ) ) ) + 2 * alpha2 * ( 54 + ( -126 * q2 + ( -3213 * ( q2 )**( 2 ) + ( 15344 * ( q2 )**( 3 ) + ( -23885 * ( q2 )**( 4 ) + 12348 * ( q2 )**( 5 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) )**( 1/2 ) ) ) ) ) ) ) ) )**( -1/3 ) + ( 2 )**( 2/3 ) * ( ( 2 * ( alpha2 )**( 9 ) * ( beta2 )**( 3 ) * ( -1 + ( 21 * beta2 + ( -120 * ( beta2 )**( 2 ) + 154 * ( beta2 )**( 3 ) ) ) ) + ( 3 * ( alpha2 )**( 8 ) * ( beta2 )**( 3 ) * ( -6 + ( beta2 * ( 49 + -106 * q2 ) + ( 14 * q2 + ( -49 * ( beta2 )**( 2 ) * ( -1 + 4 * q2 ) + 2 * ( beta2 )**( 3 ) * ( -91 + 279 * q2 ) ) ) ) ) + ( 3 * ( alpha2 )**( 7 ) * ( beta2 )**( 3 ) * ( -4 + ( -14 * q2 + ( 28 * ( q2 )**( 2 ) + ( ( beta2 )**( 3 ) * ( 56 + 90 * q2 ) + ( 3 * ( beta2 )**( 2 ) * ( 32 + ( -399 * q2 + 528 * ( q2 )**( 2 ) ) ) + -2 * beta2 * ( 56 + ( -451 * q2 + 602 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( 70 * ( beta2 )**( 4 ) + ( 27 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + ( 3 * ( beta2 )**( 3 ) * ( -19 + ( -189 * q2 + 360 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 90 + ( -1596 * q2 + ( 4971 * ( q2 )**( 2 ) + -4228 * ( q2 )**( 3 ) ) ) ) + 3 * ( beta2 )**( 2 ) * ( 49 + ( 216 * q2 + ( -987 * ( q2 )**( 2 ) + 756 * ( q2 )**( 3 ) ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 5 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 3 ) * ( -46 + 70 * q2 ) + ( 12 * q2 * ( 6 + ( -35 * q2 + ( 67 * ( q2 )**( 2 ) + -42 * ( q2 )**( 3 ) ) ) ) + ( ( beta2 )**( 2 ) * ( -28 + ( 394 * q2 + ( -903 * ( q2 )**( 2 ) + 540 * ( q2 )**( 3 ) ) ) ) + beta2 * ( -26 + ( 238 * q2 + ( -1170 * ( q2 )**( 2 ) + ( 2401 * ( q2 )**( 3 ) + -1656 * ( q2 )**( 4 ) ) ) ) ) ) ) ) + ( ( alpha2 )**( 3 ) * beta2 * ( -54 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 2 * ( beta2 )**( 2 ) * ( -8 + ( 42 * q2 + ( -69 * ( q2 )**( 2 ) + 35 * ( q2 )**( 3 ) ) ) ) + 9 * beta2 * q2 * ( 12 + ( -70 * q2 + ( 144 * ( q2 )**( 2 ) + ( -119 * ( q2 )**( 3 ) + 30 * ( q2 )**( 4 ) ) ) ) ) ) ) + ( 3 * ( alpha2 )**( 4 ) * ( beta2 )**( 2 ) * ( ( beta2 )**( 2 ) * ( 28 + ( -92 * q2 + 70 * ( q2 )**( 2 ) ) ) + ( beta2 * ( 12 + ( -238 * q2 + ( 845 * ( q2 )**( 2 ) + ( -1015 * ( q2 )**( 3 ) + 360 * ( q2 )**( 4 ) ) ) ) ) + -3 * q2 * ( 18 + ( -231 * q2 + ( 868 * ( q2 )**( 2 ) + ( -1295 * ( q2 )**( 3 ) + 678 * ( q2 )**( 4 ) ) ) ) ) ) ) + 3 * ( 3 )**( 1/2 ) * ( -1 * ( alpha2 )**( 6 ) * ( beta2 )**( 2 ) * ( ( 6 * ( alpha2 )**( 3 ) * ( beta2 )**( 3 ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) * ( -7 + 18 * q2 ) + ( q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) )**( 2 ) * ( -108 * ( q2 )**( 2 ) * ( ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) )**( 2 ) + ( 4 * beta2 * q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) * ( 27 + ( 63 * q2 + ( -396 * ( q2 )**( 2 ) + ( 343 * ( q2 )**( 3 ) + ( -9 * ( alpha2 )**( 2 ) * ( -3 + 7 * q2 ) + -54 * alpha2 * ( 1 + ( -7 * q2 + 9 * ( q2 )**( 2 ) ) ) ) ) ) ) ) + ( ( beta2 )**( 4 ) * ( 4 + ( 25 * ( alpha2 )**( 6 ) + ( 4 * alpha2 * ( -29 + 56 * q2 ) + ( 2 * ( alpha2 )**( 5 ) * ( -88 + 217 * q2 ) + ( ( alpha2 )**( 2 ) * ( 425 + ( -2114 * q2 + 2473 * ( q2 )**( 2 ) ) ) + ( ( alpha2 )**( 4 ) * ( 482 + ( -2534 * q2 + 3193 * ( q2 )**( 2 ) ) ) + ( alpha2 )**( 3 ) * ( -644 + ( 3990 * q2 + ( -9554 * ( q2 )**( 2 ) + 8232 * ( q2 )**( 3 ) ) ) ) ) ) ) ) ) ) + ( -2 * ( beta2 )**( 3 ) * ( 7 * ( alpha2 )**( 6 ) + ( ( alpha2 )**( 5 ) * ( -42 + 103 * q2 ) + ( -4 * ( 7 + ( -29 * q2 + 28 * ( q2 )**( 2 ) ) ) + ( 7 * ( alpha2 )**( 4 ) * ( 12 + ( -61 * q2 + 77 * ( q2 )**( 2 ) ) ) + ( alpha2 * ( 84 + ( -1142 * q2 + ( 3171 * ( q2 )**( 2 ) + -2473 * ( q2 )**( 3 ) ) ) ) + ( -1 * ( alpha2 )**( 3 ) * ( 42 + ( 365 * q2 + ( -1533 * ( q2 )**( 2 ) + 1249 * ( q2 )**( 3 ) ) ) ) + ( alpha2 )**( 2 ) * ( -63 + ( 1715 * q2 + ( -9667 * ( q2 )**( 2 ) + ( 19108 * ( q2 )**( 3 ) + -12348 * ( q2 )**( 4 ) ) ) ) ) ) ) ) ) ) ) + ( beta2 )**( 2 ) * ( -28 + ( ( alpha2 )**( 6 ) + ( -336 * q2 + ( 2284 * ( q2 )**( 2 ) + ( -4228 * ( q2 )**( 3 ) + ( 2473 * ( q2 )**( 4 ) + ( -2 * ( alpha2 )**( 4 ) * ( 15 + ( -84 * q2 + 103 * ( q2 )**( 2 ) ) ) + ( -18 * ( alpha2 )**( 3 ) * ( -6 + ( 70 * q2 + ( -201 * ( q2 )**( 2 ) + 168 * ( q2 )**( 3 ) ) ) ) + ( ( alpha2 )**( 2 ) * ( -159 + ( 1680 * q2 + ( -7154 * ( q2 )**( 2 ) + ( 13076 * ( q2 )**( 3 ) + -8471 * ( q2 )**( 4 ) ) ) ) ) + 2 * alpha2 * ( 54 + ( -126 * q2 + ( -3213 * ( q2 )**( 2 ) + ( 15344 * ( q2 )**( 3 ) + ( -23885 * ( q2 )**( 4 ) + 12348 * ( q2 )**( 5 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) )**( 1/2 ) ) ) ) ) ) ) ) )**( 1/3 ) ) ) )**( ( alpha2 )**( -1 ) ) )
C_ts = K2*c_ts
c_tjm = ( ( ( 4 * ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) + ( 2 * q2 * ( -1 + 2 * q2 ) + 2 * alpha2 * beta2 * ( -1 + 4 * q2 ) ) ) )**( -1 ) * ( ( alpha2 )**( 2 ) * beta2 + ( 2 * q2 * ( -1 + 2 * q2 ) + -1 * alpha2 * ( beta2 + ( -4 * beta2 * q2 + ( beta2 )**( 1/2 ) * ( ( 4 * q2 * ( -1 + 2 * q2 ) + beta2 * ( 1 + ( -2 * alpha2 + ( ( alpha2 )**( 2 ) + 8 * alpha2 * q2 ) ) ) ) )**( 1/2 ) ) ) ) ) )**( ( alpha2 )**( -1 ) )
C_tjm = K2*c_tjm
c_tjp = ( ( ( 4 * ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) + ( 2 * q2 * ( -1 + 2 * q2 ) + 2 * alpha2 * beta2 * ( -1 + 4 * q2 ) ) ) )**( -1 ) * ( ( alpha2 )**( 2 ) * beta2 + ( 2 * q2 * ( -1 + 2 * q2 ) + alpha2 * ( -1 * beta2 + ( 4 * beta2 * q2 + ( beta2 )**( 1/2 ) * ( ( 4 * q2 * ( -1 + 2 * q2 ) + beta2 * ( 1 + ( -2 * alpha2 + ( ( alpha2 )**( 2 ) + 8 * alpha2 * q2 ) ) ) ) )**( 1/2 ) ) ) ) ) )**( ( alpha2 )**( -1 ) )
C_tjp = K2*c_tjp
alpha2=params['alpha']
q2=params['q']
beta2=params['beta']
tc = f_imp(C_tc,C02,beta2,q2,alpha2,r2,K2)
tjm = f_imp(C_tjm,C02,beta2,q2,alpha2,r2,K2)
tjp = f_imp(C_tjp,C02,beta2,q2,alpha2,r2,K2)
ts = f_imp(C_ts,C02,beta2,q2,alpha2,r2,K2)
#print(ts,C02)
#global stage
if x[-1]>tjm and x[-1]<=tc:
stage = 'early intermediate'
if x[-1]>tc and x[-1]<=tjp:
stage = 'late intermediate'
if x[-1]>tjp and x[-1]<=ts:
stage = 'early saturation'
if x[-1]>ts:
stage = 'late saturation'
test_vec = result.covar
test_vec1 = [alpha2.stderr/alpha2.value, q2.stderr/q2.value, r2.stderr/r2.value, K2.stderr/K2.value]
for i in range(0, len(test_vec1)):
if test_vec1[i]>0.2:
raise ValueError
Method1 = 'GRM'
return Method1, Country, Type, name2, x, kind1, alpha2, q2, r2, beta2, K2, C02
except:
Method1 = 'RM'
Country, Type, name2, x, kind1, K, r, alpha, C0 = RMPlot(Country,Type,name2,time_fit)
return Method1, Country, Type, name2, x, kind1, alpha, 1, r, 1, K, C0
def BLM_plot(Country, Type, name2, x, kind1, alpha2, q2, r2, beta2, K2, C02, Forecast, interv_time, interv_param, togg1, togg2, x_axis, y_axis, toggdetails1):
with warnings.catch_warnings():
warnings.simplefilter("ignore")
name = Country
data = kind1
t0 = x
args2 = (alpha2, q2, r2,beta2, K2)
def f_ric1(t,y,*args2):
return r2*(y)**q2*((1-(y/K2)**alpha2))**beta2
def f(y,beta,q,alpha):
return hyp2f1(beta,(1-q)/alpha,1+(1-q)/alpha,y)
def f_imp(C,C0,beta,q,alpha,r,K):
return (C**(1-q)*f((C/K)**alpha,beta,q,alpha) - (C0)**(1-q)*f((C0/K)**alpha,beta,q,alpha))/(r*(1-q))
#find the K that maps to the time of the last data point
def timefit(C_curve):
return (f_imp(C_curve*K2,C02,beta2,q2,alpha2,r2,K2) - x[-1])**2
def timefit1(C_curve):
return (f_imp(C_curve*K2,C02,beta2,q2,alpha2,r2,K2) - (x[-1]+Forecast))**2
res = minimize_scalar(timefit, bounds=(0, 1), method='bounded')
C_plot = res.x
res1 = minimize_scalar(timefit1, bounds=(0, 1), method='bounded')
C_plot1 = res1.x
const = 0.9999 - C_plot
t_forecast1 = const*Forecast/100
c_forecast = np.arange(C02,C_plot*K2+1+Forecast,0.1)
alpha3 = alpha2
q3 = q2
beta3 = beta2
alpha2 = np.float(alpha2.value)
q2 = np.float(q2.value)
beta2 = np.float(beta2.value)
c_tc = ( ( q2 )**( -1 ) * ( alpha2 * beta2 + q2 ) )**( -1 * ( alpha2 )**( -1 ) )
C_tc = K2*c_tc
c_ts = ( 6 )**( -1 * ( alpha2 )**( -1 ) ) * ( ( ( 6 * ( alpha2 )**( 3 ) * ( beta2 )**( 3 ) + ( ( alpha2 )**( 2 ) * ( beta2 )**( 2 ) * ( -7 + 18 * q2 ) + ( q2 * ( 2 + ( -7 * q2 + 6 * ( q2 )**( 2 ) ) ) + 2 * alpha2 * beta2 * ( 1 + ( -7 * q2 + 9 * ( q2 )