pytorch-ts/examples/Time-Grad2-Electricity.ipynb

In [1]:

# import warnings
# warnings.simplefilter("ignore")

# autoreload import your package
%load_ext autoreload
%autoreload 2

In [2]:

%matplotlib inline
import matplotlib.pyplot as plt

import numpy as np
import pandas as pd

import torch

In [3]:

from gluonts.dataset.multivariate_grouper import MultivariateGrouper
from gluonts.dataset.repository.datasets import dataset_recipes, get_dataset
from gluonts.evaluation.backtest import make_evaluation_predictions
from gluonts.evaluation import MultivariateEvaluator

In [4]:

from pts.model.tempflow import TempFlowEstimator
from pts.model.time_grad2 import TimeGradEstimator2
from pts.model.transformer_tempflow import TransformerTempFlowEstimator
from pts import Trainer

In [5]:

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

In [6]:

def plot(target, forecast, prediction_length, prediction_intervals=(50.0, 90.0), color='g', fname=None):
    label_prefix = ""
    rows = 4
    cols = 4
    fig, axs = plt.subplots(rows, cols, figsize=(24, 24))
    axx = axs.ravel()
    seq_len, target_dim = target.shape
    
    ps = [50.0] + [
            50.0 + f * c / 2.0 for c in prediction_intervals for f in [-1.0, +1.0]
        ]
        
    percentiles_sorted = sorted(set(ps))
    
    def alpha_for_percentile(p):
        return (p / 100.0) ** 0.3
        
    for dim in range(0, min(rows * cols, target_dim)):
        ax = axx[dim]

        target[-2 * prediction_length :][dim].plot(ax=ax)
        
        ps_data = [forecast.quantile(p / 100.0)[:,dim] for p in percentiles_sorted]
        i_p50 = len(percentiles_sorted) // 2
        
        p50_data = ps_data[i_p50]
        p50_series = pd.Series(data=p50_data, index=forecast.index)
        p50_series.plot(color=color, ls="-", label=f"{label_prefix}median", ax=ax)
        
        for i in range(len(percentiles_sorted) // 2):
            ptile = percentiles_sorted[i]
            alpha = alpha_for_percentile(ptile)
            ax.fill_between(
                forecast.index,
                ps_data[i],
                ps_data[-i - 1],
                facecolor=color,
                alpha=alpha,
                interpolate=True,
            )
            # Hack to create labels for the error intervals.
            # Doesn't actually plot anything, because we only pass a single data point
            pd.Series(data=p50_data[:1], index=forecast.index[:1]).plot(
                color=color,
                alpha=alpha,
                linewidth=10,
                label=f"{label_prefix}{100 - ptile * 2}%",
                ax=ax,
            )

    legend = ["observations", "median prediction"] + [f"{k}% prediction interval" for k in prediction_intervals][::-1]    
    axx[0].legend(legend, loc="upper left")
    
    if fname is not None:
        plt.savefig(fname, bbox_inches='tight', pad_inches=0.05)

In [7]:

print(f"Available datasets: {list(dataset_recipes.keys())}")

Available datasets: ['constant', 'exchange_rate', 'solar-energy', 'electricity', 'traffic', 'exchange_rate_nips', 'electricity_nips', 'traffic_nips', 'solar_nips', 'wiki-rolling_nips', 'taxi_30min', 'kaggle_web_traffic_with_missing', 'kaggle_web_traffic_without_missing', 'kaggle_web_traffic_weekly', 'm1_yearly', 'm1_quarterly', 'm1_monthly', 'nn5_daily_with_missing', 'nn5_daily_without_missing', 'nn5_weekly', 'tourism_monthly', 'tourism_quarterly', 'tourism_yearly', 'cif_2016', 'london_smart_meters_without_missing', 'wind_farms_without_missing', 'car_parts_without_missing', 'dominick', 'fred_md', 'pedestrian_counts', 'hospital', 'covid_deaths', 'kdd_cup_2018_without_missing', 'weather', 'm3_monthly', 'm3_quarterly', 'm3_yearly', 'm3_other', 'm4_hourly', 'm4_daily', 'm4_weekly', 'm4_monthly', 'm4_quarterly', 'm4_yearly', 'm5', 'uber_tlc_daily', 'uber_tlc_hourly', 'airpassengers']

In [8]:

# exchange_rate_nips, electricity_nips, traffic_nips, solar_nips, wiki-rolling_nips, ## taxi_30min is buggy still
dataset = get_dataset("electricity_nips", regenerate=False)

In [9]:

dataset.metadata

Out [9]:

MetaData(freq='H', target=None, feat_static_cat=[CategoricalFeatureInfo(name='feat_static_cat_0', cardinality='370')], feat_static_real=[], feat_dynamic_real=[], feat_dynamic_cat=[], prediction_length=24)

In [10]:

train_grouper = MultivariateGrouper(max_target_dim=min(2000, int(dataset.metadata.feat_static_cat[0].cardinality)))

test_grouper = MultivariateGrouper(
    num_test_dates=int(len(dataset.test)/len(dataset.train)*2),
#     num_test_dates=7,
                                   max_target_dim=min(2000, int(dataset.metadata.feat_static_cat[0].cardinality)))

In [11]:

dataset_train = train_grouper(dataset.train)
dataset_test = test_grouper(dataset.test)

/home/wassname/miniforge3/envs/glounts/lib/python3.9/site-packages/gluonts/dataset/multivariate_grouper.py:191: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  return {FieldName.TARGET: np.array([funcs(data) for data in dataset])}

In [12]:

int(dataset.metadata.feat_static_cat[0].cardinality)

Out [12]:

370

In [13]:

int(dataset.metadata.feat_static_cat[0].cardinality)

Out [13]:

370

In [57]:

estimator = TimeGradEstimator2(
    target_dim=int(dataset.metadata.feat_static_cat[0].cardinality),
    prediction_length=dataset.metadata.prediction_length,
    context_length=dataset.metadata.prediction_length,
    cell_type='GRU',
    input_size=1484,
    freq=dataset.metadata.freq,
    loss_type='l2',
    scaling=False,
    diff_steps=100,
    beta_end=0.1,
    beta_schedule="linear",
    trainer=Trainer(device=device,
                    epochs=20,
                    learning_rate=1e-3,
                    num_batches_per_epoch=100,
                    batch_size=64,)
)

In [58]:

predictor = estimator.train(dataset_train, num_workers=0)
# predictor = estimator.train(dataset_train, num_workers=8)

cond_length 100

  0%|          | 0/99 [00:00<?, ?it/s]

shapes torch.Size([64, 370, 48]) torch.Size([64]) torch.Size([64, 100, 48])
cond -> cond_up torch.Size([64, 100, 48]) torch.Size([64, 370, 48])

/media/wassname/SGIronWolf/projects5/timeseries/pytorch-ts/pts/model/time_grad2/gaussian_diffusion_ou.py:283: UserWarning: Using a target size (torch.Size([64, 1, 48])) that is different to the input size (torch.Size([64, 370, 46])). This will likely lead to incorrect results due to broadcasting. Please ensure they have the same size.
  loss = F.mse_loss(x_recon, noise_rand)

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Input In [58], in <cell line: 1>()
----> 1 predictor = estimator.train(dataset_train, num_workers=0)

File /media/wassname/SGIronWolf/projects5/timeseries/pytorch-ts/pts/model/estimator.py:179, in PyTorchEstimator.train(self, training_data, validation_data, num_workers, prefetch_factor, shuffle_buffer_length, cache_data, **kwargs)
    169 def train(
    170     self,
    171     training_data: Dataset,
   (...)
    177     **kwargs,
    178 ) -> PyTorchPredictor:
--> 179     return self.train_model(
    180         training_data,
    181         validation_data,
    182         num_workers=num_workers,
    183         prefetch_factor=prefetch_factor,
    184         shuffle_buffer_length=shuffle_buffer_length,
    185         cache_data=cache_data,
    186         **kwargs,
    187     ).predictor

File /media/wassname/SGIronWolf/projects5/timeseries/pytorch-ts/pts/model/estimator.py:151, in PyTorchEstimator.train_model(self, training_data, validation_data, num_workers, prefetch_factor, shuffle_buffer_length, cache_data, **kwargs)
    133     validation_iter_dataset = TransformedIterableDataset(
    134         dataset=validation_data,
    135         transform=transformation
   (...)
    139         cache_data=cache_data,
    140     )
    141     validation_data_loader = DataLoader(
    142         validation_iter_dataset,
    143         batch_size=self.trainer.batch_size,
   (...)
    148         **kwargs,
    149     )
--> 151 self.trainer(
    152     net=trained_net,
    153     train_iter=training_data_loader,
    154     validation_iter=validation_data_loader,
    155 )
    157 return TrainOutput(
    158     transformation=transformation,
    159     trained_net=trained_net,
   (...)
    162     ),
    163 )

