SIR Real Data - COVID - Brazil Model Analyser.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Os dados\n",
    "\n",
    "Para coletar os dados de outros países, foram utilizadas as APIs dos sites:\n",
    "\n",
    "- [About Corona](https://about-corona.net/dashboard/)\n",
    "- [Europe RestFul](https://restcountries.eu/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import requests\n",
    "import pandas as pd\n",
    "\n",
    "covid_api = 'https://corona-api.com/countries/'\n",
    "rest_countries = 'https://restcountries.eu/rest/v2/alpha/'\n",
    "country = 'BR' # Alpha-2 ISO3166\n",
    "\n",
    "data_json =  requests.get(covid_api + country).json()\n",
    "country = requests.get(covid_api + country).json()\n",
    "\n",
    "N = country['data']['population']\n",
    "\n",
    "print(country['data']['name'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = [1,2,3,4]\n",
    "b = []\n",
    "b += a\n",
    "b += a\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_json[\"data\"][\"timeline\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Organizando os dados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from datetime import datetime\n",
    "\n",
    "df = pd.DataFrame(data_json['data']['timeline'])\n",
    "df = df.sort_values('date').reset_index()\n",
    "\n",
    "from datetime import datetime, timedelta\n",
    "df['date'] = [datetime.fromisoformat(f) for f in df['date']]\n",
    "df = df.drop_duplicates(subset='date', keep = 'last')\n",
    "\n",
    "# Criando o vetor de tempo\n",
    "first_date = df['date'].iloc[0]\n",
    "size_days = (df['date'].iloc[-1] - df['date'].iloc[0]).days\n",
    "date_vec = [first_date + timedelta(days=k) for k in range(size_days)]\n",
    "\n",
    "new_df = pd.DataFrame(date_vec, columns=['date'])\n",
    "new_df = pd.merge(new_df, df, how='left', on='date')\n",
    "new_df = new_df.drop(columns= ['index',  'updated_at', 'is_in_progress'])\n",
    "\n",
    "for col in new_df.columns[1:]:\n",
    "    new_df[col] = new_df[col].interpolate(method='polynomial', order=1)\n",
    "df = new_df.dropna()\n",
    "\n",
    "df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.iloc[135,:] = [df.iloc[135, 0]] + [None]*7\n",
    "df = df.interpolate(method ='linear', limit_direction ='forward')\n",
    "df = df.where(df.active != 0.0).dropna()\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"active\"].dtype==np.float64"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizando os dados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "# Creating the time vector for plotly\n",
    "datetime_64 = df[\"date\"].values\n",
    "ts = (datetime_64 - np.datetime64('1970-01-01T00:00:00Z')) / np.timedelta64(1, 's')\n",
    "time = [datetime.utcfromtimestamp(t) for t in ts]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import plotly.graph_objects as go\n",
    "\n",
    "first = -1\n",
    "\n",
    "figure = go.Figure()\n",
    "figure.add_trace(go.Scatter(\n",
    "                    x=time[:first], \n",
    "                    y=df[\"deaths\"].values[:first],\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    name = \"Mortes\",\n",
    "                    line = dict(color='#e91e63', width=3)))\n",
    "figure.add_trace(go.Scatter(\n",
    "                    x=time[:first], \n",
    "                    y=df[\"confirmed\"].values[:first],\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    name = \"Infectados acumulados\",\n",
    "                    line = dict(color='#03a9f4', width=3, dash=\"dot\")))\n",
    "figure.add_trace(go.Scatter(\n",
    "                    x=time[:first], \n",
    "                    y=df[\"recovered\"].values[:first],\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    name = \"Recuperados\",\n",
    "                    line = dict(color='#76d275', width=3)))\n",
    "figure.add_trace(go.Scatter(\n",
    "                    x=time[:first], \n",
    "                    y=df[\"active\"].values[:first],\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    name = \"Infectados ativos\",\n",
    "                    line = dict(color='#ffca28', width=3)))\n",
    "\n",
    "figure.update_layout(template='plotly_dark',\n",
    "                  xaxis_title='Data',\n",
    "                  yaxis_title='Indivíduos',\n",
    "                  title_text=\"Valores iniciais de infectados - \" + country['data']['name'])\n",
    "# 'ggplot2', 'seaborn', 'simple_white', 'plotly',\n",
    "# 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff',\n",
    "# 'ygridoff', 'gridon', 'none'\n",
    "\n",
    "figure.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verificando os dados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "dif_I = np.diff(df['active'])\n",
    "cum = []\n",
    "cum.append(dif_I[0])\n",
    "\n",
    "for k, i in enumerate(dif_I):\n",
    "    cum.append(cum[-1] + i)\n",
    "    \n",
    "cum = np.array(cum)\n",
    "cum += df['deaths'].to_numpy() + df['recovered'].to_numpy()\n",
    "\n",
    "print(\"Erro entre casos acumulados e valores de confirmados: {}\".format(\n",
    "        round(sum((cum - df['confirmed'].values)**2 / len(cum)),2) ) ) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Encontrando o início da epidemia\n",
    "\n",
    "Encontrando o início da epidemia para eliminar dados não informativos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Moment where cases are bigger than 2000...\n",
    "start_moment = np.argmax(df[\"active\"].to_numpy() >= 2000)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Criando os dados SIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "time_ref = time[start_moment:]\n",
    "I = df['active'].to_numpy()[start_moment:]\n",
    "R = df['recovered'].to_numpy()[start_moment:]\n",
    "M = df['deaths'].to_numpy()[start_moment:]\n",
    "S = N - R - I\n",
    "\n",
    "# Creating the time vector\n",
    "t = np.linspace(0, len(I)-1, len(I))\n",
    "\n",
    "Sd, Id, Md, Rd, td = S, I, M, R, t\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ro = (np.log(S[0]/N) - np.log(S[-1]/N))/(S[0]/N - S[-1]/N)\n",
    "print(\"Calculo do Ro baseado nos dados:\", Ro)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimando utilizando todos os dados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "from models import *\n",
    "\n",
    "dataset = dict(S=Sd, I=Id, R=Rd)\n",
    "\n",
    "# Create the model\n",
    "sir_model = ss.SIR(pop=N, focus=[\"S\", \"I\", \"R\"])\n",
    "\n",
    "# Adjust the parameters\n",
    "sir_model.fit(dataset, td,  \n",
    "              search_pop=True,\n",
    "              pop_sens=[0.00001, 0.05],\n",
    "              Ro_sens=[0.8, 15], \n",
    "              D_sens=[5, 40])\n",
    "\n",
    "\n",
    "# Predict the model\n",
    "initial = (Sd[0], Id[0], Rd[0])\n",
    "tsim = np.linspace(0, len(I)+120-1, len(I)+120)\n",
    "time_ref_sim = [time_ref[0] + timedelta(days=i) for i in tsim]\n",
    "\n",
    "sim_res = sir_model.predict(initial, tsim)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "print(\"Parâmetros estimados: \", sir_model.parameters)\n",
    "print(\"Suposto Ro: \", sir_model.parameters[0])\n",
    "print(\"Valores de beta: {}, e r: {}\".format(\n",
    "    sir_model.parameters[0]/sir_model.parameters[1], 1/sir_model.parameters[1]))\n",
    "print(\"Dias contaminados: \", sir_model.parameters[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "figure = go.Figure()\n",
    "figure.add_trace(go.Scatter(\n",
    "                    x=time_ref, \n",
    "                    y=I,\n",
    "                    name=\"Infectados\",\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    line = dict(color='#ffca28', width=3)))\n",
    "figure.add_trace(go.Scatter(\n",
    "                    x=time_ref,\n",
    "                    y=R,\n",
    "                    name=\"Recuperados\",\n",
    "                    mode=\"lines\",\n",
    "                    line_shape=\"spline\",\n",
    "                    line=dict(color='#76d275', width=3)))\n",
    "figure.add_trace(go.Scatter(\n",
    "                    x=time_ref_sim,\n",
    "                    y=sim_res[0],\n",
    "                    name=\"Suscetíveis - modelo\",\n",
    "                    mode=\"lines\",\n",
    "                    line_shape=\"spline\",\n",
    "                    line=dict(color=\"#1e88e5\", width=4, dash=\"dash\")))\n",
    "figure.add_trace(go.Scatter(\n",
    "                    x=time_ref_sim,\n",
    "                    y=sim_res[2],\n",
    "                    name=\"Recuperados - modelo\",\n",
    "                    mode=\"lines\",\n",
    "                    line_shape=\"spline\",\n",
    "                    line=dict(color=\"#43a047\", width=4, dash=\"dash\")))\n",
    "figure.add_trace(go.Scatter(\n",
    "                    x=time_ref_sim,\n",
    "                    y=sim_res[1],\n",
    "                    name=\"Infectados - modelo\",\n",
    "                    mode=\"lines\",\n",
    "                    line_shape=\"spline\",\n",
    "                    line=dict(color=\"#f4511e\", width=4, dash=\"dash\")))\n",
    "figure.update_layout(template='plotly_dark',\n",
    "                  xaxis_title='Data', yaxis_title='Indivíduos',\n",
    "                  title_text=\"Valores iniciais de infectados - \" + country['data']['name'])\n",
    "# 'ggplot2', 'seaborn', 'simple_white', 'plotly',\n",
    "# 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff',\n",
    "# 'ygridoff', 'gridon', 'none'\n",
    "figure.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monte Carlo\n",
    "\n",
    "\n",
    "Nesta parte, faremos um teste aumentando a quantidade de amostras de treinamento e prevendo o momento do pico da epidemia a medida que mais dias são utilizados para treinamento. Esse estudo vai possibilitar a análise da certeza da previsão do pico da epidemia antes desse acontecer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "saved_param = {'r':[], 'beta':[], 'pop':[]}\n",
    "saved_prediction = []\n",
    "\n",
    "start_day = 10\n",
    "pred_t = np.array(range(int(td[-1])))\n",
    "\n",
    "for i in range(start_day, len(I)):\n",
    "\n",
    "    dataset = dict(S=Sd[:i], I=Id[:i], R=Rd[:i])\n",
    "    td_ = td[:i]\n",
    "\n",
    "    # Create the model\n",
    "    sir_model = ss.SIR(pop=N, focus=[\"S\", \"I\", \"R\"], verbose = False)\n",
    "\n",
    "    # Adjust the parameters\n",
    "    sir_model.fit(dataset, td_, \n",
    "                  search_pop=True,\n",
    "                  pop_sens=[0.00001,0.05],\n",
    "                  Ro_sens=[0.8,15], \n",
    "                  D_sens=[5,40])\n",
    "    \n",
    "    saved_param['beta'].append(sir_model.parameters[0])\n",
    "    saved_param['r'].append(sir_model.parameters[1])\n",
    "    saved_param['pop'].append(sir_model.parameters[2])\n",
    "    \n",
    "    saved_prediction.append(sir_model.predict((Sd[0],Id[0], Rd[0]), pred_t))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "with open(\"./China_II_mc_runs.pickle\", \"wb\") as handle:\n",
    "    pickle.dump({\"pars\":saved_param, \"pred\":saved_prediction}, handle)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Análise do uso do sistema de saúde\n",
    "\n",
    "Nesta análise, mostramos o erro percentual do quanto antes do pico, conseguimos prever a quantidade de pessoas que realmente serão identificadas como infectadas, uma vez que $R(\\infty)$ é a quantidade de pessoas recuperadas totais, daquelas que foram noficadas como infectadas no sistema de saúde. Desta forma segue o erro proporcional do erro a medida em que novos dados diários foram incluidos no modelo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "x = range(start_day, len(I))\n",
    "usage_error = [ 100 * abs(p*N - Rd[-1]) / Rd[-1] for p in saved_param['pop']]\n",
    "\n",
    "fig9 = go.Figure()\n",
    "\n",
    "fig9.add_trace(go.Scatter(\n",
    "                    x=td[start_day:peak_pos], \n",
    "                    y=usage_error[:peak_pos-start_day],\n",
    "                    mode='lines',\n",
    "                    name='Antes do pico',\n",
    "                    line_shape='spline',\n",
    "                    line = dict(color='royalblue', width=3),\n",
    "                    hovertemplate=\"ε(%) = %{y:.0f}, <br> com %{x:.0f} dias de dados.\"))\n",
    "fig9.add_trace(go.Scatter(\n",
    "                    x=td[peak_pos:-1], \n",
    "                    y=usage_error[peak_pos-start_day:-1],\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    name='Depois do pico',\n",
    "                    line = dict(color='royalblue', width=3, dash='dot'),\n",
    "                    hovertemplate=\"ε(%) = %{y:.0f}, <br> com %{x:.0f} dias de dados.\"))\n",
    "fig9.add_trace(go.Scatter(\n",
    "                    mode=\"markers\", x=[peak_pos], y=[usage_error[peak_pos-start_day]], \n",
    "                    marker_symbol=\"hexagram-dot\", name=\"Momento do pico\",\n",
    "                    marker_line_color=\"midnightblue\", marker_color=\"lightskyblue\", \n",
    "                    marker_line_width=2, marker_size=15, \n",
    "                    hovertemplate=\"Pico no dia %{x}, com um ε(%) = %{y:.0f}.\"))\n",
    "fig9.add_trace(go.Scatter(\n",
    "                    mode=\"markers\", x=[td[-1]], y=[usage_error[-1]], \n",
    "                    marker_symbol=\"triangle-right-dot\", name=\"Valor atual\",\n",
    "                    marker_line_color=\"#a00037\", marker_color=\"#ff5c8d\", \n",
    "                    marker_line_width=2, marker_size=15, \n",
    "                    hovertemplate=\"No dia %{x} da epidemia, <br> convergindo para ε(%) = %{y:.4f}.\"))\n",
    "fig9.update_layout(template='xgridoff', yaxis_range=[-1,100],\n",
    "                  legend_orientation=\"h\", legend=dict(x=0.10, y=1.05),\n",
    "                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',\n",
    "                  yaxis_title='Erro (%) do R(∞) - notificados',\n",
    "                  title_text=\"Erro da previsão do uso do sistema de saúde - 🇨🇳 \" + country['data']['name'])\n",
    "\n",
    "fig9.show(renderer=\"png\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizando as previsões de $I(t)$\n",
    "\n",
    "Vamos analisar as previsões quando somente os dados antes do pico são fornecidos ao modelo, e as previões utilizando os dados após o pico:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "visual_peak = 30\n",
    "\n",
    "p2 = figure(plot_height=500,\n",
    "           plot_width=600, \n",
    "           tools=\"\", \n",
    "           toolbar_location=None,\n",
    "           title=\"Evolução do COVID - \" + country['data']['name'])\n",
    "\n",
    "# Preparando o estilo\n",
    "p2.grid.grid_line_alpha = 0\n",
    "p2.ygrid.band_fill_color = \"olive\"\n",
    "p2.ygrid.band_fill_alpha = 0.1\n",
    "p2.yaxis.axis_label = \"Indivíduos\"\n",
    "p2.xaxis.axis_label = \"Dias\"\n",
    "\n",
    "# Incluindo as curvas\n",
    "p2.line(td, Id,\n",
    "       legend_label=\"Infectados\", \n",
    "       line_cap=\"round\", line_width=3, color=\"#c62828\")\n",
    "\n",
    "for data in saved_prediction[visual_peak-start_day:]:\n",
    "    p2.line(pred_t[:len(td)], -data[1][:len(td)],\n",
    "           legend_label=\"Previsão Infectados - Depois do pico\", \n",
    "           line_cap=\"round\", line_dash=\"dashed\", line_width=4, color=\"#ffa000\", line_alpha = 0.1)\n",
    "\n",
    "for data in saved_prediction[:visual_peak-start_day]:\n",
    "    p2.line(pred_t[:len(td)], -data[1][:len(td)],\n",
    "           legend_label=\"Previsão Infectados - Antes do pico\", \n",
    "           line_cap=\"round\", line_dash=\"dashed\", line_width=4, color=\"#42a5f5\", line_alpha = 0.3)\n",
    "\n",
    "# Colocando as legendas\n",
    "p2.legend.click_policy=\"hide\"\n",
    "p2.legend.location = \"top_right\"\n",
    "\n",
    "show(p2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizando os ajustes dos grupos\n",
    "\n",
    "\n",
    "Aqui vamos analisar as previsões obtidas para cada grupo $S(t)$, $I(t)$ e $R(t)$, a medida que mais dias foram informados ao modelo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "p3 = figure(plot_height=350,\n",
    "           plot_width=600, \n",
    "           tools=\"\", \n",
    "           toolbar_location=None,\n",
    "           title=\"Evolução do COVID - I(t) e R(t) - \" + country['data']['name'])\n",
    "\n",
    "p4 = figure(plot_height=350,\n",
    "           plot_width=600, \n",
    "           tools=\"\", \n",
    "           toolbar_location=None,\n",
    "           title=\"Evolução do COVID - S(t) - \" + country['data']['name'])\n",
    "\n",
    "plot_all = True\n",
    "\n",
    "# Preparando o estilo\n",
    "p3.grid.grid_line_alpha = 0\n",
    "p3.ygrid.band_fill_color = \"olive\"\n",
    "p3.ygrid.band_fill_alpha = 0.1\n",
    "p3.yaxis.axis_label = \"Indivíduos\"\n",
    "p3.xaxis.axis_label = \"Dias\"\n",
    "\n",
    "p4.grid.grid_line_alpha = 0\n",
    "p4.ygrid.band_fill_color = \"olive\"\n",
    "p4.ygrid.band_fill_alpha = 0.1\n",
    "p4.yaxis.axis_label = \"Indivíduos\"\n",
    "p4.xaxis.axis_label = \"Dias\"\n",
    "\n",
    "# Incluindo as curvas\n",
    "for data in saved_prediction[20:]:\n",
    "    p3.line(pred_t, data[1],\n",
    "           legend_label=\"Previsão Infectados\", \n",
    "           line_cap=\"round\", line_dash = 'dashed', \n",
    "           line_width=4, color=\"#42a5f5\", line_alpha = 0.1)\n",
    "    \n",
    "    p3.line(pred_t, data[2],\n",
    "           legend_label=\"Previsão Recuperados\", \n",
    "           line_cap=\"round\", line_dash = 'dashed', line_width=4, \n",
    "           color=\"#9c27b0\", line_alpha = 0.07)\n",
    "    \n",
    "    \n",
    "p3.line(td, Id,\n",
    "       legend_label=\"Infectados\", \n",
    "       line_cap=\"round\", line_width=3, color=\"#005cb2\")\n",
    "\n",
    "p3.line(td, Rd,\n",
    "           legend_label=\"Recuperados\", \n",
    "           line_cap=\"round\", line_width=3, color=\"#5e35b1\")\n",
    "\n",
    "if plot_all:\n",
    "    for data in saved_prediction:\n",
    "        p4.line(pred_t, data[0] + N*(1-saved_param['pop'][-1]),\n",
    "               legend_label=\"Previsão Suscetiveis\", \n",
    "               line_cap=\"round\", line_dash = 'dashed', \n",
    "               line_width=4, color=\"#ff5722\", line_alpha = 0.07)\n",
    "    p4.line(td, Sd,\n",
    "           legend_label=\"Suscetiveis\", \n",
    "           line_cap=\"round\", line_width=3, color=\"#b71c1c\")\n",
    "\n",
    "# Colocando as legendas\n",
    "p3.legend.click_policy=\"hide\"\n",
    "p3.legend.location = \"top_left\"\n",
    "\n",
    "p4.legend.click_policy=\"hide\"\n",
    "p4.legend.location = \"top_right\"\n",
    "\n",
    "show(column(p3,p4))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Análise de variação do $R_0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = range(start_day, len(I))\n",
    "beta_norm = [b*p for b,p in zip(saved_param['beta'],saved_param['pop'])]\n",
    "Ro = [b*p/r for b,r,p in zip(saved_param['beta'],saved_param['r'],saved_param['pop'])]\n",
    "\n",
    "\n",
    "fig1 = go.Figure()\n",
    "\n",
    "fig1.add_trace(go.Scatter(\n",
    "                    x=td[start_day:peak_pos], \n",
    "                    y=beta_norm[:peak_pos-start_day],\n",
    "                    mode='lines',\n",
    "                    name='Antes do pico',\n",
    "                    line_shape='spline',\n",
    "                    line = dict(color='royalblue', width=3),\n",
    "                    hovertemplate=\"β = %{y:.4f}, <br> com %{x:.0f} dias de dados.\"))\n",
    "fig1.add_trace(go.Scatter(\n",
    "                    x=td[peak_pos:-1], \n",
    "                    y=beta_norm[peak_pos-start_day:-1],\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    name='Depois do pico',\n",
    "                    line = dict(color='royalblue', width=3, dash='dot'),\n",
    "                    hovertemplate=\"β = %{y:.4f}, <br> com %{x:.0f} dias de dados.\"))\n",
    "fig1.add_trace(go.Scatter(\n",
    "                    mode=\"markers\", x=[peak_pos], y=[beta_norm[peak_pos-start_day]], \n",
    "                    marker_symbol=\"hexagram-dot\", name=\"Momento do pico\",\n",
    "                    marker_line_color=\"midnightblue\", marker_color=\"lightskyblue\", \n",
    "                    marker_line_width=2, marker_size=15, \n",
    "                    hovertemplate=\"Pico no dia %{x}, com um β = %{y:.4f}.\"))\n",
    "fig1.add_trace(go.Scatter(\n",
    "                    mode=\"markers\", x=[td[-1]], y=[beta_norm[-1]], \n",
    "                    marker_symbol=\"triangle-right-dot\", name=\"Valor atual\",\n",
    "                    marker_line_color=\"#a00037\", marker_color=\"#ff5c8d\", \n",
    "                    marker_line_width=2, marker_size=15, \n",
    "                    hovertemplate=\"No dia %{x} da epidemia, <br> convergindo para β = %{y:.4f}.\"))\n",
    "fig1.update_layout(template='xgridoff',\n",
    "                  legend_orientation=\"h\", legend=dict(x=0.10, y=1.05),\n",
    "                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',\n",
    "                  yaxis_title='Parâmetro β',\n",
    "                  title_text=\"Estimativa do parâmetro β - 🇨🇳 \" + country['data']['name'])\n",
    "\n",
    "fig1.show(renderer=\"png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fig2 = go.Figure()\n",
    "\n",
    "fig2.add_trace(go.Scatter(\n",
    "                    x=td[start_day:peak_pos], \n",
    "                    y=saved_param['r'][:peak_pos-start_day],\n",
    "                    mode='lines',\n",
    "                    name='Antes do pico',\n",
    "                    line_shape='spline',\n",
    "                    line = dict(color='royalblue', width=3),\n",
    "                    hovertemplate=\"r = %{y:.4f}, <br> com %{x:.0f} dias de dados.\"))\n",
    "fig2.add_trace(go.Scatter(\n",
    "                    x=td[peak_pos:-1], \n",
    "                    y=saved_param['r'][peak_pos-start_day:-1],\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    name='Depois do pico',\n",
    "                    line = dict(color='royalblue', width=3, dash='dot'),\n",
    "                    hovertemplate=\"r = %{y:.4f}, <br> com %{x:.0f} dias de dados.\"))\n",
    "fig2.add_trace(go.Scatter(\n",
    "                    mode=\"markers\", x=[peak_pos], y=[saved_param['r'][peak_pos-start_day]], \n",
    "                    marker_symbol=\"hexagram-dot\", name=\"Momento do pico\",\n",
    "                    marker_line_color=\"midnightblue\", marker_color=\"lightskyblue\", \n",
    "                    marker_line_width=2, marker_size=15, \n",
    "                    hovertemplate=\"Pico no dia %{x}, com um r = %{y:.4f}.\"))\n",
    "fig2.add_trace(go.Scatter(\n",
    "                    mode=\"markers\", x=[td[-1]], y=[saved_param['r'][-1]], \n",
    "                    marker_symbol=\"triangle-right-dot\", name=\"Valor atual\",\n",
    "                    marker_line_color=\"#a00037\", marker_color=\"#ff5c8d\", \n",
    "                    marker_line_width=2, marker_size=15, \n",
    "                    hovertemplate=\"No dia %{x} da epidemia, <br> convergindo para r = %{y:.4f}.\"))\n",
    "fig2.update_layout(template='xgridoff',\n",
    "                  legend_orientation=\"h\", legend=dict(x=0.07, y=1.06),\n",
    "                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',\n",
    "                  yaxis_title='Parâmetro r',\n",
    "                  title_text=\"Estimativa do parâmetro r - 🇨🇳 \" + country['data']['name'])\n",
    "\n",
    "fig2.show(renderer=\"png\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fig3 = go.Figure()\n",
    "\n",
    "fig3.add_trace(go.Scatter(\n",
    "                    x=td[start_day:peak_pos], \n",
    "                    y=Ro[:peak_pos-start_day],\n",
    "                    mode='lines',\n",
    "                    name='Antes do pico',\n",
    "                    line_shape='spline',\n",
    "                    line = dict(color='royalblue', width=3),\n",
    "                    hovertemplate=\"r = %{y:.4f}, <br> com %{x:.0f} dias de dados.\"))\n",
    "fig3.add_trace(go.Scatter(\n",
    "                    x=td[peak_pos:-1], \n",
    "                    y=Ro[peak_pos-start_day:-1],\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    name='Depois do pico',\n",
    "                    line = dict(color='royalblue', width=3, dash='dot'),\n",
    "                    hovertemplate=\"Ro = %{y:.4f}, <br> com %{x:.0f} dias de dados.\"))\n",
    "fig3.add_trace(go.Scatter(\n",
    "                    mode=\"markers\", x=[peak_pos], y=[Ro[peak_pos-start_day]], \n",
    "                    marker_symbol=\"hexagram-dot\", name=\"Momento do pico\",\n",
    "                    marker_line_color=\"midnightblue\", marker_color=\"lightskyblue\", \n",
    "                    marker_line_width=2, marker_size=15, \n",
    "                    hovertemplate=\"Pico no dia %{x}, com um Ro = %{y:.4f}.\"))\n",
    "fig3.add_trace(go.Scatter(\n",
    "                    mode=\"markers\", x=[td[-1]], y=[Ro[-1]], \n",
    "                    marker_symbol=\"triangle-right-dot\", name=\"Valor atual\",\n",
    "                    marker_line_color=\"#a00037\", marker_color=\"#ff5c8d\", \n",
    "                    marker_line_width=2, marker_size=15, \n",
    "                    hovertemplate=\"No dia %{x} da epidemia, <br> convergindo para Ro = %{y:.4f}.\"))\n",
    "fig3.update_layout(template='xgridoff',\n",
    "                  legend_orientation=\"h\", legend=dict(x=0.07, y=1.06),\n",
    "                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',\n",
    "                  yaxis_title='Parâmetro Ro',\n",
    "                  title_text=\"Estimativa do parâmetro Ro - 🇨🇳 \" + country['data']['name'])\n",
    "\n",
    "fig3.show(renderer=\"png\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Análise de confiança"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from PyAstronomy import pyasl\n",
    "\n",
    "dI = np.gradient(pyasl.smooth(I, 13, \"hamming\"))\n",
    "t = np.linspace(0, len(dI), len(dI))\n",
    "\n",
    "signal = np.array([di >= 0 for di in dI[::-1]])\n",
    "\n",
    "\n",
    "import plotly.graph_objects as go\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(\n",
    "                    x=t[signal], \n",
    "                    y=dI[::-1][signal],\n",
    "                    mode='lines',\n",
    "                    name='Antes do pico - Derivada positiva',\n",
    "                    line_shape='spline',\n",
    "                    line = dict(color='#512da8', width=3)))\n",
    "fig.add_trace(go.Scatter(\n",
    "                    x=t[~signal], \n",
    "                    y=dI[::-1][~signal],\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    name='Depois do pico - Derivada negativa',\n",
    "                    line = dict(color='#d8576b', width=3)))\n",
    "\n",
    "fig.update_layout(template='xgridoff',\n",
    "                  xaxis_title='Dias (invertido)',\n",
    "                  yaxis_title='Taxa de variação de indivíduos',\n",
    "                  title_text=\"Derivada da curva de Infectados - \" + country['data']['name'])\n",
    "# 'ggplot2', 'seaborn', 'simple_white', 'plotly',\n",
    "# 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff',\n",
    "# 'ygridoff', 'gridon', 'none'\n",
    "\n",
    "fig.show(renderer=\"png\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "peak_pos = np.argmax(~signal[::-1])\n",
    "print(\"O pico da epidemia acontece no dia\", peak_pos)\n",
    "\n",
    "final_day = int(td[-1])\n",
    "\n",
    "estimated_peaks = []\n",
    "for data in saved_prediction:\n",
    "    dI = np.gradient(data[1][:final_day])\n",
    "    signal_pred = np.array([di >= 0 for di in dI[::-1]])\n",
    "    estimated_peaks.append(len(Id) - np.argmax(signal_pred))\n",
    "estimated_peaks = np.array(estimated_peaks)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import plotly.graph_objects as go\n",
    "\n",
    "peak_error = np.abs(estimated_peaks - peak_pos)\n",
    "\n",
    "fig1 = go.Figure()\n",
    "fig1.add_trace(go.Scatter(\n",
    "                    x=td[start_day:peak_pos], \n",
    "                    y=peak_error[:peak_pos-start_day],\n",
    "                    mode='lines',\n",
    "                    name='Antes do pico',\n",
    "                    line_shape='spline',\n",
    "                    line = dict(color='royalblue', width=3),\n",
    "                    hovertemplate=\"Erro de %{y} dias, <br> com %{x:.0f} dias de dados.\"))\n",
    "fig1.add_trace(go.Scatter(\n",
    "                    x=td[peak_pos:], \n",
    "                    y=peak_error[peak_pos-start_day:],\n",
    "                    mode='lines',\n",
    "                    line_shape='spline',\n",
    "                    name='Depois do pico',\n",
    "                    line = dict(color='royalblue', width=3, dash='dot'),\n",
    "                    hovertemplate=\"Erro de %{y} dias, <br> com %{x:.0f} dias de dados.\"))\n",
    "fig1.add_trace(go.Scatter(\n",
    "                    mode=\"markers\", x=[peak_pos], y=[peak_error[peak_pos-start_day]], \n",
    "                    marker_symbol=\"hexagram-dot\", name=\"Momento do pico\",\n",
    "                    marker_line_color=\"midnightblue\", marker_color=\"lightskyblue\", \n",
    "                    marker_line_width=2, marker_size=15, \n",
    "                    hovertemplate=\"Pico no dia %{x} depois do começo da epidemia.\"))\n",
    "\n",
    "fig1.update_layout(yaxis_range=[-1,31], \n",
    "                  template='xgridoff',\n",
    "                  legend_orientation=\"h\", legend=dict(x=0.20, y=1.0),\n",
    "                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',\n",
    "                  yaxis_title='Erro (em dias) da estimativa do pico',\n",
    "                  title_text=\"Erro da estimativa do pico - \" + country['data']['name'])\n",
    "# 'ggplot2', 'seaborn', 'simple_white', 'plotly',\n",
    "# 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff',\n",
    "# 'ygridoff', 'gridon', 'none'\n",
    "\n",
    "fig1.show(renderer=\"png\")\n"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}