Clustering de series temporales para mejorar modelos de forecasting

# Data management
# ==============================================================================
import numpy as np
import pandas as pd
from datetime import datetime
from skforecast.datasets import fetch_dataset

# Plots
# ==============================================================================
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
plt.style.use('seaborn-v0_8-darkgrid')

# Forecasting
# ==============================================================================
import lightgbm
from lightgbm import LGBMRegressor
from sklearn.preprocessing import StandardScaler
import skforecast
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import backtesting_forecaster_multiseries

# Feature extraction
# ==============================================================================
import tsfresh
from tsfresh import extract_features
from tsfresh.feature_extraction.settings import ComprehensiveFCParameters
from tsfresh.feature_extraction.settings import from_columns

# Clustering
# ==============================================================================
import sklearn
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
import sktime
from sktime.clustering.k_means import TimeSeriesKMeans

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

color = '\033[1m\033[38;5;208m' 
print(f"{color}Versión skforecast: {skforecast.__version__}")
print(f"{color}Versión scikit-learn: {sklearn.__version__}")
print(f"{color}Versión lightgbm: {lightgbm.__version__}")
print(f"{color}Versión tsfresh: {tsfresh.__version__}")
print(f"{color}Versión sktime: {sktime.__version__}")
print(f"{color}Versión pandas: {pd.__version__}")
print(f"{color}Versión numpy: {np.__version__}")

# Lectura de datos
# ==============================================================================
data = fetch_dataset('ashrae_daily')
data.head()

# Imputación de valores faltantes de air_temperature y wind_speed usando fowrard y backward fill
# ==============================================================================
# La imputación debe hacerse de forma independiente para cada serie temporal
data = data.sort_values(by=['building_id', 'timestamp'])
data['air_temperature'] = data.groupby('building_id')['air_temperature'].ffill().bfill()
data['wind_speed'] = data.groupby('building_id')['wind_speed'].ffill().bfill()
data = data.sort_index()

print(
    f"Rage of dates available : {data.index.min()} --- {data.index.max()}  "
    f"(n_days={(data.index.max() - data.index.min()).days})"
)

# Número de edificios y tipo de edificios basado en el uso principal
# ==============================================================================
n_building = data['building_id'].nunique()
n_type_building = data['primary_use'].nunique()
range_datetime = [data.index.min(), data.index.max()]

print(f"Number of buildings: {n_building}")
print(f"Number of building types: {n_type_building}")
display(data.drop_duplicates(subset=['building_id'])['primary_use'].value_counts())

# Tipos de edificios (primary use) con menos de 50 muestras se agrupan como "Other".
# ==============================================================================
infrequent_categories = (
    data
    .drop_duplicates(subset=['building_id'])['primary_use']
    .value_counts()
    .loc[lambda x: x < 50]
    .index
    .tolist()
)
print("Infrequent categories:")
print("======================")
print('\n'.join(infrequent_categories))

data['primary_use'] = np.where(
    data['primary_use'].isin(infrequent_categories),
    'Other',
    data['primary_use']
)

# Serie temporal para 1 edificio seleccionado al azar por grupo
# ==============================================================================
fig, axs = plt.subplots(nrows=3, ncols=2, figsize=(8, 5.5), sharex=True, sharey=False)
sample_ids = (
    data
    .groupby('primary_use')['building_id']
    .apply(lambda x: x.sample(1, random_state=333))
    .tolist()
)
axs = axs.flatten()

for i, building_id in enumerate(sample_ids):
    data_sample = data[data['building_id'] == building_id]
    building_type = data_sample['primary_use'].unique()[0]
    data_sample.plot(
        y        = 'meter_reading',
        ax       = axs[i],
        legend   = False,
        title    = f"Building: {building_id}, type: {building_type}",
        fontsize = 8
    )
    axs[i].set_xlabel("")
    axs[i].set_ylabel("")
    # Scientific notation for y axis
    axs[i].ticklabel_format(axis='y', style='sci', scilimits=(0, 0))
    axs[i].title.set_size(9)

fig.suptitle('Energy consumption for 6 random buildings', fontsize=12)
fig.tight_layout()
plt.show()

# Consumo energético por tipo de edificio (una línea gris por edificio)
# ==============================================================================
fig, axs = plt.subplots(nrows=3, ncols=2, figsize=(8, 5.5), sharex=True, sharey=False)
axs = axs.flatten()

for i, building_type in enumerate(data['primary_use'].unique()):
    data_sample = data[data['primary_use'] == building_type]
    data_sample = data_sample.pivot_table(
                      index   = 'timestamp',
                      columns = 'building_id',
                      values  = 'meter_reading',
                      aggfunc = 'mean'
                  )
    data_sample.plot(
        legend   = False,
        title    = f"Type: {building_type}",
        color    = 'gray',
        alpha    = 0.2,
        fontsize = 8,
        ax       = axs[i]
    )
    data_sample.mean(axis=1).plot(
        ax       = axs[i],
        color    = 'blue',
        fontsize = 8
    )
    axs[i].set_xlabel("")
    axs[i].set_ylabel("")
    # Scientific notation for y axis
    axs[i].ticklabel_format(axis='y', style='sci', scilimits=(0,0))
    axs[i].title.set_size(9)

    # Limit the axis to 5 times the maximum mean value to improve visualization
    axs[i].set_ylim(
        bottom = 0,
        top    = 5 * data_sample.mean(axis=1).max()
    )

fig.tight_layout()
fig.suptitle('Energy consumption by type of building', fontsize=12)
plt.subplots_adjust(top=0.9)
plt.show()

# Configuración por defecto de tsfresh
# ==============================================================================
default_features = ComprehensiveFCParameters()
print("Name of the features extracted by tsfresh:")
print("=========================================")
print('\n'.join(list(default_features)))

# Configuración por defecto para "partial_autocorrelation"
# ==============================================================================
default_features['partial_autocorrelation']

# Configuración por defecto para todas las características extraídas
# ==============================================================================
# from pprint import pprint
# pprint(default_features)

# Feature extraction
# ==============================================================================
ts_features = extract_features(
    timeseries_container  = data[['building_id', 'meter_reading']].reset_index(),
    column_id             = "building_id",
    column_sort           = "timestamp",
    column_value          = "meter_reading",
    default_fc_parameters = default_features,
    impute_function       = tsfresh.utilities.dataframe_functions.impute,
    n_jobs                = 4
)

print("Shape of ts_features:", ts_features.shape)
ts_features.head(3)

# Eliminar características con valores faltantes
# ==============================================================================
ts_features = ts_features.dropna(axis=1, how='any')
print(f"Number of features after selection and removing missing values: {ts_features.shape[1]}")

# Diccionario con las características seleccionadas y su configuración
# ==============================================================================
ts_features_info = from_columns(ts_features)
# pprint(ts_features_info['meter_reading'])

# Escalar características para que todas tengan media 0 y desviación estándar 1
# ==============================================================================
scaler = StandardScaler().set_output(transform="pandas")
ts_features_scaled = scaler.fit_transform(ts_features)
ts_features_scaled.head(2)

# Gráfico de reducción de varianza en función del número de componentes PCA
# ==============================================================================
pca = PCA()
pca.fit(ts_features_scaled)
n_components = np.argmax(np.cumsum(pca.explained_variance_ratio_) > 0.85) + 1

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 2.5))
ax.plot(np.cumsum(pca.explained_variance_ratio_))
ax.set_title('Variance explained by PCA components')
ax.set_xlabel('Number of components')
ax.set_ylabel('Cumulative explained variance')
ax.axhline(y=0.85, color='red', linestyle='--')
ax.axvline(x=n_components, color='red', linestyle='--')
ax.text(
    x=n_components + 10,
    y=0.4,
    s=f"n_components = {n_components}",
    color='red',
    verticalalignment='center'
)

plt.show();

# PCA con tantos componentes como sea necesario para explicar el 85% de la varianza
# ==============================================================================
pca = PCA(n_components=0.85)
pca_projections = pca.fit_transform(ts_features_scaled)

# Create a data frame with the projections. Each column is a principal component
# and each row is the id of the building.
pca_projections = pd.DataFrame(
    pca_projections,
    index   = ts_features_scaled.index,
    columns = [f"PC{i}" for i in range(1, pca.n_components_ + 1)]
)

print(f"Number of components selected: {pca.n_components_}")
print(f"Explained variance: {pca.explained_variance_ratio_.sum():.3f}")
pca_projections.head(3)

# Número óptimo de clusters usando el método del codo
# ==============================================================================
range_n_clusters = range(1, 30)
inertias = []
size_of_clusters = []

for n_clusters in range_n_clusters:
    kmeans = KMeans(
                 n_clusters   = n_clusters, 
                 n_init       = 20, 
                 random_state = 963852
             )
    kmeans.fit(pca_projections)
    inertias.append(kmeans.inertia_)
    size_of_clusters.append(np.unique(kmeans.labels_, return_counts=True)[1].tolist())

inertias = pd.Series(inertias, index=range_n_clusters)
dicrease_in_inertia = inertias.pct_change() * -100

fig, axs = plt.subplots(1, 2, figsize=(9, 2.5))
inertias.plot(marker='o', ax=axs[0])
axs[0].xaxis.set_major_locator(ticker.MaxNLocator(integer=True))
axs[0].set_title("Cluster inertia vs number of clusters")
axs[0].set_xlabel('Number of clusters')
axs[0].set_ylabel('Intra-cluster variance (inertia)')

dicrease_in_inertia.plot(kind='bar', ax=axs[1])
axs[1].set_title("Decrease in inertia")
axs[1].set_xlabel('Number of clusters')
axs[1].set_ylabel('Decrease in inertia (%)')
fig.tight_layout()
plt.show();

# Distribución de tamaños de los clusters
# ==============================================================================
for n_clusters, sizes in zip(range_n_clusters, size_of_clusters):
    print(f"Size of clusters (n = {n_clusters}): {sizes}")

# Entrenar el modelo de clustering con 7 clusters y asignar cada edificio a un cluster
# ==============================================================================
kmeans = KMeans(n_clusters=7, n_init=20, random_state=963852)
kmeans.fit(pca_projections)
clusters = kmeans.predict(pca_projections)
clusters = pd.DataFrame({
    'building_id': pca_projections.index,
    'cluster_base_on_features': clusters.astype(str)
})

# Merge cluster minor clusters into a single cluster
threshold = 20
cluster_size = clusters['cluster_base_on_features'].value_counts()
cluster_size = cluster_size[cluster_size < threshold].index.tolist()
clusters['cluster_base_on_features'] = np.where(
    clusters['cluster_base_on_features'].isin(cluster_size),
    'Other',
    clusters['cluster_base_on_features']
)
display(clusters.head(3))
clusters['cluster_base_on_features'].value_counts().to_frame(name='Number of series')

# Añadir el cluster predicho al dataframe con la información de los edificios
# ==============================================================================
data = pd.merge(
           data.reset_index(),  # To avoid losing the index
           clusters,
           on       ='building_id',
           how      ='left',
           validate = 'm:1'
       ).set_index('timestamp')

data.head(3)

# Porcentaje de tipos de edificios por cluster
# ==============================================================================
primary_usage_per_cluster = (
    data
    .groupby('cluster_base_on_features')['primary_use']
    .value_counts(normalize=True)
    .unstack()
    .fillna(0)
)

primary_usage_per_cluster = (100 * primary_usage_per_cluster).round(1)
primary_usage_per_cluster

# Convertir los datos a formato long con un multiíndice: (building_id, timestamp)
# ==============================================================================
data_long = (
    data
    .reset_index()
    .loc[:, ['timestamp', 'building_id', 'meter_reading']]
    .set_index(['building_id', 'timestamp'])
    .sort_index(ascending=True)
)
data_long

# Fit clustering
# ==============================================================================
model = TimeSeriesKMeans(
            n_clusters   = 4,
            metric       = "dtw",
            max_iter     = 10,
            random_state = 123456
        )

model.fit(data_long)

# Predocción de clusters
# ==============================================================================
clusters = model.predict(data_long)
clusters = pd.DataFrame({
               'building_id': data_long.index.get_level_values('building_id').drop_duplicates(),
               'cluster_base_on_dtw': clusters.astype(str),
           })

# Size of each cluster
clusters['cluster_base_on_dtw'].value_counts().sort_index().to_frame(name='Number of series')

# Añadir el cluster predicho al dataframe con la información de los edificios
# ==============================================================================
data = pd.merge(
           data.reset_index(),  # To avoid losing the index
           clusters,
           on       = 'building_id',
           how      = 'left',
           validate = 'm:1'
       ).set_index('timestamp')

data.head(3)

# Guardar en disco los datos para modelado para evitar repetir los pasos anteriores
# ==============================================================================
data.to_parquet('data_modelling.parquet')

# Lectura de datos para modelado
# ==============================================================================
data = pd.read_parquet('data_modelling.parquet')

# División de los datos en entrenamiento, validación y test
# ==============================================================================
end_train = '2016-07-31 23:59:00'
end_validation = '2016-09-30 23:59:00'
data_train = data.loc[: end_train, :].copy()
data_val   = data.loc[end_train:end_validation, :].copy()
data_test  = data.loc[end_validation:, :].copy()

print(
    f"Rage of dates available : {data.index.min()} --- {data.index.max()} "
    f"(n_days={(data.index.max() - data.index.min()).days})"
)
print(
    f"  Dates for training    : {data_train.index.min()} --- {data_train.index.max()} "
    f"(n_days={(data_train.index.max() - data_train.index.min()).days})"
)
print(
    f"  Dates for validation  : {data_val.index.min()} --- {data_val.index.max()} "
    f"(n_days={(data_val.index.max() - data_val.index.min()).days})"
)
print(
    f"  Dates for test        : {data_test.index.min()} --- {data_test.index.max()} "
    f"(n_days={(data_test.index.max() - data_test.index.min()).days})"
)

# Tabla de resultados
# ==============================================================================
table_results = pd.DataFrame(columns=['model', 'mae', 'abs_error', 'bias', 'elapsed_time'])
table_results = table_results.set_index('model')
table_results = table_results.astype({'mae': float, 'abs_error': float, 'bias': float, 'elapsed_time': object})

# Función para añadir información del día de la semana
# ==============================================================================
def day_of_week_cyclical_encoding(data):
    """
    Calculate cyclical encoding for the day of week using a 
    DataFrame with a datetime index.
    """
    data['day_of_week'] = data.index.dayofweek
    data['sin_day_of_week'] = np.sin(2*np.pi*data['day_of_week']/7)
    data['cos_day_of_week'] = np.cos(2*np.pi*data['day_of_week']/7)
    
    return data

# Variables exógenas incluidas en el modelo
# ==============================================================================
exog_features = ['sin_day_of_week', 'cos_day_of_week', 'air_temperature', 'wind_speed']

# Transformar los datos para tener una columna por edificio
# ==============================================================================
data_pivot = data.pivot_table(
                 index    = 'timestamp',
                 columns  = 'building_id',
                 values   = 'meter_reading',
                 aggfunc  = 'mean' 
             )
data_pivot.columns = data_pivot.columns.astype(str)
data_pivot = data_pivot.asfreq('D').sort_index()

# Add calendar features
data_pivot = day_of_week_cyclical_encoding(data_pivot)

# Add temperature and wind speed as exogenous variables
data_pivot = data_pivot.merge(
                 data[['air_temperature', 'wind_speed']].resample('D').mean(),
                 left_index  = True,
                 right_index = True,
                 how         = 'left',
                 validate    = '1:m'
             )

data_pivot.head(3)

# Forecaster multiseries para modelar todos los edificios a la vez
# ==============================================================================
exog_features = ['sin_day_of_week', 'cos_day_of_week', 'air_temperature', 'wind_speed']

params_lgbm = {
    'n_estimators': 500,
    'learning_rate': 0.01,
    'max_depth': 10,
    'random_state': 8520,
    'verbose': -1
}
forecaster = ForecasterRecursiveMultiSeries(
                 regressor = LGBMRegressor(**params_lgbm),
                 lags      = 31,
                 encoding  = 'ordinal_category'
             )

start = datetime.now()

cv = TimeSeriesFold(
        steps              = 7,
        initial_train_size = len(data_pivot.loc[:end_validation, :]),
        refit              = False
     )
metric, predictions = backtesting_forecaster_multiseries(
                          forecaster            = forecaster,
                          series                = data_pivot.filter(regex='^id_'), # Select only buildings
                          exog                  = data_pivot[exog_features],
                          cv                    = cv,
                          metric                = 'mean_absolute_error',
                          add_aggregated_metric = False,
                          verbose               = False,
                          show_progress         = True
                      )

end = datetime.now()

mean_metric_all_buildings = metric['mean_absolute_error'].mean()
errors_all_buildings = predictions - data_pivot.loc[predictions.index, predictions.columns]
sum_abs_errors_all_buildings = errors_all_buildings.abs().sum().sum()
sum_bias_all_buildings = errors_all_buildings.sum().sum()

table_results.loc['Global model', 'mae'] = mean_metric_all_buildings
table_results.loc['Global model', 'abs_error'] = sum_abs_errors_all_buildings
table_results.loc['Global model', 'bias'] = sum_bias_all_buildings
table_results.loc['Global model', 'elapsed_time'] = end - start

print("")
print(f"Average mean absolute error for all buildings: {mean_metric_all_buildings:.0f}")
print(f"Sum of absolute errors for all buildings (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")

# Gráfico de predicciones vs valor real para 2 edificios seleccionados al azar
# ==============================================================================
rng = np.random.default_rng(147)
selected_buildings = rng.choice(data['building_id'].unique(), size=2, replace=False)

fig, axs = plt.subplots(2, 1, figsize=(6, 4), sharex=True)
axs = axs.flatten()

for i, building in enumerate(selected_buildings):
    data_test.query("building_id == @building")['meter_reading'].plot(ax=axs[i], label='test')
    predictions[building].plot(ax=axs[i], label='One model per building')
    axs[i].set_title(f"Building {building}", fontsize=10)
    axs[i].set_xlabel("")
    axs[i].legend()

fig.tight_layout()
plt.show();

# Forecaster multi-series para los edificios agrupados por uso principal
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []

params_lgbm = {
    'n_estimators': 500,
    'learning_rate': 0.01,
    'max_depth': 10,
    'random_state': 8520,
    'verbose': -1
}
forecaster = ForecasterRecursiveMultiSeries(
                 regressor = LGBMRegressor(**params_lgbm),
                 lags      = 31,
                 encoding  = "ordinal_category"
             )

start = datetime.now()

for primary_usage in data['primary_use'].unique():

    print(
        f"Training and testing model for primary usage: {primary_usage} "
        f"(n = {data[data['primary_use'] == primary_usage]['building_id'].nunique()})"
    )

    # Create a subset based on primary use clusters
    data_subset = data[data['primary_use'] == primary_usage]
    data_subset = data_subset.pivot_table(
                      index   = 'timestamp',
                      columns = 'building_id',
                      values  = 'meter_reading',
                      aggfunc = 'mean'
                  )
    data_subset.columns = data_subset.columns.astype(str)
    data_subset = data_subset.asfreq('D').sort_index()

    # Add calendar features
    data_subset = day_of_week_cyclical_encoding(data_subset)
    data_subset = data_subset.merge(
                      data[['air_temperature', 'wind_speed']].resample('D').mean(),
                      left_index  = True,
                      right_index = True,
                      how         = 'left',
                      validate    = '1:m'
                  )

    cv = TimeSeriesFold(
            steps              = 7,
            initial_train_size = len(data_subset.loc[:end_validation, :]),
            refit              = False
         )
    metric, predictions = backtesting_forecaster_multiseries(
                              forecaster            = forecaster,
                              series                = data_subset.filter(regex='^id_'),
                              exog                  = data_subset[['sin_day_of_week', 'cos_day_of_week']],
                              cv                    = cv,
                              metric                = 'mean_absolute_error',
                              add_aggregated_metric = False,
                              verbose               = False,
                              show_progress         = False
                          )
    predictions_all_buildings.append(predictions)
    metrics_all_buildings.append(metric)

end = datetime.now()

predictions_all_buildings = pd.concat(predictions_all_buildings, axis=1)
metrics_all_buildings = pd.concat(metrics_all_buildings, axis=0)
errors_all_buildings = (
    predictions_all_buildings - data_pivot.loc[predictions_all_buildings.index, predictions_all_buildings.columns]
)
mean_metric_all_buildings = metrics_all_buildings['mean_absolute_error'].mean()
sum_abs_errors_all_buildings = errors_all_buildings.abs().sum().sum()
sum_bias_all_buildings = errors_all_buildings.sum().sum()

table_results.loc['Global model per primary usage', 'mae'] = mean_metric_all_buildings
table_results.loc['Global model per primary usage', 'abs_error'] = sum_abs_errors_all_buildings
table_results.loc['Global model per primary usage', 'bias'] = sum_bias_all_buildings
table_results.loc['Global model per primary usage', 'elapsed_time'] = end - start

print("")
print(f"Average mean absolute error for all buildings: {mean_metric_all_buildings:.0f}")
print(f"Sum of absolute errors for all buildings (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")

# Añadir predicciones al gráfico ya existente (no mostrar el gráfico todavía)
# ==============================================================================
for i, building in enumerate(selected_buildings):
    predictions_all_buildings[building].plot(ax=axs[i], label='Global model per primary usage')
    axs[i].legend()

# Forecaster multiserie para los edificios agrupados por características de las series
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []

params_lgbm = {
    'n_estimators': 500,
    'learning_rate': 0.01,
    'max_depth': 10,
    'random_state': 8520,
    'verbose': -1
}
forecaster = ForecasterRecursiveMultiSeries(
                 regressor = LGBMRegressor(**params_lgbm),
                 lags      = 31,
                 encoding  = "ordinal_category"
             )

start = datetime.now()

for cluster in data['cluster_base_on_features'].unique():
    
    print(
        f"Training and testing model for cluster: {cluster} "
        f"(n = {data[data['cluster_base_on_features'] == cluster]['building_id'].nunique()})"
    )

    # Create subset based on features clusters
    data_subset = data[data['cluster_base_on_features'] == cluster]
    data_subset = data_subset.pivot_table(
                        index   = 'timestamp',
                        columns = 'building_id',
                        values  = 'meter_reading',
                        aggfunc = 'mean'
                  )
    data_subset.columns = data_subset.columns.astype(str)
    data_subset = data_subset.asfreq('D').sort_index()

    # Add calendar features
    data_subset = day_of_week_cyclical_encoding(data_subset)
    data_subset = data_subset.merge(
                      data[['air_temperature', 'wind_speed']].resample('D').mean(),
                      left_index  = True,
                      right_index = True,
                      how         = 'left',
                      validate    = '1:m'
                  )
    cv = TimeSeriesFold(
            steps              = 7,
            initial_train_size = len(data_subset.loc[:end_validation, :]),
            refit              = False
         )
    metric, predictions = backtesting_forecaster_multiseries(
                              forecaster            = forecaster,
                              series                = data_subset.filter(regex='^id_'),
                              exog                  = data_subset[['sin_day_of_week', 'cos_day_of_week']],
                              cv                    = cv,
                              metric                = 'mean_absolute_error',
                              add_aggregated_metric = False,
                              verbose               = False,
                              show_progress         = False
                          )
    predictions_all_buildings.append(predictions)
    metrics_all_buildings.append(metric)

end = datetime.now()

predictions_all_buildings = pd.concat(predictions_all_buildings, axis=1)
metrics_all_buildings = pd.concat(metrics_all_buildings, axis=0)
errors_all_buildings = (
    predictions_all_buildings - data_pivot.loc[predictions_all_buildings.index, predictions_all_buildings.columns]
)
mean_metric_all_buildings = metrics_all_buildings['mean_absolute_error'].mean()
sum_abs_errors_all_buildings = errors_all_buildings.abs().sum().sum()
sum_bias_all_buildings = errors_all_buildings.sum().sum()

table_results.loc['Global model per cluster (features)', 'mae'] = mean_metric_all_buildings
table_results.loc['Global model per cluster (features)', 'abs_error'] = sum_abs_errors_all_buildings
table_results.loc['Global model per cluster (features)', 'bias'] = sum_bias_all_buildings
table_results.loc['Global model per cluster (features)', 'elapsed_time'] = end - start

print("")
print(f"Average mean absolute error for all buildings: {mean_metric_all_buildings:.0f}")
print(f"Sum of absolute errors for all buildings (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")

# Añaadir predicciones al gráfico ya existente (no mostrar el gráfico todavía)
# ==============================================================================
for i, building in enumerate(selected_buildings):
    predictions_all_buildings[building].plot(ax=axs[i], label='Global model per cluster (features)')
    axs[i].legend()

# Forecaster multiserie para los edificios agrupados por DTW
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []

params_lgbm = {
    'n_estimators': 500,
    'learning_rate': 0.01,
    'max_depth': 10,
    'random_state': 8520,
    'verbose': -1
}
forecaster = ForecasterRecursiveMultiSeries(
                 regressor = LGBMRegressor(**params_lgbm),
                 lags      = 31,
                 encoding  = 'ordinal_category'
             )

start = datetime.now()

for cluster in data['cluster_base_on_dtw'].unique():

    print(
        f"Training and testing model for cluster: {cluster} "
        f"(n = {data[data['cluster_base_on_dtw'] == cluster]['building_id'].nunique()})"
    )

    # Create subset based on DTW clusters
    data_subset = data[data['cluster_base_on_dtw'] == cluster]
    data_subset = data_subset.pivot_table(
                      index   = 'timestamp',
                      columns = 'building_id',
                      values  = 'meter_reading',
                      aggfunc = 'mean'
                  )
    data_subset.columns = data_subset.columns.astype(str)
    data_subset = data_subset.asfreq('D').sort_index()

    # Add calendar features
    data_subset = day_of_week_cyclical_encoding(data_subset)
    data_subset = data_subset.merge(
                      data[['air_temperature', 'wind_speed']].resample('D').mean(),
                      left_index  = True,
                      right_index = True,
                      how         = 'left',
                      validate    = '1:m'
                  )
    cv = TimeSeriesFold(
            steps              = 7,
            initial_train_size = len(data_subset.loc[:end_validation, :]),
            refit              = False
        )
    metric, predictions = backtesting_forecaster_multiseries(
                              forecaster            = forecaster,
                              series                = data_subset.filter(regex='^id_'),
                              exog                  = data_subset[['sin_day_of_week', 'cos_day_of_week']],
                              cv                    = cv,
                              metric                = 'mean_absolute_error',
                              add_aggregated_metric = False,
                              verbose               = False,
                              show_progress         = False
                          )
    predictions_all_buildings.append(predictions)
    metrics_all_buildings.append(metric)

end = datetime.now()

predictions_all_buildings = pd.concat(predictions_all_buildings, axis=1)
metrics_all_buildings = pd.concat(metrics_all_buildings, axis=0)
errors_all_buildings = (
    predictions_all_buildings - data_pivot.loc[predictions_all_buildings.index, predictions_all_buildings.columns]
)
mean_metric_all_buildings = metrics_all_buildings['mean_absolute_error'].mean()
sum_abs_errors_all_buildings = errors_all_buildings.abs().sum().sum()
sum_bias_all_buildings = errors_all_buildings.sum().sum()

table_results.loc['Global model per cluster (DTW)', 'mae'] = mean_metric_all_buildings
table_results.loc['Global model per cluster (DTW)', 'abs_error'] = sum_abs_errors_all_buildings
table_results.loc['Global model per cluster (DTW)', 'bias'] = sum_bias_all_buildings
table_results.loc['Global model per cluster (DTW)', 'elapsed_time'] = end - start

print("")
print(f"Average mean absolute error for all buildings: {mean_metric_all_buildings:.0f}")
print(f"Sum of absolute errors for all buildings (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")

# Añaadir predicciones al gráfico ya existente (no mostrar el gráfico todavía)
# ==============================================================================
for i, building in enumerate(selected_buildings):
    predictions_all_buildings[building].plot(ax=axs[i], label='Global model per cluster (DTW)')
    axs[i].legend()

# Tabla de resultados
# ==============================================================================
table_results['elapsed_time'] = table_results['elapsed_time'].astype(str).str[:7]
table_results.style.highlight_min(axis=0, color='green').format(precision=0)

# Plot predictions vs real value for 2 random buildings for all models
# ==============================================================================
handles, labels = axs[0].get_legend_handles_labels()
fig.legend(handles, labels, loc='upper center', bbox_to_anchor=(0.5, 0.05), ncol=2)
for ax in axs:
    ax.legend().remove()
fig

import session_info
session_info.show(html=False)