# Procesado de datos
# ==============================================================================
import numpy as np
import pandas as pd
from skforecast.datasets import fetch_dataset

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from skforecast.plot import set_dark_theme

# Forecastig
# ==============================================================================
import skforecast
import sklearn
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from feature_engine.datetime import DatetimeFeatures
from feature_engine.creation import CyclicalFeatures
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import backtesting_forecaster
from skforecast.preprocessing import RollingFeatures

# Warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')

print('Versión skforecast:', skforecast.__version__)
print('Versión sklearn:', sklearn.__version__)

# Descarga de datos
# ==============================================================================
data = fetch_dataset("ett_m1")
print(f"Data shape: {data.shape}")
data.head(3)

# Selección de datos
# ==============================================================================
data = data.loc["20170101":, ["OT"]]

# Anomalías sintéticas
# ==============================================================================
rng = np.random.default_rng(4631789)
n_anomalies = 20
anomalies_loc = data.index[rng.uniform(low=0, high=len(data), size=n_anomalies).astype(int)]
data.loc[anomalies_loc, "OT"] = data.loc[anomalies_loc, "OT"] * rng.choice([0.5, 2], size=n_anomalies)

# Gráfico
# ==============================================================================
set_dark_theme()
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(data["OT"])
ax.scatter(anomalies_loc, data.loc[anomalies_loc, "OT"], color='orange', label='Anomalies', zorder=2)
ax.set_title('OT')
ax.legend();

# Partición de datos
# ==============================================================================
end_train = '2018-04-03 23:59:00'
data_train = data.loc[: end_train, :]
data_test  = data.loc[end_train:, :]

anomalies_loc_train = [date for date in anomalies_loc if date <= pd.to_datetime(end_train)]
anomalies_loc_test  = [date for date in anomalies_loc if date > pd.to_datetime(end_train)]

print(f"Dates train : {data_train.index.min()} --- {data_train.index.max()}  (n={len(data_train)})")
print(f"    Anomalies train: {len(anomalies_loc_train)}")
print(f"Dates test  : {data_test.index.min()} --- {data_test.index.max()}  (n={len(data_test)})")
print(f"    Anomalies test : {len(anomalies_loc_test)}")

# Plot partitions
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(data_train['OT'], label='Train')
ax.plot(data_test['OT'], label='Test')
ax.scatter(anomalies_loc, data.loc[anomalies_loc, "OT"], color='orange', label='Anomalies', zorder=3)
ax.set_title('OT')
ax.legend();

# Variables de calendario
# ==============================================================================
features_to_extract = [
    'year',
    'month',
    'week',
    'day_of_week',
    'hour'
]
calendar_transformer = DatetimeFeatures(
    variables           = 'index',
    features_to_extract = features_to_extract,
    drop_original       = False,
)

# Codificación cíclica
# ==============================================================================
features_to_encode = [
    "month",
    "week",
    "day_of_week",
    "hour",
]
max_values = {
    "month": 12,
    "week": 52,
    "day_of_week": 7,
    "hour": 24,
}
cyclical_encoder = CyclicalFeatures(
                        variables     = features_to_encode,
                        max_values    = max_values,
                        drop_original = True
                   )

exog_transformer = make_pipeline(
                        calendar_transformer,
                        cyclical_encoder
                   )
display(exog_transformer)

data = exog_transformer.fit_transform(data)
# Remove rows with NaNs created by lag features
data = data.dropna()
exog_features = data.columns.difference(['OT']).tolist()
display(data.head(3))

# Forecaster
# ==============================================================================
window_features = RollingFeatures(stats=['mean', 'min', 'max'], window_sizes=24)
forecaster = ForecasterRecursive(
                 regressor       = HistGradientBoostingRegressor(random_state=15926, loss='quantile', quantile=0.5),
                 lags            = [1, 3, 4, 5, 6, 13, 14, 15, 16, 17],
                 differentiation = 1,
                 transformer_y   = StandardScaler(),
                 window_features = window_features,
             )
forecaster.fit(
    y    = data.loc[:end_train, 'OT'],
    exog = data.loc[:end_train, exog_features]
)
forecaster

