<jupyter_start><jupyter_text>Forecasting models<jupyter_code>import sys
import os
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from qgis.core import *
from qgis.PyQt.QtGui import *
from qgis.PyQt.QtCore import *
from IPython.display import Image
QgsApplication.setPrefixPath(r'C:\\OSGeo4W64\\apps\\qgis', True)
qgs = QgsApplication([], True)
qgs.initQgis()
sys.path.append(r'C:\OSGeo4W64\apps\qgis\python\plugins')
project = QgsProject.instance()
path = "C:\\OSGeo4W64\\bin\\SIG\\Projeto_Italy\\Mapa\\ITA_adm1.shp"
map_layer = QgsVectorLayer(path, 'Italy map', 'ogr')

if not map_layer.isValid():
    print("Failed to load the layer!")
else:
    project.addMapLayer(map_layer)
    print("Sucess")
csv_path = "file:///C:/OSGeo4W64/bin/SIG/Projeto_Italy/Dataset/covid_italy.csv?delimiter=,'"
csv_layer = QgsVectorLayer(csv_path, 'Data', 'delimitedtext')

if not csv_layer.isValid():
    print('Layer failed to load!')
else:
    project.addMapLayer(csv_layer)
    print("Sucess")
joinName = 'name_region'
targetName = 'NAME_1'

joinObject = QgsVectorLayerJoinInfo()
joinObject.setJoinFieldName(joinName)
joinObject.setTargetFieldName(targetName)
joinObject.setJoinLayerId(csv_layer.id())
joinObject.setUsingMemoryCache(True)
joinObject.setJoinLayer(csv_layer)
flag = map_layer.addJoin(joinObject)
## For data
import pandas as pd
import numpy as np
## For plotting
import matplotlib.pyplot as plt
## For parametric fitting
from scipy import optimize<jupyter_output><empty_output><jupyter_text>**Lombardia** Most affected region of Italy<jupyter_code>values_case = []
values_deaths = []
values_new = []
datas = []

lowerTotal = 0
i = 0


for feature in csv_layer.getFeatures():
    if feature['name_region']=="Lombardia":
        if feature['total_positive']>0:
            values_case.append(feature['total_case'])
            values_deaths.append(feature['deaths'])
            values_new.append(feature['new_positive'])
            datas.append(feature['data'])


fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True, figsize=(13,7))
y_pos = np.arange(len(datas))

ax[0].scatter(datas, values_case, color="indigo")
plt.xticks(y_pos, datas, rotation='vertical')
ax[0].set(title="Total cases")
ax[1].bar(datas, values_deaths, color="rebeccapurple")
ax[1].set(title="Deaths")
plt.savefig('../Imagens/Prevision_lombardia_totalcase_deaths.png', dpi=300, format='png')
plt.show()<jupyter_output><empty_output><jupyter_text>**Number of cases**<jupyter_code>'''
Linear function: f(x) = a + b*x
'''
def f(x):
    return 10 + 1500*x

y_linear = f(x=np.arange(len(values_case)))

'''
Exponential function: f(x) = a + b^x
'''
def f(x):
    return 10 + 1.18**x

y_exponential = f(x=np.arange(len(values_case)))


'''
Logistic function: f(x) = a / (1 + e^(-b*(x-c)))
'''
def f(x): 
    return 90000 / (1 + np.exp(-0.5*(x-20)))

y_logistic = f(x=np.arange(len(values_case)))
fig, ax = plt.subplots(figsize=(13,5))
y_pos = np.arange(len(datas))


ax.scatter(y_pos, values_case, alpha=0.5, color='black')
plt.xticks(y_pos, datas, rotation='vertical')
ax.plot(y_pos, y_linear, label="linear", color="darkgreen")
ax.plot(y_pos, y_exponential, label="exponential", color="seagreen")
ax.plot(y_pos, y_logistic, label="logistic", color="limegreen")
plt.ylabel('Number of cases')
plt.xlabel('Dates')
plt.title('Number of cases in Lombardia')
plt.savefig('../Imagens/Prevision_lombardia_numberofcases.png', dpi=300, format='png')
plt.show()<jupyter_output><empty_output><jupyter_text>**Number of deaths**<jupyter_code>'''
Linear function: f(x) = a + b*x
'''
def f(x):
    return 10 + 1500*x

