#nbi:hide_in
# Cell to define ALL the widgets that appear in the code

# Preamble
import pandas
import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.style.use('ggplot')
!pip install -q lmfit
import lmfit
!pip install -q numdifftools
import numdifftools
from timeit import default_timer as timer
from lmfit import Parameters, Minimizer, report_fit
import ipywidgets as widgets
from ipywidgets import interact, Dropdown
from scipy.optimize import curve_fit, minimize_scalar, fsolve
from scipy.special import hyp2f1
from scipy.interpolate import InterpolatedUnivariateSpline
from scipy.integrate import solve_ivp
from IPython.display import display
from IPython.display import FileLink, FileLinks
from contextlib import suppress
from IPython.display import HTML
import base64
!pip install -q fpdf
from fpdf import FPDF
import time
import warnings
warnings.filterwarnings("ignore")
from datetime import datetime, timezone, timedelta
!pip install -q sigfig
from sigfig import round
!pip install -q tsmoothie
from tsmoothie.smoother import LowessSmoother
from scipy.signal import find_peaks
from scipy.integrate import odeint

global R2_generalc
global stage
stage = ' '
global lang
lang = 'Por'
global index
index = 0
global data_segunda_onda
data_segunda_onda = '0'
global date_first
date_first = '0'

button_pt = widgets.Button(description='Português')
button_en = widgets.Button(description='English')
display(widgets.HBox([button_pt, button_en]))
html_display_widget_intro = widgets.HTML(
    value='<hr><p>Este aplicativo faz ajustes matemáticos das curvas epidêmicas de países, bem como de estados e cidades do Brasil e dos EUA.</p> <p> Esta versão implementa um <a href="https://doi.org/10.1101/2021.01.31.21250867"> modelo com duas ondas</a> epidêmicas. Os ajustes são mais complexos e, portanto, mais demorados do que os da <a href="http://fisica.ufpr.br/modinterv/onewave.html"> versão anterior</a>, com o <a href="https://www.nature.com/articles/s41598-021-84165-1"> modelo de uma onda</a>. </p> <hr> <p>O aplicativo foi desenvolvido pela Rede de Pesquisa Modinterv-Covid19 formada por pesquisadores das Universidades Federais do Paraná, Pernambuco e Sergipe. Este aplicativo está disponibilizado sob licença <a href="https://creativecommons.org/licenses/by-nc-nd/4.0/">Creative Commons CC BY-NC-ND 4.0</a>.</p> <hr>')
display(html_display_widget_intro)
html_widget_lang = widgets.HTML(value='Por')

data_e_hora_atuais = datetime.now()
diferenca = timedelta(hours=-3)
fuso_horario = timezone(diferenca)
data_e_hora_sao_paulo = data_e_hora_atuais.astimezone(fuso_horario)
data_e_hora_sao_paulo_em_texto = data_e_hora_sao_paulo.strftime('%d/%m/%Y %H:%M')
#data_e_hora_sao_paulo_em_texto = '21/08/2020 14:10'
data_e_hora_sao_paulo_em_texto_en = data_e_hora_sao_paulo.strftime('%m/%d/%Y %H:%M')
data_e_hora_papers = data_e_hora_sao_paulo.strftime('%b %e, %Y %H:%M')

time_dic = {'Eng':data_e_hora_sao_paulo_em_texto_en,'Por':data_e_hora_sao_paulo_em_texto}

fase_en = ['increasing acceleration', 'decreasing acceleration', 'increasing deceleration', 'transition to saturation', 'saturation']
fase_br = ['aceleração crescente', 'aceleração decrescente','desaceleração crescente', 'transição para saturação', 'saturação'] 
fase = {'Eng':fase_en,'Por':fase_br}

estagio = {'Eng':'Epidemic stage:','Por':'Estágio da epidemia:'}

label_x = {'Eng':{'Deaths':'Days since first death','Cases':'Days since first case'},'Por':{'Deaths':'Dias desde o primeiro óbito','Cases':'Dias desde o primeiro caso'}}
label_y = {'Eng':{'Deaths':'Accumulated deaths','Cases':'Accumulated cases'},'Por':{'Deaths':'Óbitos acumulados','Cases':'Casos acumulados'}}
label_y_daily = {'Eng':{'Deaths':'Daily deaths','Cases':'Daily cases'},'Por':{'Deaths':'Óbitos diários','Cases':'Casos diários'}}