File /media/wassname/SGIronWolf/projects5/timeseries/pytorch-ts/pts/trainer.py:67, in Trainer.__call__(self, net, train_iter, validation_iter)
     64 optimizer.zero_grad()
     66 inputs = [v.to(self.device) for v in data_entry.values()]
---> 67 output = net(*inputs)
     69 if isinstance(output, (list, tuple)):
     70     loss = output[0]

File ~/miniforge3/envs/glounts/lib/python3.9/site-packages/torch/nn/modules/module.py:1190, in Module._call_impl(self, *input, **kwargs)
   1186 # If we don't have any hooks, we want to skip the rest of the logic in
   1187 # this function, and just call forward.
   1188 if not (self._backward_hooks or self._forward_hooks or self._forward_pre_hooks or _global_backward_hooks
   1189         or _global_forward_hooks or _global_forward_pre_hooks):
-> 1190     return forward_call(*input, **kwargs)
   1191 # Do not call functions when jit is used
   1192 full_backward_hooks, non_full_backward_hooks = [], []

File /media/wassname/SGIronWolf/projects5/timeseries/pytorch-ts/pts/model/time_grad2/time_grad_network.py:407, in TimeGradTrainingNetwork2.forward(self, target_dimension_indicator, past_time_feat, past_target_cdf, past_observed_values, past_is_pad, future_time_feat, future_target_cdf, future_observed_values)
    405 target = target.permute(0, 2, 1)
    406 distr_args = distr_args.permute(0, 2, 1)