y_linear = f(x=np.arange(len(values_deaths)))

'''
Exponential function: f(x) = a + b^x
'''
def f(x):
    return 10 + 1.18**x

y_exponential = f(x=np.arange(len(values_deaths)))


'''
Logistic function: f(x) = a / (1 + e^(-b*(x-c)))
'''
def f(x): 
    return 90000 / (1 + np.exp(-0.5*(x-20)))

y_logistic = f(x=np.arange(len(values_deaths)))
fig, ax = plt.subplots(figsize=(13,5))
y_pos = np.arange(len(datas))


ax.scatter(y_pos, values_case, alpha=0.5, color='black')
plt.xticks(y_pos, datas, rotation='vertical')
ax.plot(y_pos, y_linear, label="linear", color="peru")
ax.plot(y_pos, y_exponential, label="exponential", color="goldenrod")
ax.plot(y_pos, y_logistic, label="logistic", color="khaki")
plt.ylabel('Number of deaths')
plt.xlabel('Dates')
plt.title('Number of deaths in Lombardia')
plt.savefig('../Imagens/Prevision_lombardia_numberofdeaths.png', dpi=300, format='png')
plt.show()<jupyter_output><empty_output><jupyter_text>**Parametric Fitting**<jupyter_code>df=pd.read_csv(open("../Dataset/lombardia_last.csv"))
df.index = df['data']
date = df['data'].tolist()
xvalues=np.arange(len(date))

totalcases = df['total_case'].tolist()

print(date)
print(totalcases)

'''
Logistic function: f(x) = capacity / (1 + e^-k*(x - midpoint) )
'''
def logistic_f(xdata, c, k, m):
    ydata = c / (1 + np.exp(-k*(xdata-m)))
    return ydata
## optimize from scipy
logistic_model, cov = optimize.curve_fit(logistic_f, xdata=xvalues, ydata=totalcases, maxfev=100000, p0=[np.max(totalcases), 1, 1])
## print the parameters
logistic_model
'''
Plot parametric fitting.
'''
def utils_plot_parametric(dtf, zoom=30, figsize=(15,5)):
    ## interval
    dtf["residuals"] = dtf["ts"] - dtf["model"]
    dtf["conf_int_low"] = dtf["forecast"] - 1.96*dtf["residuals"].std()
    dtf["conf_int_up"] = dtf["forecast"] + 1.96*dtf["residuals"].std()
    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=figsize)
    
    ## entire series
    dtf["ts"].plot(marker=".", linestyle='None', ax=ax[0], title="Parametric Fitting", color="black")
    dtf["model"].plot(ax=ax[0], color="maroon")
    dtf["forecast"].plot(ax=ax[0], grid=True, color="salmon")
    ax[0].fill_between(x=dtf.index, y1=dtf['conf_int_low'], y2=dtf['conf_int_up'], color='b', alpha=0.3)
    
    plt.show()
    return dtf[["ts","model","residuals","conf_int_low","forecast","conf_int_up"]]