# Backtesting con datos de entrenamiento
# ==============================================================================
cv = TimeSeriesFold(
        initial_train_size = None, # Utilizar todos los datos asumiendo que el modelo ya está entrenado
        steps              = 24,   # Todas las horas del siguiente día
        differentiation    = 1,
     )

metric_train, predictions_train = backtesting_forecaster(
                                    forecaster = forecaster,
                                    y          = data.loc[:end_train, 'OT'],
                                    exog       = data.loc[:end_train, exog_features],
                                    cv         = cv,
                                    metric     = 'mean_absolute_error'
                                )

results_train = predictions_train.copy()
results_train['y_true'] = data.loc[:end_train, 'OT']
results_train['residual'] = results_train['y_true'] - results_train['pred']
results_train['residual_abs'] = np.abs(results_train['residual'])
results_train.head(3)

# Gráfico predicciones, valor real y residuos
# ==============================================================================
fig, axs = plt.subplots(
    2, 1, figsize=(7, 5), gridspec_kw={"height_ratios": [2, 1]}, sharex=True
)
axs[0].plot(data.loc[:end_train, "OT"], label="Real value")
axs[0].plot(predictions_train["pred"], label="Prediction")
axs[0].scatter(
    anomalies_loc_train,
    data.loc[anomalies_loc_train, "OT"],
    zorder=2,
    color="orange",
    label="Anomalies",
)
axs[0].set_title("Backtesting on train data")
axs[0].legend()
axs[1].plot(results_train["residual_abs"], label="Absolute residuals", zorder=2)
axs[1].set_title("Absolute residuals")
top_residuals_train = (
    results_train.sort_values("residual_abs", ascending=False).head(len(anomalies_loc_train))
)
top_residuals_train["is_anomaly"] = top_residuals_train.index.isin(anomalies_loc_train)
axs[1].scatter(
    anomalies_loc_train,
    results_train.loc[anomalies_loc_train, "residual_abs"],
    zorder=2,
    color="orange",
    label="Anomalies",
)
plt.tight_layout()

display(top_residuals_train)
top_residuals_train["is_anomaly"].value_counts()

print(f"Percentage of true positives: {100 * top_residuals_train['is_anomaly'].mean():.2f}%")
print(f"Percentage of false positives: {100 * (~top_residuals_train['is_anomaly']).mean():.2f}%")

# Backtesting en test
# ==============================================================================
cv = TimeSeriesFold(
        initial_train_size = len(data_train),
        steps              = 24,  # all hours of next day
        refit              = False,
        differentiation    = 1,
     )

metric_test, predictions_test = backtesting_forecaster(
                                    forecaster = forecaster,
                                    y          = data.loc[:, 'OT'],
                                    exog       = data.loc[:, exog_features],
                                    cv         = cv,
                                    metric     = 'mean_absolute_error'
                                )
results_test = predictions_test.copy()
results_test['y_true'] = data.loc[end_train:, 'OT']
results_test['residual'] = results_test['y_true'] - results_test['pred']
results_test['residual_abs'] = np.abs(results_test['residual'])
results_test.head(3)

# Gráfico predicciones, valor real y residuos
# ==============================================================================
fig, axs = plt.subplots(
    2, 1, figsize=(8, 5), gridspec_kw={"height_ratios": [2, 1]}, sharex=True
)
axs[0].plot(data.loc[end_train:, "OT"], label="Real value")
axs[0].plot(results_test["pred"], label="Prediction")
axs[0].scatter(
    anomalies_loc_test,
    data.loc[anomalies_loc_test, "OT"],
    color="orange",
    zorder=2,
    label="Anomalies",
)
axs[0].set_title("Backtesting on test data")
axs[0].legend()
axs[1].plot(results_test["residual_abs"], label="Absolute residuals")
axs[1].set_title("Absolute residuals")
axs[1].scatter(
    anomalies_loc_test,
    results_test.loc[anomalies_loc_test, "residual_abs"],
    color="orange",
    zorder=2,
    label="Anomalies",
)
plt.tight_layout()

top_residuals_test = results_test.sort_values('residual_abs', ascending=False).head(len(anomalies_loc_test))
top_residuals_test['is_anomaly'] = top_residuals_test.index.isin(anomalies_loc_test)
display(top_residuals_test)
top_residuals_test['is_anomaly'].value_counts()

import session_info
session_info.show(html=False)