subs = {'Eng':{'Deaths':'Deaths','Cases':'Cases'},'Por':{'Deaths':'Mortes','Cases':'Casos'}}
methods_name = {'Eng':['q-exponential fit (1-Wave)','RM fit (1-Wave)','GRM fit (1-Wave)','BLM fit (1-Wave)','BLM fit (2-Waves)'],'Por':['ajuste q-exponencial (1-Onda)','ajuste MR (1-Onda)','ajuste MRG (1-Onda)','ajuste MLB (1-Onda)','ajuste MLB (2-Ondas)']}
methods_name1 = {'Eng':['q-exponential fit (1-Wave)','RM fit (1-Wave)','GRM fit (1-Wave)','BLM fit (1-Wave)','BLM fit (2-Waves)'],'Por':['Ajuste q-exponencial (1-Onda)','Ajuste MR (1-Onda)','Ajuste MRG (1-Onda)','Ajuste MLB (1-Onda)','Ajuste MLB (2-Ondas)']}
R2_message = {'Eng':'The methods did not provide reliable fits for this location, possibly because the curve has unusual behavior. The least-worst fit can be seen by clicking the checkbox below.','Por':'Os métodos não produziram ajustes confiáveis para esta localidade, possivelmente porque a curva tem um comportamento anômalo. Não obstante, o melhor ajuste dentre os métodos pode ser visto clicando na caixa de seleção abaixo.'}
dic = {'Linear':'linear', 'Logarithmic':'log'}
dic1 = {'Linear':'linear', 'Logarithmic':'log'}
wave_date = {'Eng':'Start 2nd wave: ', 'Por':r'Início $2^{{\rm a}}$ onda: '}
first_death = {'Eng':{'Deaths':'1st death: ','Cases':'1st case: '}, 'Por':{'Deaths':r'$1^{{\rm o}}$ óbito: ','Cases':r'$1^{{\rm o}}$ caso: '}}
error_message = {'Eng':'Warning: fit not entirely reliable because of large errors in some of the parameters.','Por':'Aviso: ajuste não totalmente confiável por causa de erros grandes em algum dos parâmetros.'}
first_plateau = {'Eng':r'Plateau 1st wave: $K_1$= ','Por':r'Platô $1^{{\rm a}}$ onda: $K_1$= '}
second_plateau = {'Eng':r'Plateau 2nd wave: $K_2$= ','Por':r'Platô $2^{{\rm a}}$ onda: $K_2$= '}
html_loading_message = {'Eng':'<html><head><style>body {background-color: #ffffff;}h1 {background: url("http://i.stack.imgur.com/tQTRW.gif");background-repeat:no-repeat;background-size:120px 30px;font-size:14pt}</style></head><body><h1>Computing...</h1></body></html>','Por':'<html><head><style>body {background-color: #ffffff;}h1 {background: url("http://i.stack.imgur.com/tQTRW.gif");background-repeat:no-repeat;background-size:120px 30px;font-size:14pt}</style></head><body><h1>Calculando...</h1></body></html>'}
first_peak = {'Eng':'Peak 1st wave: ','Por':r'Pico $1^{{\rm a}}$ onda: '}
second_peak = {'Eng':'Peak 2nd wave: ','Por':r'Pico $2^{{\rm a}}$ onda: '}
date_last_message = {'Eng':'Date of last point: ', 'Por':'Data do último ponto: '}
forecast_message = {'Eng':{'7':'Forecast 7d: ','14':'Forecast 14d: ','21':'Forecast 21d: ','28':'Forecast 28d: '}, 'Por':{'7':'Previsão 7d: ','14':'Previsão 14d: ','21':'Previsão 21d: ','28':'Previsão 28d: '}}

#nbi:hide_in
# Importing the training database for the AI and training the AI.
url = 'https://gist.githubusercontent.com/ArthurAraujoBrum/12d867cc198f1ac2e3730ad067eed46d/raw/0cddaf965f421b1c9a223381d29b1a43206667e8/training_base.csv'
df1 = pandas.read_csv(url, error_bad_lines=False)
feature_names = ['0p', '5p', '10p', '15p', '20p', '25p', '30p', '35p', '40p', '45p','50p', '55p', '60p', '65p', '70p', '75p', '80p', '85p', '90p', '95p']
X = df1[feature_names]
y = df1['2-wave']

lang = html_widget_lang.value
html_display_widget_loading = widgets.HTML(value='<html><head><style>body {background-color: #ffffff;}h1 {background: url("http://i.stack.imgur.com/tQTRW.gif");background-repeat:no-repeat;background-size:120px 30px;font-size:14pt}</style></head><body><h1>Carregando...</h1></body></html>')
display(html_display_widget_loading)

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

from sklearn.tree import DecisionTreeClassifier
clf = DecisionTreeClassifier().fit(X_train, y_train)
#print('Accuracy of Decision Tree classifier on training set: {:.2f}'
     #.format(clf.score(X_train, y_train)))
#print('Accuracy of Decision Tree classifier on test set: {:.2f}'
     #.format(clf.score(X_test, y_test)))

from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier()
knn.fit(X_train, y_train)
#print('Accuracy of K-NN classifier on training set: {:.2f}'
     #.format(knn.score(X_train, y_train)))
#print('Accuracy of K-NN classifier on test set: {:.2f}'
     #.format(knn.score(X_test, y_test)))

from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import make_classification
RF = RandomForestClassifier(n_estimators=1000, random_state=0)
RF.fit(X_train, y_train)
#print('Accuracy of RF classifier on training set: {:.2f}'
     #.format(RF.score(X_train, y_train)))
#print('Accuracy of RF classifier on test set: {:.2f}'
     #.format(RF.score(X_test, y_test)))

html_display_widget_loading.value=' '

#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 )