'''
Forecast unknown future.
:parameter
    :param ts: pandas series
    :param f: function
    :param model: list of optim params
    :param pred_ahead: number of observations to forecast (ex. pred_ahead=30)
    :param freq: None or str - 'B' business day, 'D' daily, 'W' weekly, 'M' monthly, 'A' annual, 'Q' quarterly
    :param zoom: for plotting
'''
def forecast_curve(ts, f, model, pred_ahead=None, zoom=30, figsize=(15,5)):
    ## fit
    X = np.arange(len(ts))
    fitted = f(X, model[0], model[1], model[2])
    dtf = ts.to_frame(name="ts")
    dtf["model"] = fitted
    
    ## index
    start=0
    index = pd.date_range(start=start,periods=pred_ahead)
    index = index[1:]
    ## forecast
    Xnew = np.arange(len(ts)+1, len(ts)+1+len(index))
    preds = f(Xnew, model[0], model[1], model[2])
    dtf = dtf.append(pd.DataFrame(data=preds, index=index, columns=["forecast"]))
    
    ## plot
    utils_plot_parametric(dtf, zoom=zoom)
    return dtf
preds = forecast_curve(df["total_case"], logistic_f, logistic_model, pred_ahead=30, zoom=7)<jupyter_output><empty_output><jupyter_text>Molise **Less affected region of Italy**<jupyter_code>values_case = []
values_deaths = []
values_new = []
datas = []

lowerTotal = 0
i = 0


for feature in csv_layer.getFeatures():
    if feature['name_region']=="Molise":
        if feature['total_positive']>0:
            values_case.append(feature['total_case'])
            values_deaths.append(feature['deaths'])
            values_new.append(feature['new_positive'])
            datas.append(feature['data'])
fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True, figsize=(13,7))
y_pos = np.arange(len(datas))

ax[0].scatter(datas, values_case, color="darkcyan")
plt.xticks(y_pos, datas, rotation='vertical')
ax[0].set(title="Total cases")
ax[1].bar(datas, values_deaths, color="c")
ax[1].set(title="Deaths")
plt.savefig('../Imagens/Prevision__Molise_totalcase_deaths.png', dpi=300, format='png')
plt.show()<jupyter_output><empty_output><jupyter_text>**Number of cases**<jupyter_code>'''
Linear function: f(x) = a + b*x
'''
def f(x):
    return 10 + 1500*x

y_linear = f(x=np.arange(len(values_case)))

'''
Exponential function: f(x) = a + b^x
'''
def f(x):
    return 10 + 1.18**x

y_exponential = f(x=np.arange(len(values_case)))


'''
Logistic function: f(x) = a / (1 + e^(-b*(x-c)))
'''
def f(x): 
    return 90000 / (1 + np.exp(-0.5*(x-20)))

y_logistic = f(x=np.arange(len(values_case)))
fig, ax = plt.subplots(figsize=(13,5))
y_pos = np.arange(len(datas))


ax.scatter(y_pos, values_case, alpha=0.5, color='black')
plt.xticks(y_pos, datas, rotation='vertical')
ax.plot(y_pos, y_linear, label="linear", color="red")
ax.plot(y_pos, y_exponential, label="exponential", color="crimson")
ax.plot(y_pos, y_logistic, label="logistic", color="firebrick")
plt.ylabel('Number of cases')
plt.xlabel('Dates')
plt.title('Number of cases in Molise')
plt.savefig('../Imagens/Prevision_molise_numberofcases.png', dpi=300, format='png')
plt.show()<jupyter_output><empty_output><jupyter_text>**Number of deaths**<jupyter_code>'''
Linear function: f(x) = a + b*x
'''
def f(x):
    return 10 + 1500*x

y_linear = f(x=np.arange(len(values_deaths)))

'''
Exponential function: f(x) = a + b^x
'''
def f(x):
    return 10 + 1.18**x

y_exponential = f(x=np.arange(len(values_deaths)))


'''
Logistic function: f(x) = a / (1 + e^(-b*(x-c)))
'''
def f(x): 
    return 90000 / (1 + np.exp(-0.5*(x-20)))

y_logistic = f(x=np.arange(len(values_deaths)))
fig, ax = plt.subplots(figsize=(13,5))
y_pos = np.arange(len(datas))


