#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

global R2_generalc
global stage
stage = ' '
global lang
lang = 'Por'
global index
index = 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 foi projetado para gerar ajustes matemáticos das curvas epidêmicas da Covid-19 para países, bem como para estados/cidades do Brasil e dos EUA. Para uma descrição dos modelos, vide tutorial <a href="http://fisica.ufpr.br/modinterv/Tutorial_app_english.pdf"> aqui</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','RM fit','GRM fit','BLM fit'],'Por':['ajuste q-exponencial','ajuste MR','ajuste MRG','ajuste MLB']}
methods_name1 = {'Eng':['q-exponential fit','RM fit','GRM fit','BLM fit'],'Por':['Ajuste q-exponencial','Ajuste MR','Ajuste MRG','Ajuste MLB']}
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.'}

#nbi:hide_in

def RMPlot(Country,Type,name2):
  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_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)]
  active = confirmed - recovered - deaths
 
 
  if Type == 'Deaths':
    kind1 = deaths
  if Type == 'Cases':
    kind1 = confirmed

  kind1 = kind1[0:len(kind1)-1-index]
 
  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):

  name = Country
  t_forecast = 1 + Forecast/100
 
  t = np.arange(0, int(t_forecast*len(kind1))+1, 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 ) ) ) ))

  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)

  dic = {'Linear':'linear', 'Logarithmic':'log'}

  if Type=='Deaths':
    tipo = 'Mortes'
    tipo1 = 'morte'
    tipo2 = 'a primeira morte'
  if Type=='Cases':
    tipo = 'Casos'
    tipo1 = 'caso'
    tipo2 = 'o primeiro caso'

  #global lang
  lang = html_widget_lang.value
  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=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$2020 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$')
    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):
 
  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_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)]
  active = confirmed - recovered - deaths
 
 
  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):
        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
        
          res = minimize_scalar(timefit, bounds=(0, 1), method='bounded')
          C_plot = res.x

          const = 0.9999 - C_plot
          t_forecast1 = const*Forecast/100
          c_forecast = np.arange(C02,(C_plot+t_forecast1)*K2+1,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)

          dic = {'Linear':'linear', 'Logarithmic':'log'}

          #global lang
          lang = html_widget_lang.value

          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=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=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$')
            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):
  with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    name = Country
    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 

    kind1 = kind1[0:len(kind1)-1-index]

    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)
        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):
        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
        
          res = minimize_scalar(timefit, bounds=(0, 1), method='bounded')
          C_plot = res.x

          const = 0.9999 - C_plot
          t_forecast1 = const*Forecast/100
          c_forecast = np.arange(C02,(C_plot+t_forecast1)*K2+1,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