# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ### Carga das bibliotecas

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split

pd.set_option('display.max_columns', 100) # ajusta a exibicao para um numero de colunas superior a quantidade padrao

# ### Carga dos dados de treino e teste

dataset = pd.read_csv('../../data/FINAL/lista_dados_sage_FINAL.csv', sep=';')

dataset

pd.plotting.lag_plot(dataset.loc[dataset.cidade.str.contains('Salvador'), 'total_confirmado_abs'])
plt.show()


# +
def create_lag_columns(df, n_in=1, n_out=1, dropnan=True):
    """
    Frame a time series as a supervised learning dataset.
    Arguments:
        df: Input Pandas DataFrame.
        n_in: Number of lag observations as input (X).
        n_out: Number of observations as output (y).
        dropnan: Boolean whether or not to drop rows with NaN values.
    Returns:
        Pandas DataFrame of series framed for supervised learning.
    """
    n_vars = df.shape[1]
    cols, names = list(), list()
    # input sequence (t-n, ... t-1)
    for i in range(n_in, 0, -1):
        cols.append(df.shift(i))
        names += [('%s_ant%d' % (df.columns.to_list()[j], i)) for j in range(n_vars)]
    """
    # forecast sequence (t, t+1, ... t+n)
    for i in range(0, n_out):
        cols.append(df.shift(-i))
        if i == 0:
            names += df.columns.to_list()
        else:
            names += [('%s_pos%d' % (df.columns.to_list()[j], i)) for j in range(n_vars)]
    """
    # put it all together
    agg = pd.concat(cols, axis=1)
    agg.columns = names
    # drop rows with NaN values
    if dropnan:
        agg.dropna(inplace=True)
    return agg

lag_cols = ['evaporacao_piche', 'insolacao', 'precipitacao', 'temp_comp_media', 'temp_maxima', 'temp_minima',
            'umidade_relativa_media', 'velocidade_do_vento_media', 'google_trends_score', 'total_confirmado_abs']
df_lag_cols = create_lag_columns(dataset[lag_cols], 11, 1, False)
#df_lag_cols.corr().loc[df_lag_cols.corr()['evaporacao_piche_ant11'] > 0.7]
df_lag_cols
# -

df_parcial = pd.concat([dataset[['ibge', 'uf', 'cidade', 'n_ano', 'n_mes',]],
                        df_lag_cols,
                        dataset[['t_agua', 't_banagua', 't_dens', 't_lixo', 't_luz', 'agua_esgoto', 'evaporacao_piche',
                                 'insolacao', 'precipitacao', 'temp_comp_media', 'temp_maxima', 'temp_minima',
                                 'umidade_relativa_media', 'velocidade_do_vento_media', 'google_trends_score',
                                 'total_confirmado_abs']]], axis=1)

# +
# arquivo parcial gerado para o vinicius trabalhar uma versao de modelo em cima dele
# sera identico a nova versao final do arquivo lista_dados_sage_FINAL.csv, que sera ajustado no notebook Juncao_Final_Dados_SAGE

# csv
df_parcial.to_csv('../../data/FINAL/lista_dados_sage_PARCIAL.csv', index=False, sep=';')
# csv
df1 = pd.read_csv('../../data/FINAL/lista_dados_sage_PARCIAL.csv', sep=';')
df1
# -

dataset.info()

# ### Análise Descritiva

dataset.shape

dataset.describe(include='all')

dataset

# ### Feature Engineering

# Serão criadas as seguintes colunas:
# * **`total_confirmado_abs_1m`**: Equivale ao conteúdo da coluna **`total_confirmado_abs`** do mês anterior dentro do respectivo município;
# * **`total_confirmado_acum_2m`**: Acumulado dos 2 meses anteriores (exclui o atual) para a coluna **`total_confirmado_abs`** dentro do respectivo município;
# * **`total_confirmado_acum_3m`**: Acumulado dos 3 meses anteriores (exclui o atual) para a coluna **`total_confirmado_abs`** dentro do respectivo município;
# * **`estacao_ano`**: Estação do ano (verao = 1; outono = 2; inverno = 3; primavera = 4);
# * **`n_trimestre`**: Trimestre do mês;
#
# Para os totais acumulados, será criado um Dataframe auxiliar (**`df_ibge_ano_mes`**) contendo todos os 12 meses para cada municipio, de forma a calcular corretamente os meses acumulados.

# +
# este dataframe serve para, ao juntar com a base original, preencher as lacunas de meses inexistentes na base original.
# isso e necessario para gerar os acumulados dos ultimos 2, 3 e 4 meses