--> 407 likelihoods = self.diffusion.log_prob(target, distr_args).unsqueeze(-1)
    409 # assert_shape(likelihoods, (-1, seq_len, 1))
    411 past_observed_values = torch.min(
    412     past_observed_values, 1 - past_is_pad.unsqueeze(-1)
    413 )

File /media/wassname/SGIronWolf/projects5/timeseries/pytorch-ts/pts/model/time_grad2/gaussian_diffusion_ou.py:298, in GaussianDiffusionOU.log_prob(self, x, cond, *args, **kwargs)
    295 B, T, _ = x.shape
    297 time = torch.randint(0, self.num_timesteps, (B,), device=x.device).long()
--> 298 loss = self.p_losses(
    299     x, cond, time, *args, **kwargs
    300 )
    302 return loss

File /media/wassname/SGIronWolf/projects5/timeseries/pytorch-ts/pts/model/time_grad2/gaussian_diffusion_ou.py:283, in GaussianDiffusionOU.p_losses(self, x_start, cond, t)
    281     loss = F.l1_loss(x_recon, noise_rand)
    282 elif self.loss_type == "l2":
--> 283     loss = F.mse_loss(x_recon, noise_rand)
    284 elif self.loss_type == "huber":
    285     loss = F.smooth_l1_loss(x_recon, noise_rand)

File ~/miniforge3/envs/glounts/lib/python3.9/site-packages/torch/nn/functional.py:3291, in mse_loss(input, target, size_average, reduce, reduction)
   3288 if size_average is not None or reduce is not None:
   3289     reduction = _Reduction.legacy_get_string(size_average, reduce)
-> 3291 expanded_input, expanded_target = torch.broadcast_tensors(input, target)
   3292 return torch._C._nn.mse_loss(expanded_input, expanded_target, _Reduction.get_enum(reduction))

File ~/miniforge3/envs/glounts/lib/python3.9/site-packages/torch/functional.py:74, in broadcast_tensors(*tensors)
     72 if has_torch_function(tensors):
     73     return handle_torch_function(broadcast_tensors, tensors, *tensors)
---> 74 return _VF.broadcast_tensors(tensors)

RuntimeError: The size of tensor a (46) must match the size of tensor b (48) at non-singleton dimension 2

In [ ]:

%debug

> /home/wassname/miniforge3/envs/glounts/lib/python3.9/site-packages/torch/functional.py(74)broadcast_tensors()
     72     if has_torch_function(tensors):
     73         return handle_torch_function(broadcast_tensors, tensors, *tensors)
---> 74     return _VF.broadcast_tensors(tensors)  # type: ignore[attr-defined]
     75 
     76 

ipdb> u
> /home/wassname/miniforge3/envs/glounts/lib/python3.9/site-packages/torch/nn/functional.py(3291)mse_loss()
   3289         reduction = _Reduction.legacy_get_string(size_average, reduce)
   3290 
