SIR Real Data - COVID - China Model with Cabral LS.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 = 'CN' # 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": "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.head()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizando os dados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from bokeh.models   import Legend, ColumnDataSource, RangeTool, LinearAxis, Range1d, HoverTool\n",
    "from bokeh.palettes import brewer, Inferno256\n",
    "from bokeh.plotting import figure, show\n",
    "from bokeh.layouts  import column\n",
    "from bokeh.io       import output_notebook\n",
    "\n",
    "output_notebook()\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "# Criando os valores para legenda no plot\n",
    "year =  [str(int(d.year)) for d in df['date'] ]\n",
    "month = [(\"0\"+str(int(d.month)))[-2:] for d in df['date'] ]\n",
    "day =   [(\"0\"+str(int(d.day)))[-2:] for d in df['date'] ]\n",
    "\n",
    "# Criando a fonte de dados\n",
    "source = ColumnDataSource(data={\n",
    "    'Data'       : df['date'].values,\n",
    "    'd': day, 'm': month, 'y': year,\n",
    "    'Infectados Acc' : df['confirmed'].values,\n",
    "    'Mortes'     : df['deaths'].values,\n",
    "    'Ativo'     : df['active'].values,\n",
    "    'Recuperados': df['recovered'].values\n",
    "})\n",
    "\n",
    "\n",
    "# Criando a figura\n",
    "p = figure(plot_height=500,\n",
    "           plot_width=600,\n",
    "           x_axis_type=\"datetime\", \n",
    "           tools=\"\", \n",
    "           #y_axis_type=\"log\",\n",
    "           toolbar_location=None,\n",
    "           title=\"Evolução do COVID - \" + country['data']['name'])\n",
    "\n",
    "# Preparando o estilo\n",
    "p.grid.grid_line_alpha = 0\n",
    "p.ygrid.band_fill_color = \"olive\"\n",
    "p.ygrid.band_fill_alpha = 0.1\n",
    "p.yaxis.axis_label = \"Indivíduos\"\n",
    "p.xaxis.axis_label = \"Dias\"\n",
    "\n",
    "# Incluindo as curvas\n",
    "i_p = p.line(x='Data', y='Ativo',\n",
    "             legend_label=\"Infectados\", \n",
    "             line_cap=\"round\", line_width=5, color=\"#c62828\", source=source)\n",
    "m_p = p.line(x='Data', y='Mortes',\n",
    "             legend_label=\"Mortes\",\n",
    "             line_cap=\"round\", line_width=5, color=\"#512da8\", source=source)\n",
    "c_p = p.line(x='Data', y='Infectados Acc',\n",
    "             legend_label=\"Infectados Acc\",\n",
    "             line_cap=\"round\", line_width=5, color=\"#0288d1\", source=source)\n",
    "r_p = p.line(x='Data', y='Recuperados',\n",
    "             legend_label=\"Recuperados\",\n",
    "             line_cap=\"round\", line_width=5, color=\"#388e3c\", source=source)\n",
    "\n",
    "# Colocando as legendas\n",
    "p.legend.click_policy=\"hide\"\n",
    "# p.legend.location = \"top_left\"\n",
    "p.legend.location = \"top_left\"\n",
    "\n",
    "# Incluindo a ferramenta de hover\n",
    "p.add_tools(HoverTool(\n",
    "    tooltips=[\n",
    "        ( 'Indivíduos', '$y{i}'),\n",
    "        ( 'Data',       '@d/@m/@y' ),\n",
    "    ],\n",
    "    renderers=[\n",
    "        m_p, i_p, c_p, r_p\n",
    "    ]\n",
    "))\n",
    "\n",
    "show(p)\n"
   ]
  },
  {
   "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": [
    "## Criando os dados SIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "I = df['active'].to_numpy()\n",
    "R = df['recovered'].to_numpy()\n",
    "M = df['deaths'].to_numpy()\n",
    "S = N - R - I\n",
    "\n",
    "# Creating the time vector\n",
    "t = np.linspace(0, len(I)-1, len(I))\n",
    "\n",
    "until = -1\n",
    "\n",
    "Sd, Id, Md, Rd, td = S[:until], I[:until], M[:until], R[:until], t[:until]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CABRAL - Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import scipy.signal as scs\n",
    "from scipy.interpolate import Rbf, UnivariateSpline\n",
    "\n",
    "\n",
    "spl_I = UnivariateSpline(t, I)\n",
    "spl_R = UnivariateSpline(t, R)\n",
    "spl_S = UnivariateSpline(t, S)\n",
    "\n",
    "time = np.linspace(0, len(I), 10*len(I))\n",
    "\n",
    "I = spl_I(time)\n",
    "R = spl_R(time)\n",
    "S = spl_S(time)\n",
    "\n",
    "until = 40*len(I)\n",
    "dt = 1 / 40\n",
    "\n",
    "Sd, Id, Md, Rd, td = S[:until], I[:until], M[:until], R[:until], t[:until]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(Sd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "size = len(Id) - 1\n",
    "Y = np.zeros((3*size, 1))       # Y   = [ 3*n x 1 ]\n",
    "Phi = np.zeros((3*size, 3))     # Phi = [ 3*n x 2 ]\n",
    " \n",
    "k = 0        # S = [1,2,3,4,5]  e.g. s(k-1) = [1,2,3,4] and s(k) = [2,3,4,5]  \n",
    "for s, i, r, s_, i_, r_ in zip(Sd[:-1], Id[:-1], Rd[:-1], Sd[1:], Id[1:], Rd[1:]):\n",
    "    \n",
    "    # s <- s(k-1)    i <- i(k-1)   r <- r(k-1)\n",
    "    # s_ <- s(k)     i_ <- i(k)    r_ <- r(k) \n",
    "    y_k = np.zeros((3,1))\n",
    "    y_k[0] = s_            # s(k)   ->  b * (- s(k-1) i(k-1) / N) + s(k-1)\n",
    "    y_k[1] = i_            # i(k)   ->  b * (s(k-1) i(k-1) / N) - r i(k-1) + i(k-1)\n",
    "    y_k[2] = r_            # r(k)   ->  r i(k-1) + r(k-1) \n",
    "    \n",
    "    Y[k:k+3,:] = y_k[:,:]     # Y = [ y(0), ..., y(N) ]ˆ\\top\n",
    "    \n",
    "    phi_k = np.zeros((3,3))\n",
    "    phi_k[0,0] = - dt * s * i / N  # - s(k-1) i(k-1) / N   | - s(k-1) i(k-1) / N,     0    , s(k-1)|\n",
    "    phi_k[1,0] = dt * s * i / N    # s(k-1) i(k-1) / N     |   s(k-1) i(k-1) / N,  -i(k-1) , i(k-1)|\n",
    "    phi_k[1,1] = - dt * i          # - i(k-1)              |          0         ,  i(k-1)  , r(k-1)|\n",
    "    phi_k[2,1] = dt * i            # i(k-1)\n",
    "    phi_k[0,2] = s\n",
    "    phi_k[1,2] = i\n",
    "    phi_k[2,2] = r\n",
    "                              #                         | b |\n",
    "    Phi[k:k+3,:] = phi_k[:,:] #                         | r |\n",
    "                              #                         | 1 |\n",
    "    k += 3\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "size = len(I) - 1\n",
    "Phi_sim = np.zeros((3*size, 3))     # Phi = [ 3*n x 2 ]\n",
    " \n",
    "k = 0        # S = [1,2,3,4,5]  e.g. s(k-1) = [1,2,3,4] and s(k) = [2,3,4,5]\n",
    "for s, i, r, s_, i_, r_ in zip(S[:-1], I[:-1], R[:-1], S[1:], I[1:], R[1:]):\n",
    "    \n",
    "    phi_k = np.zeros((3,3))\n",
    "    phi_k[0,0] = - dt * s * i / N  # - s(k-1) i(k-1) / N     | - s(k-1) i(k-1) / N,     0    |\n",
    "    phi_k[1,0] = dt * s * i / N    # s(k-1) i(k-1) / N       |   s(k-1) i(k-1) / N,  -i(k-1) |\n",
    "    phi_k[1,1] = - dt * i          # - i(k-1)                |          0         ,  i(k-1)  |\n",
    "    phi_k[2,1] = dt * i             # i(k-1)\n",
    "    phi_k[0,2] = s\n",
    "    phi_k[1,2] = i\n",
    "    phi_k[2,2] = r\n",
    "                                  #                         | b |\n",
    "    Phi_sim[k:k+3,:] = phi_k[:,:] #                         | r |\n",
    "    \n",
    "    k += 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Y <- Phi * theta\n",
    "# pi Phi Y <- theta\n",
    "\n",
    "theta = np.linalg.pinv(Phi) @ Y\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta[0,0] / theta[1,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "Y_pred = Phi_sim @ theta\n",
    "\n",
    "ds_pred = [ ]\n",
    "di_pred = [ ]\n",
    "dr_pred = [ ]\n",
    "\n",
    "ds_sim = [ ]\n",
    "di_sim = [ ]\n",
    "dr_sim = [ ]\n",
    "\n",
    "phi_k = np.zeros((3,3))\n",
    "phi_k[:,:] = Phi_sim[0:3,:]\n",
    "for k in range(len(I)-1): \n",
    "    # k = [   0  ,   0  ,   0  ,   1  ,   1  ,   1  ,  ]\n",
    "    # i = [   0  ,   1  ,   2  ,   3  ,   4  ,   5  ,  ]\n",
    "    # Y = [ ds(0), di(0), dr(0), ds(1), di(1), dr(1),  ]\n",
    "    ds_pred.append(Y_pred[3*k+0])\n",
    "    di_pred.append(Y_pred[3*k+1])\n",
    "    dr_pred.append(Y_pred[3*k+2])\n",
    "    \n",
    "    # Focus simulation\n",
    "    y_k = phi_k @ theta\n",
    "\n",
    "    ds_sim.append(y_k[0])\n",
    "    di_sim.append(y_k[1])\n",
    "    dr_sim.append(y_k[2])\n",
    "    \n",
    "    ds = [-dt*y_k[0]*y_k[1]/N,         0.0, y_k[0]]\n",
    "    di = [ dt*y_k[0]*y_k[1]/N , -dt*y_k[1], y_k[1]]\n",
    "    dr = [0.0,                   dt*y_k[1], y_k[2]]\n",
    "    \n",
    "    for i in range(3):\n",
    "        phi_k[0,i] = ds[i]\n",
    "        phi_k[1,i] = di[i]\n",
    "        phi_k[2,i] = dr[i]\n",
    "\n",
    "di_pred_2 = N - np.array(ds_pred) - np.array(dr_pred)\n",
    "di_pred_2 = np.reshape(np.array(di_pred_2), (len(I)-1,))\n",
    "ds_pred = np.reshape(np.array(ds_pred), (len(I)-1,))\n",
    "di_pred = np.reshape(np.array(di_pred), (len(I)-1,))\n",
    "dr_pred = np.reshape(np.array(dr_pred), (len(I)-1,))\n",
    "\n",
    "ds_sim = np.reshape(np.array(ds_sim), (len(I)-1,))\n",
    "di_sim = np.reshape(np.array(di_sim), (len(I)-1,))\n",
    "dr_sim = np.reshape(np.array(dr_sim), (len(I)-1,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "p1 = figure(y_range=(-20000,80000))\n",
    "\n",
    "time = np.linspace(0, len(I), len(I))\n",
    "\n",
    "p1.line(time, I, line_dash=\"dashed\", line_width = 3)\n",
    "p1.line(time, R, line_dash=\"dashed\", line_width = 3)\n",
    "\n",
    "p1.line(time, di_pred, color=\"red\", line_width = 2)\n",
    "#p1.line(time, di_pred_2, color=\"green\")\n",
    "p1.line(time, dr_pred, color=\"red\", line_width = 2)\n",
    "\n",
    "p1.scatter(time, di_sim, color=\"green\")\n",
    "\n",
    "show(p1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "p = figure(plot_height=500,\n",
    "           plot_width=600, \n",
    "           toolbar_location=None)\n",
    "\n",
    "time = range(len(I))\n",
    "\n",
    "p.line(time, S, line_width = 2)\n",
    "p.line(time, ds_pred, color=\"green\", line_width = 3)\n",
    "p.line(time, ds_sim, color=\"green\", line_dash=\"dotted\", line_width = 3)\n",
    "\n",
    "show(p)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "# Creating the time vector\n",
    "\n",
    "Ro, beta_par, r_par = [], [], []\n",
    "for until in range(5, len(I), 1):\n",
    "    \n",
    "    I = df['active'].to_numpy()\n",
    "    R = df['recovered'].to_numpy()\n",
    "    M = df['deaths'].to_numpy()\n",
    "    S = N - R - I\n",
    "    \n",
    "    t = np.linspace(0, len(I), len(I))\n",
    "\n",
    "    Sd, Id, Md, Rd, td = S[:until], I[:until], M[:until], R[:until], t[:until]\n",
    "\n",
    "\n",
    "    size = len(Id) - 1\n",
    "    Y = np.zeros((3*size, 1))       # Y   = [ 3*n x 1 ]\n",
    "    Phi = np.zeros((3*size, 3))     # Phi = [ 3*n x 2 ]\n",
    "\n",
    "    k = 0        # S = [1,2,3,4,5]  e.g. s(k-1) = [1,2,3,4] and s(k) = [2,3,4,5]  \n",
    "    for s, i, r, s_, i_, r_ in zip(Sd[:-1], Id[:-1], Rd[:-1], Sd[1:], Id[1:], Rd[1:]):\n",
    "\n",
    "        # s <- s(k-1)    i <- i(k-1)   r <- r(k-1)\n",
    "        # s_ <- s(k)     i_ <- i(k)    r_ <- r(k) \n",
    "        y_k = np.zeros((3,1))\n",
    "        y_k[0] = s_            # s(k)   ->  b * (- s(k-1) i(k-1) / N) + s(k-1)\n",
    "        y_k[1] = i_            # i(k)   ->  b * (s(k-1) i(k-1) / N) - r i(k-1) + i(k-1)\n",
    "        y_k[2] = r_            # r(k)   ->  r i(k-1) + r(k-1) \n",
    "\n",
    "        Y[k:k+3,:] = y_k[:,:]     # Y = [ y(0), ..., y(N) ]ˆ\\top\n",
    "\n",
    "        phi_k = np.zeros((3,3))\n",
    "        phi_k[0,0] = - dt * s * i / N  # - s(k-1) i(k-1) / N   | - s(k-1) i(k-1) / N,     0    , s(k-1)|\n",
    "        phi_k[1,0] = dt * s * i / N    # s(k-1) i(k-1) / N     |   s(k-1) i(k-1) / N,  -i(k-1) , i(k-1)|\n",
    "        phi_k[1,1] = - dt * i          # - i(k-1)              |          0         ,  i(k-1)  , r(k-1)|\n",
    "        phi_k[2,1] = dt * i            # i(k-1)\n",
    "        phi_k[0,2] = s\n",
    "        phi_k[1,2] = i\n",
    "        phi_k[2,2] = r\n",
    "                                  #                         | b |\n",
    "        Phi[k:k+3,:] = phi_k[:,:] #                         | r |\n",
    "                                  #                         | 1 |\n",
    "        k += 3\n",
    "        \n",
    "    theta = np.linalg.pinv(Phi) @ Y\n",
    "    beta_par.append(theta[0,0])\n",
    "    r_par.append(theta[1,0])\n",
    "    Ro.append(theta[0,0]/ theta[1,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = figure(plot_height=300,plot_width=600,y_range=[-120, 120])\n",
    "p2 = figure(plot_height=300,plot_width=600)\n",
    "p3 = figure(plot_height=300,plot_width=600)\n",
    "\n",
    "time = range(5, len(Ro), 1)\n",
    "\n",
    "p1.line(time, Ro, color=\"green\", line_width = 3)\n",
    "p2.line(time, beta_par, color=\"red\", line_width = 3)\n",
    "p3.line(time, r_par, color=\"blue\", line_width = 3)\n",
    "\n",
    "show(column(p1,p2,p3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Criando a superfície custo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from models import *\n",
    "\n",
    "S = N - R - I\n",
    "dataset = dict(S=S, I=I, R=R)\n",
    "\n",
    "beta_grid = np.logspace(-3, 2, 1000)\n",
    "r_grid = np.logspace(-3, 2, 1000)\n",
    "\n",
    "s_error_grid = np.zeros((beta_grid.shape[0], r_grid.shape[0]))\n",
    "i_error_grid = np.zeros((beta_grid.shape[0], r_grid.shape[0]))\n",
    "r_error_grid = np.zeros((beta_grid.shape[0], r_grid.shape[0]))\n",
    "\n",
    "for i, beta in enumerate(beta_grid):\n",
    "    for j, r in enumerate(r_grid):\n",
    "        # Build the model with desired parameters\n",
    "        sir_model = ss.SIR(pop=N, focus=[\"S\", \"I\", \"R\"])\n",
    "        sir_model.parameters = [beta, r]\n",
    "        # Predict the model next step\n",
    "        sim_res = sir_model.predict((S[0],I[0],R[0]), t)\n",
    "        # Compute the model error\n",
    "        s_error_grid[i,j] = (1/np.mean(S))*np.mean(abs(sim_res[0] - S))\n",
    "        i_error_grid[i,j] = (1/np.mean(I))*np.mean(abs(sim_res[1] - I))\n",
    "        r_error_grid[i,j] = (1/np.mean(R))*np.mean(abs(sim_res[2] - R))\n",
    "        #print(\"Step i => \", i, \" and j => \", j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "%matplotlib notebook\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import ticker, cm\n",
    "plt.style.use('dark_background')\n",
    "\n",
    "X, Y = np.meshgrid(r_grid, beta_grid)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9,7))\n",
    "cs = ax.contourf(X, Y, np.log10(r_error_grid), levels=150)\n",
    "cbar = fig.colorbar(cs)\n",
    "ax.set_title(\"Recovered Error\")\n",
    "ax.set_xlabel('r values')\n",
    "ax.set_ylabel('Beta values')\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import ticker, cm\n",
    "plt.style.use('dark_background')\n",
    "\n",
    "X, Y = np.meshgrid(r_grid, beta_grid)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9,7))\n",
    "cs = ax.contourf(X, Y, np.log10(i_error_grid), levels=150)\n",
    "cbar = fig.colorbar(cs)\n",
    "ax.set_title(\"Infected Error\")\n",
    "ax.set_xlabel('r values')\n",
    "ax.set_ylabel('Beta values')\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import ticker, cm\n",
    "plt.style.use('dark_background')\n",
    "\n",
    "X, Y = np.meshgrid(r_grid, beta_grid)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9,7))\n",
    "cs = ax.contourf(X, Y, s_error_grid, levels=150)\n",
    "cbar = fig.colorbar(cs)\n",
    "ax.set_title(\"Suceptible Error\")\n",
    "ax.set_xlabel('r values')\n",
    "ax.set_ylabel('Beta values')\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import ticker, cm\n",
    "plt.style.use('dark_background')\n",
    "\n",
    "X, Y = np.meshgrid(r_grid, beta_grid)\n",
    "total_error = i_error_grid + r_error_grid + s_error_grid\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9,7))\n",
    "cs = ax.contourf(X, Y, np.log10(total_error), levels=150)#, locator=ticker.LogLocator())\n",
    "cbar = fig.colorbar(cs)\n",
    "ax.set_title(\"Total Error\")\n",
    "ax.set_xlabel('r values')\n",
    "ax.set_ylabel('Beta values')\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "S = N - R - I\n",
    "\n",
    "dataset = dict(S=S, I=I, R=R)\n",
    "\n",
    "# Create the model\n",
    "sir_model = ss.SIR(pop=N, focus=[\"S\", \"I\", \"R\"], simulation=\"ivp_continuous\")\n",
    "\n",
    "# Adjust the parameters\n",
    "sir_model.fit(dataset, t, \n",
    "              search_pop=True,\n",
    "              pop_sens=[0.00001,0.05],\n",
    "              Ro_sens=[0.8,15], \n",
    "              D_sens=[5,40])\n",
    "\n",
    "# Predict the model\n",
    "sim_res = sir_model.predict((S[0],I[0],R[0]), t)\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(sir_model.parameters[0]/(sir_model.parameters[1]*sir_model.parameters[2]), 1/sir_model.parameters[1]))\n",
    "print(\"Dias contaminados: \", sir_model.parameters[1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_res = sir_model.predict((S[0],I[0],R[0]), t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "p = 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",
    "p.grid.grid_line_alpha = 0\n",
    "p.ygrid.band_fill_color = \"olive\"\n",
    "p.ygrid.band_fill_alpha = 0.1\n",
    "p.yaxis.axis_label = \"Indivíduos\"\n",
    "p.xaxis.axis_label = \"Dias\"\n",
    "\n",
    "p.line(t, I,  \n",
    "       legend_label=\"Infectados\", color=\"#c2185b\", line_width=4)\n",
    "p.line(t, R,  \n",
    "       legend_label=\"Removidos\", color=\"#0288d1\", line_width=4)\n",
    "\n",
    "# Show the results\n",
    "p.line(t, sim_res[1],  \n",
    "       legend_label=\"Infectados - Modelo\", line_dash=\"dashed\", color=\"#5e35b1\", line_width=3)\n",
    "# p.line(t, sim_res[2],  \n",
    "#        legend_label=\"Infectados - Notificados - Modelo\", line_dash=\"dashed\", color=\"#e91e63\", line_width=3)\n",
    "p.line(t, sim_res[2],  \n",
    "    legend_label=\"Removidos - Modelo\", line_dash=\"dashed\", color=\"#f57f17\", line_width=3)\n",
    "p.line(td, sim_res[0],  \n",
    "      legend_label=\"Susceptible - Modelo\", line_dash=\"dashed\", color=\"#43a047\", line_width=3)\n",
    "\n",
    "show(p)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 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",
    "p.grid.grid_line_alpha = 0\n",
    "p.ygrid.band_fill_color = \"olive\"\n",
    "p.ygrid.band_fill_alpha = 0.1\n",
    "p.yaxis.axis_label = \"Indivíduos\"\n",
    "p.xaxis.axis_label = \"Dias\"\n",
    "\n",
    "p.line(t, I,  \n",
    "       legend_label=\"Infectados\", color=\"#c2185b\", line_width=4)\n",
    "p.line(t, R,  \n",
    "       legend_label=\"Removidos\", color=\"#0288d1\", line_width=4)\n",
    "\n",
    "show(p)"
   ]
  },
  {
   "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",
    "                  beta_sens=[100,10000], \n",
    "                  r_sens=[100,100])\n",
    "    \n",
    "    saved_param['r'].append(sir_model.parameters[0])\n",
    "    saved_param['beta'].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": "markdown",
   "metadata": {},
   "source": [
    "# Visualizando as previsões\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_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=500,\n",
    "           plot_width=600, \n",
    "           tools=\"\", \n",
    "           toolbar_location=None,\n",
    "           title=\"Evolução do COVID - \" + 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",
    "# Incluindo as curvas\n",
    "for data in saved_prediction:\n",
    "    p3.line(pred_t, data[1],\n",
    "           legend_label=\"Previsão Infectados\", \n",
    "           line_cap=\"round\", line_width=4, color=\"#42a5f5\", line_alpha = 0.1)\n",
    "\n",
    "    p3.line(pred_t, data[0],\n",
    "       legend_label=\"Previsão Suscetiveis\", \n",
    "       line_cap=\"round\", line_width=4, color=\"#ff5722\", line_alpha = 0.07)\n",
    "    \n",
    "    p3.line(pred_t, data[2],\n",
    "       legend_label=\"Previsão Recuperados\", \n",
    "       line_cap=\"round\", line_width=4, 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\", line_dash = 'dashed')\n",
    "\n",
    "if plot_all:\n",
    "    p3.line(td, Rd,\n",
    "           legend_label=\"Recuperados\", \n",
    "           line_cap=\"round\", line_width=3, color=\"#5e35b1\", line_dash = 'dashed')\n",
    "\n",
    "    p3.line(td, N*saved_param['pop'][-1] - Rd - Id,\n",
    "           legend_label=\"Suscetiveis\", \n",
    "           line_cap=\"round\", line_width=3, color=\"#b71c1c\", line_dash = 'dashed')\n",
    "\n",
    "# Colocando as legendas\n",
    "p3.legend.click_policy=\"hide\"\n",
    "p3.legend.location = \"top_right\"\n",
    "\n",
    "show(p3)\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 = len(Id) - np.argmax(signal)\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1.write_html(\"./docs/media_content/peak_estimate_CN.html\")"
   ]
  }
 ],
 "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
}