# lista com os anos distintos da base original
n_ano_list = dataset.n_ano.unique().tolist()

# lista fixa com os numeros dos meses de um ano (de 1 a 12)
n_mes_list = list(range(1, 13))

# lista com os codigos ibge distintos da base original
ibge_list = dataset.ibge.unique().tolist()

# gera uma lista aninhada em 3 niveis, combinando cada elemento das 3 listas anteriores para formar um dicionario
# o objetivo aqui e gerar uma linha para cada combinacao distinta de {n_ano, n_mes, ibge}
list_of_list = [[[{'n_ano': i, 'n_mes': j, 'ibge': k} for i in n_ano_list] for j in n_mes_list] for k in ibge_list]

# transforma a lista aninhada em uma lista simples de dicionarios
list_of_dict = [item for l3 in list_of_list for l2 in l3 for item in l2]

# gera o dataframe a partir da lista de dicionarios
df_ibge_ano_mes = pd.DataFrame(list_of_dict)
df_ibge_ano_mes

# +
# junta (merge) os dataframes df_ibge_ano_mes e dataset, gerando linhas com NaN para os meses em que
# a base original nao possui registros
dataset2 = pd.merge(df_ibge_ano_mes, dataset, on = ['ibge', 'n_ano', 'n_mes'], how = 'left')

# ordena a base para garantir que os totais acumulados serao correspondentes aos meses adjacentes
dataset2.sort_values(by = ['ibge', 'n_ano', 'n_mes'], inplace = True)

# trataremos aqui apenas a coluna total_confirmado_abs, pois essas linhas geradas serao removidas apos
# o calculo dos dados acumulados
dataset2.loc[dataset2.total_confirmado_abs.isnull(), 'total_confirmado_abs'] = 0
dataset2.total_confirmado_abs = dataset2.total_confirmado_abs.astype(int)

# cria as colunas total_confirmado_abs_1m, total_confirmado_acum_2m e total_confirmado_acum_3m
df_2m = dataset2[['total_confirmado_abs']].shift(periods=1, fill_value=0).rename(columns = {'total_confirmado_abs': 'total_confirmado_abs_1m'})
df_3m = dataset2[['total_confirmado_abs']].shift(periods=2, fill_value=0).rename(columns = {'total_confirmado_abs': 'total_confirmado_acum_2m'})
df_4m = dataset2[['total_confirmado_abs']].shift(periods=3, fill_value=0).rename(columns = {'total_confirmado_abs': 'total_confirmado_acum_3m'})

# junta as colunas criadas ao dataframe original
dataset2 = pd.concat([dataset2, df_2m, df_3m, df_4m], axis = 1)

# zera os acumulados que envolvem dados de outros municipios, deslocados pela funcao shift
dataset2.loc[dataset2.ibge != dataset2.ibge.shift(periods=1), 'total_confirmado_abs_1m'] = 0
dataset2.loc[dataset2.ibge != dataset2.ibge.shift(periods=2), 'total_confirmado_acum_2m'] = 0
dataset2.loc[dataset2.ibge != dataset2.ibge.shift(periods=3), 'total_confirmado_acum_3m'] = 0

# por fim, soma os dados absolutos para gerar os acumulados dos ultimos 3 e 2 meses (nesta ordem, para nao distorcer os acumulados)
dataset2.total_confirmado_acum_3m = dataset2.total_confirmado_abs_1m + dataset2.total_confirmado_acum_2m + dataset2.total_confirmado_acum_3m
dataset2.total_confirmado_acum_2m = dataset2.total_confirmado_abs_1m + dataset2.total_confirmado_acum_2m

# descarta as linhas inseridas apenas para o calculo dos acumulados
dataset2 = dataset2.loc[dataset2.uf.notnull()]
dataset2
# -

# formato: {mes: estacao}, onde as estacoes do ano sao:
# verao = 1
# outono = 2
# inverno = 3
# primavera = 4
month_to_season = {1: 1, 2: 1, 3: 2, 4: 2, 5: 2, 6: 3, 7: 3, 8: 3, 9: 4, 10: 4, 11: 4, 12: 1}
month_to_season

# cria as colunas estacao_ano e n_trimestre e reorganiza as colunas
dataset2['estacao_ano'] = dataset2.n_mes.map(month_to_season)
dataset2['n_trimestre'] = np.ceil(dataset2.n_mes/3).astype(int)
dataset2 = dataset2[['ibge', 'n_ano', 'n_mes', 'uf', 'cidade',
       't_agua', 't_banagua', 't_dens', 't_lixo', 't_luz', 'agua_esgoto',
       'evaporacao_piche', 'insolacao', 'precipitacao', 'temp_comp_media',
       'temp_maxima', 'temp_minima', 'umidade_relativa_media', 'velocidade_do_vento_media',
       'google_trends_score', 'estacao_ano', 'n_trimestre',
       'total_confirmado_abs_1m', 'total_confirmado_acum_2m',
       'total_confirmado_acum_3m', 'total_confirmado_abs']]