ax.scatter(y_pos, values_case, alpha=0.5, color='black')
plt.xticks(y_pos, datas, rotation='vertical')
ax.plot(y_pos, y_linear, label="linear", color="navy")
ax.plot(y_pos, y_exponential, label="exponential", color="royalblue")
ax.plot(y_pos, y_logistic, label="logistic", color="blue")
plt.ylabel('Number of deaths')
plt.xlabel('Dates')
plt.title('Number of deaths in Molise')
plt.savefig('../Imagens/Prevision_molise_numberofdeaths.png', dpi=300, format='png')
plt.show()
df=pd.read_csv(open("../Dataset/molise_last.csv"))
df.index = df['data']
date = df['data'].tolist()
xvalues=np.arange(len(date))

totalcases = df['total_case'].tolist()

print(date)
print(totalcases)

'''
Logistic function: f(x) = capacity / (1 + e^-k*(x - midpoint) )
'''
def logistic_f(xdata, c, k, m):
    ydata = c / (1 + np.exp(-k*(xdata-m)))
    return ydata
## optimize from scipy
logistic_model, cov = optimize.curve_fit(logistic_f, xdata=xvalues, ydata=totalcases, maxfev=100000, p0=[np.max(totalcases), 1, 1])
## print the parameters
logistic_model
'''
Plot parametric fitting.
'''
def utils_plot_parametric(dtf, zoom=30, figsize=(15,5)):
    ## interval
    dtf["residuals"] = dtf["ts"] - dtf["model"]
    dtf["conf_int_low"] = dtf["forecast"] - 1.96*dtf["residuals"].std()
    dtf["conf_int_up"] = dtf["forecast"] + 1.96*dtf["residuals"].std()
    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=figsize)
    
    ## entire series
    dtf["ts"].plot(marker=".", linestyle='None', ax=ax[0], title="Parametric Fitting", color="black")
    dtf["model"].plot(ax=ax[0], color="dodgerblue")
    dtf["forecast"].plot(ax=ax[0], grid=True, color="fuchsia")
    ax[0].fill_between(x=dtf.index, y1=dtf['conf_int_low'], y2=dtf['conf_int_up'], color='b', alpha=0.3)
   
    plt.show()
    return dtf[["ts","model","residuals","conf_int_low","forecast","conf_int_up"]]
'''
Forecast unknown future.
:parameter
    :param ts: pandas series
    :param f: function
    :param model: list of optim params
    :param pred_ahead: number of observations to forecast (ex. pred_ahead=30)
    :param freq: None or str - 'B' business day, 'D' daily, 'W' weekly, 'M' monthly, 'A' annual, 'Q' quarterly
    :param zoom: for plotting
'''
def forecast_curve(ts, f, model, pred_ahead=None, zoom=30, figsize=(15,5)):
    ## fit
    X = np.arange(len(ts))
    fitted = f(X, model[0], model[1], model[2])
    dtf = ts.to_frame(name="ts")
    dtf["model"] = fitted
    
    ## index
    start=0
    index = pd.date_range(start=start,periods=pred_ahead)
    index = index[1:]
    ## forecast
    Xnew = np.arange(len(ts)+1, len(ts)+1+len(index))
    preds = f(Xnew, model[0], model[1], model[2])
    dtf = dtf.append(pd.DataFrame(data=preds, index=index, columns=["forecast"]))
    
    ## plot
    utils_plot_parametric(dtf, zoom=zoom)
    return dtf
preds = forecast_curve(df["total_case"], logistic_f, logistic_model, pred_ahead=30, zoom=7)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18567168990249014, 'fraction_numerical': 0.021597227818518815, 'mean_word_length': 5.296296296296297, 'pattern_counts': {'":': 2, '<': 33, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '514092', 'license': 'CC0-1.0', 'n_tokens_mistral': 5099, 'n_tokens_neox': 4836, 'n_words': 881, 'path': 'Casos de estudo/Italy_Covid19/Projeto_Italy/Notebooks/Italy_prevision.ipynb', 'repo_name': 'kika-nogueira97/covid-pt', 'size': '12409'}