Come creare un grafico previsionale sulla diffusione del Covid-19 con Linux, Bash e Python

All’inizio dell’emergenza sanitaria da Coronavirus (nCoV-19), durante i primi giorni di lockdown, mi sono dato subito da fare mettendo in piedi una semplice webapp www.coronavirus-italy.it, per il monitoraggio, tramite API nazionali e internazionali, della pandemia.

Poi arricchendola di alcune funzionalità più specifiche, una di queste: le previsioni sulla diffusione in Italia del virus, in base ai dati forniti giorno per giorno dalla Protezione Civile.

In questo articolo voglio condividere la mia esperienza per mettere in piedi un sistema previsionale automatico sul nostro server Linux.

Iniziamo con la preparazione dell’ambiente.

Nel nostro server Linux dobbiamo accertarci di avere installato Python3 e questi moduli necessari: scipy, venv, pandas, sklearn, numpy, matplotlib.

Ora possiamo iniziare con la creazione del file Python che genererà i grafici previsionali con i dati della Protezione Civile.

Teniamo presente che i nuovi dati vengono diffusi ogni giorno intorno alle 18.00 circa, e comunque entro le 18.30.

Per fare i calcoli previsionali tramite modello matematico logistico, mi sono affidato a questo progetto: https://github.com/marcello-dev/coronavirus-forecast, facendone le opportune modifiche che condivido qui.

Possiamo inserire questo codice all’interno di un nuovo file di testo che chiameremo covid-forecast.py

import pandas as pd
import numpy as np
from datetime import datetime, timedelta
from sklearn.metrics import mean_squared_error
from scipy.optimize import curve_fit
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

url = "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale/dpc-covid19-ita-andamento-nazionale.csv"
df = pd.read_csv(url)

df = df.loc[:, ['data', 'totale_casi']]

FMT = '%Y-%m-%dT%H:%M:%S'
date = df['data']
df['data'] = date.map(lambda x2: (datetime.strptime(x2, FMT) - datetime.strptime("2020-01-01T00:00:00", FMT)).days)

def logistic_model(x3, a3, b3, c3):
    return c3/(1+np.exp(-(x3-b3)/a3))

x = list(df.iloc[:, 0])
y = list(df.iloc[:, 1])

fit = curve_fit(logistic_model, x, y, p0=[2, 100, 20000])

a = fit[0][0]
b = fit[0][1]
c = fit[0][2]

first_january = datetime.strptime('2020/01/01', "%Y/%m/%d")
a = fit[0][0]
b = fit[0][1]
c = fit[0][2]

first_january = datetime.strptime('2020/01/01', "%Y/%m/%d")
infection_peak_date = first_january + timedelta(days=int(b))
print('{"status": "ok","previsioni": [{ "piccoContagio": "',datetime.strftime(infection_peak_date,"%d/%m/%Y"),'",')
errors = [np.sqrt(fit[1][i][i]) for i in [0, 1, 2]]
print('"totaleInfetti": "{}", "minInfetti": "{}", "maxInfetti": "{}",'.format(int(c), int(c-errors[2]),int(c+errors[2])))

sol = int(fsolve(lambda x : logistic_model(x, a, b, c) - int(c),b))

infection_end_date = first_january + timedelta(days=int(sol))
print('"fineContagio": "',datetime.strftime(infection_end_date,"%d/%m/%Y"),'",')

def add_real_data(df, label, color=None):
    x = df['data'].tolist()
    y = df['totale_casi'].tolist()
    plt.scatter(x, y, label="Dati reali (" + label + ")", c=color)

pred_x = list(range(max(x), sol))
plt.rcParams['figure.figsize'] = [7, 7]
plt.rc('font', size=14)
# Real data
add_real_data(df[-1:], "oggi")
add_real_data(df[-2:-1], "ieri")
add_real_data(df[:-2], "2 giorni fa")
# Predicted curve of today
plt.plot(x+pred_x, [logistic_model(i,fit[0][0],fit[0][1],fit[0][2]) for i in x+pred_x], label="Previsione dati oggi")

# Predicted curve of yesterday
x = list(df[:-1].iloc[:, 0])
y = list(df[:-1].iloc[:, 1])
pred_x = list(range(max(x), sol))
fit = curve_fit(logistic_model, x, y, p0=[2, 100, 20000])
plt.plot(x+pred_x, [logistic_model(i,fit[0][0],fit[0][1],fit[0][2]) for i in x+pred_x],
         label="Previsione dati ieri", dashes=[4, 4])

# Predicted curve of 2 days ago curve
x = list(df[:-2].iloc[:, 0])
y = list(df[:-2].iloc[:, 1])
pred_x = list(range(max(x), sol))
fit = curve_fit(logistic_model, x, y, p0=[2, 100, 20000])
plt.plot(x+pred_x, [logistic_model(i,fit[0][0],fit[0][1],fit[0][2]) for i in x+pred_x],
         label="Previsione dati 2 giorni fa",dashes=[8, 8])

today_date = datetime.today().strftime('%Y-%m-%d')
plt.title("Previsioni per casi accertati in Italia del " + today_date)

plt.legend()
plt.xlabel("Giorni da 1 gennaio 2020")
plt.ylabel("Numero totale persone infette")
plt.ylim((min(y)*0.9,c*1.1))

filename = 'plot-' + today_date + '.png'
plt.savefig('plots/'+filename, bbox_inches="tight")
y_pred_logistic = [logistic_model(i,fit[0][0],fit[0][1],fit[0][2]) for i in x]
print('"erroreModello": "', mean_squared_error(y,y_pred_logistic),'"}]}')

Creiamo allo stesso livello del nostro file Python una cartella che chiameremo plots/.

A questo punto, prima di procedere facciamo il test, per verificare che tutti gli import e le istruzioni Python vadano a buon fine, digitiamo questi 3 comandi uno dopo l’altro, il tutto si dovrebbe eseguire senza errori.

python3 -m venv venv

source venv/bin/activate

python3 main.py > plots/previsione.json

Tutta la parte di codice di generazione del file Json assieme al grafico previsionale è una modifica che ho implementato io e non troverete nel progetto originale, ad ogni modo la condivido in quanto ritenuta utile a future implementazioni di ogni genere, rimane pur sempre un JSON aggiornato quotidianamente contenente dei dati.

Possiamo ora procedere all’automatizzazione del tutto semplicemente creando uno script Bash che chiameremo covidscript:

#!/bin/bash
cd /home/Covid/coronavirus-forecast
python3 -m venv venv
source venv/bin/activate
python3 main.py > plots/previsione.json

Ovviamente ricordate di personalizzare l'indirizzo della cartella dove risiede il vostro applicativo Python per questo progetto, dichiarato nella seconda riga di questo script appena creato.

Se vogliamo inserire lo script in un Cron automatico ricordiamoci di non inserire alcuna estensione al file Bash e di renderlo eseguibile:

chmod +x covidscript

Poi effettuiamo un link simbolico dello script alla cartella di Cron:

ln -s /home/covid/covidscript /etc/cron.daily

A questo punto ogni giorno avremo sempre aggiornati il grafico previsionale e il file Json pronto per qualsiasi utilizzo.

Se vogliamo modificare l’orario dei Cron giornalieri possiamo inserire una riga all’interno del file /etc/crontab ricordandoci di inserire un orario non inferiore alle 18.30, perché altrimenti rischiamo ogni giorno di essere troppo in anticipo rispetto alla diffusione pubblica dei dati di quel giorno:

35 18   * * *   root    test -x /usr/sbin/anacron || ( cd / && run-parts --report /etc/cron.daily )

Currently there are no comments, so be the first!