dataset2

dataset2.info()

# ### Tratamento das variáveis categóricas

# Os campos abaixo, apesar de serem do tipo numérico, são na verdade variáveis categóricas, uma vez que seu conteúdo não tem valor matemático, e sim, de identificação:
# * **`ibge`**;
# * **`n_ano`**;
# * **`n_mes`**;
# * **`estacao_ano`**;
# * **`n_trimestre`**;
#
# Entretanto, os campos abaixo serão descartados do modelo pelos seguintes motivos:
# * **`ibge`**: Representa cada município e seu respectivo estado (que já possuem os campos **`uf`** e **`cidade`** para representá-los diretamente);
# * **`n_ano`**: Possui um único valor, uma vez que a base possui dados apenas do ano de 2016;
#
# Dessa forma, apenas as variáveis abaixo serão tratadas:
# * **`uf`**;
# * **`cidade`**;
# * **`n_mes`**;
# * **`estacao_ano`**;
# * **`n_trimestre`**;

dataset2.columns

dataset3 = dataset2.drop(columns = ['ibge', 'n_ano'])

# +
'''
X = dataset3.drop(columns = ['total_confirmado_abs'])
y = dataset3.total_confirmado_abs
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=11)
'''

X = dataset.drop(columns = ['ibge', 'n_ano', 'total_confirmado_abs'])
y = dataset.total_confirmado_abs
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=11)
# -

dataset.info()

# +
from category_encoders import OneHotEncoder, BinaryEncoder

#encoder = OrdinalEncoder(cols = ['n_mes', 'uf', 'cidade', 'estacao_ano', 'n_trimestre'])
encoder = OneHotEncoder(cols = ['n_mes', 'uf'])

encoder.fit(X_train, y_train)
X_train = encoder.transform(X_train)
X_test = encoder.transform(X_test)

encoder = BinaryEncoder(cols = 'cidade')

encoder.fit(X_train, y_train)
X_train = encoder.transform(X_train)
X_test = encoder.transform(X_test)
# -

X_train.shape

X_test.shape

# ## Random Forest

# +
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score

reg = RandomForestRegressor()
reg.fit(X_train, y_train)

mse = mean_squared_error(y_test, reg.predict(X_test))
print("MSE: %.4f" % mse)
r2 = r2_score(y_test, reg.predict(X_test))
print("R2: %.4f" % r2)
# -

X_train.head()









# +
from sklearn.model_selection import GridSearchCV

tuned_parameters = [{'max_depth': [5, 10, 15],
                     'max_features': [1, 5, 10],
                     'n_estimators': [10, 30, 50]}]

clf = GridSearchCV(RandomForestRegressor(random_state=11),
                   tuned_parameters, cv=5, scoring='r2')
clf.fit(X_train, y_train)

# +
print("Best parameters set found on development set:")
print()
print(clf.best_params_)
print()
print("Grid scores on development set:")
print()
means = clf.cv_results_['mean_test_score']
stds = clf.cv_results_['std_test_score']
for mean, std, params in zip(means, stds, clf.cv_results_['params']):
    print("%0.3f (+/-%0.03f) for %r"
        % (mean, std * 2, params))
print()

print("Detailed classification report:")
print()
print("The model is trained on the full development set.")
print("The scores are computed on the full evaluation set.")
print()
mse = mean_squared_error(y_test, clf.predict(X_test))
print("MSE: %.4f" % mse)
r2 = r2_score(y_test, clf.predict(X_test))
print("R2: %.4f" % r2)

# +
tuned_parameters = {'max_depth': clf.best_params_['max_depth'],
                    'n_estimators': clf.best_params_['n_estimators'],
                    'max_features': clf.best_params_['max_features'],
                    'random_state': 11}

reg = RandomForestRegressor(**tuned_parameters)
reg.fit(X_train, y_train)

# +
# Set figure width and height
plt.rcParams["figure.figsize"] = [5, 10]