-> 3291     expanded_input, expanded_target = torch.broadcast_tensors(input, target)
   3292     return torch._C._nn.mse_loss(expanded_input, expanded_target, _Reduction.get_enum(reduction))
   3293 

ipdb> u
> /media/wassname/SGIronWolf/projects5/timeseries/pytorch-ts/pts/model/time_grad2/gaussian_diffusion_ou.py(283)p_losses()
    281             loss = F.l1_loss(x_recon, noise_rand)
    282         elif self.loss_type == "l2":
--> 283             loss = F.mse_loss(x_recon, noise_rand)
    284         elif self.loss_type == "huber":
    285             loss = F.smooth_l1_loss(x_recon, noise_rand)

ipdb> u
> /media/wassname/SGIronWolf/projects5/timeseries/pytorch-ts/pts/model/time_grad2/gaussian_diffusion_ou.py(298)log_prob()
    296 
    297         time = torch.randint(0, self.num_timesteps, (B,), device=x.device).long()
--> 298         loss = self.p_losses(
    299             x, cond, time, *args, **kwargs
    300         )

ipdb> x.shape
torch.Size([64, 370, 48])
ipdb> cond.shape
torch.Size([64, 100, 48])
ipdb> time.shape
torch.Size([64])

In [ ]:

forecast_it, ts_it = make_evaluation_predictions(dataset=dataset_test,
                                                 predictor=predictor,
                                                 num_samples=100)

In [ ]:

forecasts = list(forecast_it)
targets = list(ts_it)

In [ ]:

plot(
    target=targets[0],
    forecast=forecasts[0],
    prediction_length=dataset.metadata.prediction_length,
)
plt.show()

In [ ]:

evaluator = MultivariateEvaluator(quantiles=(np.arange(20)/20.0)[1:], 
                                  target_agg_funcs={'sum': np.sum})

In [ ]:

agg_metric, item_metrics = evaluator(targets, forecasts, num_series=len(dataset_test))

In [ ]:

print("CRPS:", agg_metric["mean_wQuantileLoss"])
print("ND:", agg_metric["ND"])
print("NRMSE:", agg_metric["NRMSE"])
print("")
print("CRPS-Sum:", agg_metric["m_sum_mean_wQuantileLoss"])
print("ND-Sum:", agg_metric["m_sum_ND"])
print("NRMSE-Sum:", agg_metric["m_sum_NRMSE"])

In [ ]:

In [ ]:

In [ ]:

In [ ]:

Scratch: debug transforms

See estimator

In [ ]:

predictor = estimator.train(dataset_train, num_workers=0)

In [ ]:

%debug

In [ ]:

from tqdm.auto import tqdm

In [ ]:

transformation = estimator.create_transformation()
training_instance_splitter = estimator.create_instance_splitter("training")

In [ ]:

g1 = trans.apply(dataset_train, is_train=True)
b = next(iter(g1))
b.keys()

In [ ]:

for _ in tqdm(g1):
    pass

In [ ]:

transform = transformation + training_instance_splitter
g2 = transform.apply(dataset_train, is_train=True)
gg = iter(g2)
b = next(gg)
b.keys()

In [ ]:

for _ in tqdm(g2):
    pass

In [ ]:

b = next(gg)
b.keys()

In [ ]:

from pts.model import get_module_forward_input_names
from gluonts.transform import SelectFields, Transformation

trained_net = estimator.create_training_network(estimator.trainer.device)
input_names = get_module_forward_input_names(trained_net)
transform = transformation + training_instance_splitter + SelectFields(input_names)
g = transform.apply(dataset_train, is_train=True)
b = next(iter(g))
b.keys()

In [ ]:

In [ ]:

from pts.dataset.loader import TransformedIterableDataset
training_data = dataset_train
training_iter_dataset = TransformedIterableDataset(
    dataset=training_data,
    transform=transformation
    + training_instance_splitter
    + SelectFields(input_names),
    is_train=True,
    shuffle_buffer_length=None,
#     cache_data=cache_data,
)
training_iter_dataset

In [ ]:

next(iter(training_iter_dataset))

In [ ]:

from torch.utils.data import DataLoader
training_data_loader = DataLoader(
    training_iter_dataset,
    batch_size=estimator.trainer.batch_size,
    num_workers=0,
#     prefetch_factor=prefetch_factor,
    pin_memory=True,
    worker_init_fn=estimator._worker_init_fn,
#         **kwargs,
)

In [ ]:

next(iter(training_data_loader))

In [ ]:

for b in tqdm(training_data_loader):
    pass

In [ ]: