Examples

Some didactic examples can be found in the mbtr/examples dir. All the examples are based on the Hierarchical Demand Forecasting Benchmark, which is downloaded at the beginning of the examples. The dataset is downloaded only once, successive calls to mbtr.ut.load_dataset() will only read the locally downloaded file.

import numpy as np
import mbtr.utils as ut
from scipy.linalg import hankel
import matplotlib.pyplot as plt
from mbtr.mbtr import MBT
from mbtr.utils import set_figure

# --------------------------- Download and format data ----------------------------------------------------------------
# download power data from "Hierarchical Demand Forecasting Benchmark for the Distribution Grid" dataset
# from https://zenodo.org/record/3463137#.XruXbx9fiV7 if needed
power_data = ut.load_dataset()

Multivariate forecasts and linear regression

The examples/multivariate_forecast.py shows an example of usage of the mbtr.losses.MSE and mbtr.losses.LinRegLoss losses.

We start creating the training and test set. We use the P_mean signal as a target, and we want to predict the next day ahead only using past values from the same time series. We start downsampling the signal to one hour, and then we embed it in a 2-days matrix. The first 24 columns refers to the previous day in time, and are the signals which will be used to predict the last 24 columns:

total_power = power_data['P_mean']['all'].values.reshape(-1,1)
# down-sample it to 1 hour, bin averages
total_power = np.mean(total_power[:len(total_power)-len(total_power) % 6].reshape(-1, 6), axis=1, keepdims=True)

# embed the signal in a 2-days matrix
total_power = hankel(total_power, np.zeros((2 * 24, 1)))

# create feature matrix and target for the training and test sets
x = total_power[:, :24]
y = total_power[:, 24:]
n_tr = int(len(x)*0.8)

x_tr, y_tr, x_te, y_te = [x[:n_tr, :], y[:n_tr, :], x[n_tr:, :], y[n_tr:, :]]

we can perform a visual check on the features and target signals:

# visual check on the first 50 samples of features and targets
fig, ax = set_figure((5,4))
y_min = np.min(y_tr[:50, :]) * 0.9
y_max = np.max(y_tr[:50, :]) * 1.1

for i in range(50):
    ax.cla()
    ax.plot(np.arange(24), x_tr[i,:], label='features')
    ax.plot(np.arange(24) + 24, y_tr[i, :], label='multivariate targets')
    ax.set_xlabel('step ahead [h]')
    ax.set_ylabel('P [kW]')
    ax.legend(loc='upper right')
    ax.set_ylim(y_min, y_max)
    plt.pause(1e-6)
plt.close('all')

Now we are ready to fit an MBT instance. We start fitting a mean squared error loss function:

print('#'*20 + '    Fitting MBT with mse loss   ' + '#'*20)
m = MBT(n_boosts=30,  min_leaf=100, lambda_weights=1e-3).fit(x_tr, y_tr, do_plot=True)
y_hat = m.predict(x_te)

As a comparison, we also fit also 24 different LightGBM instances with the utility class mbtr.ut.LightGBMMISO :

print('#'*20 + '    Fitting 24 MISO LightGBMs   ' + '#'*20)

m_lgb = ut.LightGBMMISO(30).fit(x_tr, y_tr)
y_hat_lgb = m_lgb.predict(x_te)

As a last comparison, we fit a linear response MBT, using the mbtr.losses.LinRegLoss. This class requires as additional input a set of features for fitting the linear response inside each leaf. In order to reduce the computational time, we only use the mean, maximum, minimum and the first and last values of the original regressors matrix x as features for finding the best splits of the trees.

x_build = np.hstack([np.mean(x, axis=1,keepdims=True), np.max(x,axis=1, keepdims=True),
                     np.min(x, axis=1, keepdims=True), x[:,[0,23]]])
x_build_tr, x_build_te = [x_build[:n_tr, :],x_build[n_tr:, :]]
m_lin = MBT(loss_type='linear_regression', n_boosts=30,  min_leaf=1500,
            lambda_weights=1e-3).fit(x_build_tr, y_tr,x_lr=x_tr, do_plot=True)
y_hat_lin = m_lin.predict(x_build_te, x_lr=x_te)

We can now plot some results:

fig, ax = set_figure((5,4))
y_min = np.min(y_tr[:150, :]) * 0.9
y_max = np.max(y_tr[:150, :]) * 1.1


for i in range(150):
    ax.cla()
    ax.plot(np.arange(24), y_te[i,:], label='test')
    ax.plot(y_hat[i, :], '--', label='mbtr')
    ax.plot(y_hat_lgb[i, :], '--', label='lgb')
    ax.plot(y_hat_lin[i, :], '--', label='mbtr-lin')
    ax.set_xlabel('step ahead [h]')
    ax.set_ylabel('P [kW]')
    ax.legend(loc='upper right')
    ax.set_ylim(y_min, y_max)
    plt.pause(1e-6)

and compare the models in term of RMSE:

mean_rmse = lambda x,y : np.mean(np.mean((x-y)**2, axis=1)**0.5)
rmse_mbt = mean_rmse(y_te, y_hat)
rmse_lgb = mean_rmse(y_te, y_hat_lgb)
rmse_mbt_lin = mean_rmse(y_te, y_hat_lin)

print('#'*20 + '    Mean-horizon RMSEs  ' + '#'*20 )
[print('{}: {:0.2e}'.format(n, s)) for n, s in zip(['mbtr', 'lgb', 'mbtr-lin'], [rmse_mbt, rmse_lgb, rmse_mbt_lin])]

Time smoothing and Fourier regression

The examples/fourier_and_smoothing.py shows an example of usage of the mbtr.losses.TimeSmoother and mbtr.losses.FourierLoss losses. The first part of the code is identical to the one used in the examples/multivariate_forecast.py example; we download the dataset and create the training and test sets. We now use the time_smoother loss function, which penalize the second order discrete derivative of the response function:

print('#'*20 + '    Fitting MBT with smooth loss   ' + '#'*20)
m_sm = MBT(loss_type='time_smoother', lambda_smooth=1, n_boosts=30,
           min_leaf=300, lambda_weights=1e-3).fit(x_tr, y_tr, do_plot=True)
y_hat_sm = m_sm.predict(x_te)

Keeping all the other MBT parameters unchanged, we can fit two Fourier losses with different number of harmonics:

print('#'*20 + '    Fitting MBT with Fourier loss and 3 harmonics    ' + '#'*20)

m_fou_3 = MBT(loss_type='fourier', n_harmonics=3, n_boosts=30,
              min_leaf=300, lambda_weights=1e-3).fit(x_tr, y_tr, do_plot=True)
y_hat_fou_3 = m_fou_3.predict(x_te)

print('#'*20 + '    Fitting MBT with Fourier loss and 5 harmonics    ' + '#'*20)

m_fou_5 = MBT(loss_type='fourier', n_harmonics=5, n_boosts=30,  min_leaf=300,
              lambda_weights=1e-3).fit(x_tr, y_tr, do_plot=True)
y_hat_fou_5 = m_fou_5.predict(x_te)

We can now plot some results from the different fitted losses:

fig, ax = set_figure((5,4))
y_min = np.min(y_te[:150, :]) * 0.9
y_max = np.max(y_te[:150, :]) * 1.1

for i in range(150):
    ax.cla()
    ax.plot(np.arange(24), y_te[i,:], label='test')
    ax.plot(y_hat_sm[i, :], '--', label='time-smoother')
    ax.plot(y_hat_fou_3[i, :], '--', label='fourier-3')
    ax.plot(y_hat_fou_5[i, :], '--', label='fourier-5')

    ax.set_xlabel('step ahead [h]')
    ax.set_ylabel('P [kW]')
    ax.legend(loc='upper right')
    ax.set_ylim(y_min, y_max)
    plt.pause(1e-6)

Finally, we can compare the models in term of RMSE:

mean_rmse = lambda x, y: np.mean(np.mean((x - y) ** 2, axis=1) ** 0.5)
rmse_sm = mean_rmse(y_te, y_hat_sm)
rmse_fou_3 = mean_rmse(y_te, y_hat_fou_3)
rmse_fou_5 = mean_rmse(y_te, y_hat_fou_5)


print('#' * 20 + '    Mean-horizon RMSEs  ' + '#' * 20)
[print('{}: {:0.2e}'.format(n, s)) for n, s in zip(['smoother', 'fourier-3', 'fourier-5'], [rmse_sm, rmse_fou_3, rmse_fou_5])]

Quantile loss

The examples/quantiles.py shows an example of usage of the mbtr.losses.QuantileLoss loss. In this example we aim at predicting the quantiles of the next step ahead, using the previous 24 hours of the signal as covariates. After downloading the dataset as described in the previous example, we build the training and test sets:

# embed the signal in a 2-days matrix
total_power = hankel(total_power, np.zeros((25, 1)))[:-25, :]

# create feature matrix and target for the training and test sets
x = total_power[:, :24]
y = total_power[:, 24:]
n_tr = int(len(x)*0.8)
x_tr, y_tr, x_te, y_te = [x[:n_tr, :], y[:n_tr, :], x[n_tr:, :], y[n_tr:, :]]

we plot some training instances of the features and the target to have a visual check:

fig, ax = set_figure((5,4))
y_min = np.min(y_tr[:50, :]) * 0.9
y_max = np.max(y_tr[:50, :]) * 1.1

for i in range(50):
    ax.cla()
    ax.plot(np.arange(24), x_tr[i,:], label='features')
    ax.scatter(25, y_tr[i, :], label='multivariate targets', marker='.')
    ax.set_xlabel('step ahead [h]')
    ax.set_ylabel('P [kW]')
    ax.legend(loc='upper right')
    ax.set_ylim(y_min, y_max)
    plt.pause(1e-6)
plt.close('all')

Finally, we can train a MBT instance with a mbtr.losses.QuantileLoss loss. Note that this loss requires the alphas additional parameter. This is an array of quantiles to be fitted:

alphas = np.linspace(0.05, 0.95, 7)
m = MBT(loss_type='quantile', alphas=alphas, n_boosts=40,
        min_leaf=300, lambda_weights=1e-3).fit(x_tr, y_tr, do_plot=True)

At last, we can plot some predictions for the required quantiles:

y_hat = m.predict(x_te)
fig, ax = set_figure((5,4))
y_min = np.min(y_tr[:50, :]) * 0.9
y_max = np.max(y_tr[:50, :]) * 1.1
n_q = y_hat.shape[1]
n_sa = y_te.shape[1]
n_plot = 300
colors = plt.get_cmap('plasma', int(n_q))
for fl in np.arange(np.floor(n_q / 2), dtype=int):
    q_low = np.squeeze(y_hat[:n_plot, fl])
    q_up = np.squeeze(y_hat[:n_plot, n_q - fl - 1])
    x = np.arange(len(q_low))
    ax.fill_between(x, q_low, q_up, color=colors(fl), alpha=0.1 + 0.6*fl/n_q, linewidth=0.0)
plt.plot(y_te[:n_plot], linewidth=2)
plt.xlabel('step [h]')
plt.ylabel('P [kW]')
plt.title('Quantiles on first 300 samples')