names = X_train.columns
feature_importance = reg.feature_importances_
# make importances relative to max importance
feature_importance = 100.0 * (feature_importance / feature_importance.max())
sorted_idx = np.argsort(feature_importance)
pos = np.arange(sorted_idx.shape[0]) + .5
plt.barh(pos, feature_importance[sorted_idx], align='center')
plt.yticks(pos, map(lambda x: names[x], sorted_idx))
plt.xlabel('Relative Importance')
plt.title('Variable Importance')
plt.show()
# -

# ### Learning Curve

# +
from sklearn.model_selection import learning_curve
from sklearn.model_selection import ShuffleSplit


def plot_learning_curve(estimator, title, X, y, ylim=None, cv=None,
                        n_jobs=None, train_sizes=np.linspace(.1, 1.0, 5)):
    """
    Generate a simple plot of the test and training learning curve.

    Parameters
    ----------
    estimator : object type that implements the "fit" and "predict" methods
        An object of that type which is cloned for each validation.

    title : string
        Title for the chart.

    X : array-like, shape (n_samples, n_features)
        Training vector, where n_samples is the number of samples and
        n_features is the number of features.

    y : array-like, shape (n_samples) or (n_samples, n_features), optional
        Target relative to X for classification or regression;
        None for unsupervised learning.

    ylim : tuple, shape (ymin, ymax), optional
        Defines minimum and maximum yvalues plotted.

    cv : int, cross-validation generator or an iterable, optional
        Determines the cross-validation splitting strategy.
        Possible inputs for cv are:
          - None, to use the default 3-fold cross-validation,
          - integer, to specify the number of folds.
          - :term:`CV splitter`,
          - An iterable yielding (train, test) splits as arrays of indices.

        For integer/None inputs, if ``y`` is binary or multiclass,
        :class:`StratifiedKFold` used. If the estimator is not a classifier
        or if ``y`` is neither binary nor multiclass, :class:`KFold` is used.

        Refer :ref:`User Guide <cross_validation>` for the various
        cross-validators that can be used here.

    n_jobs : int or None, optional (default=None)
        Number of jobs to run in parallel.
        ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
        ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
        for more details.

    train_sizes : array-like, shape (n_ticks,), dtype float or int
        Relative or absolute numbers of training examples that will be used to
        generate the learning curve. If the dtype is float, it is regarded as a
        fraction of the maximum size of the training set (that is determined
        by the selected validation method), i.e. it has to be within (0, 1].
        Otherwise it is interpreted as absolute sizes of the training sets.
        Note that for classification the number of samples usually have to
        be big enough to contain at least one sample from each class.
        (default: np.linspace(0.1, 1.0, 5))
    """
    plt.figure()
    plt.title(title)
    if ylim is not None:
        plt.ylim(*ylim)
    plt.xlabel("Training examples")
    plt.ylabel("Score")
    train_sizes, train_scores, test_scores = learning_curve(
        estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes, scoring='r2')
    train_scores_mean = np.mean(train_scores, axis=1)
    train_scores_std = np.std(train_scores, axis=1)
    test_scores_mean = np.mean(test_scores, axis=1)
    test_scores_std = np.std(test_scores, axis=1)
    plt.grid()

    plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
                     train_scores_mean + train_scores_std, alpha=0.1,
                     color="r")
    plt.fill_between(train_sizes, test_scores_mean - test_scores_std,
                     test_scores_mean + test_scores_std, alpha=0.1, color="g")
    plt.plot(train_sizes, train_scores_mean, 'o-', color="r",
             label="Training score")
    plt.plot(train_sizes, test_scores_mean, 'o-', color="g",
             label="Cross-validation score")

    plt.legend(loc="best")
    return plt

title = "Learning Curves"
# Cross validation with 100 iterations to get smoother mean test and train
# score curves, each time with 20% data randomly selected as a validation set.
cv = ShuffleSplit(n_splits=100, test_size=0.2, random_state=0)

estimator = RandomForestRegressor(max_depth = clf.best_params_['max_depth'],
                                  n_estimators = clf.best_params_['n_estimators'],
                                  max_features = clf.best_params_['max_features'],
                                  random_state = 11)

plot_learning_curve(estimator, title, X_train, y_train, ylim=(0, 1.01), cv=cv, n_jobs=4)

plt.show()

# +
from sklearn.model_selection import cross_val_score

estimator = RandomForestRegressor(max_depth = clf.best_params_['max_depth'],
                                  n_estimators = clf.best_params_['n_estimators'],
                                  max_features = clf.best_params_['max_features'],
                                  random_state = 11)

scores = cross_val_score(estimator, X_train, y_train, cv=5)
print("Accuracy: %0.2f (+/- %0.2f)" % (scores.mean(), scores.std() * 2))
# -



StarcoderdataJupyterScriptsDedupFiltered