<jupyter_start><jupyter_text>created by <NAME> - email: <EMAIL>  ROCKET (RandOm Convolutional KErnel Transform) is a new Time Series Classification (TSC) method that has just been released (Oct 29th, 2019), and has achieved **state-of-the-art performance on the UCR univariate time series classification datasets, surpassing HIVE-COTE (the previous state of the art since 2017) in accuracy, with exceptional speed compared to other traditional DL methods.** To achieve these 2 things at once is **VERY IMPRESSIVE**. ROCKET is certainly a new TSC method you should try.Authors:<NAME>., <NAME>., & <NAME>. (2019). ROCKET: Exceptionally fast and accurate time series classification using random convolutional kernels. arXiv preprint arXiv:1910.13051.[paper](https://arxiv.org/pdf/1910.13051)There are 2 main limitations to the original ROCKET method though:- Released code doesn't handle multivariate data- It doesn't run on a GPU, so it's slow when used with a large datasetsIn this notebook you will learn: - how you can use the original ROCKET method- you will also learn about a new ROCKET version I have developed in Pytorch, that handles both **univariate and multivariate** data, and uses **GPU**- you will see how you can integrate the ROCKET features with fastai or other classifiers  Import libraries 📚<jupyter_code>import sys
import os
ISCOLAB = 'google.colab' in sys.modules
if ISCOLAB:
    # for bleeding edge
    !pip install git+https://github.com/fastai/fastcore.git@master -q
    !pip install git+https://github.com/fastai/fastai.git@master -q
    !pip install git+https://github.com/timeseriesAI/timeseriesAI.git@master -q
    
    # for latest stable version
    # !pip install tsai -q
    
import tsai
from tsai.all import *
display(HTML("<style>.container {width:95% !important; }</style>"))
print('tsai       :', tsai.__version__)
print('fastai    :', fastai.__version__)
print('fastcore   :', fastcore.__version__)
print('torch      :', torch.__version__)
print('scipy      :', sp.__version__)
print('numpy      :', np.__version__)
print('pandas     :', pd.__version__)
print(f'Total RAM  : {bytes2GB(psutil.virtual_memory().total):5.2f} GB')
print(f'Used RAM   : {bytes2GB(psutil.virtual_memory().used):5.2f} GB')
print('n_cpus     :', cpus)
iscuda = torch.cuda.is_available()
if iscuda: print('device     : {} ({})'.format(device, torch.cuda.get_device_name(0)))
else: print('device     :', device)<jupyter_output>tsai       : 0.1.0
fastai    : 2.0.8
fastcore   : 1.0.1
torch      : 1.6.0
scipy      : 1.5.2
numpy      : 1.19.1
pandas     : 1.1.1
Total RAM  : 31.26 GB
Used RAM   : 14.42 GB
n_cpus     : 8
device     : cuda (GeForce RTX 2080 Ti)<jupyter_text>How to use the original ROCKET method? 🚀 ROCKET is applied in 2 phases:1. Generate features from each time series: ROCKET calculates 20k features from each time series, independently of the sequence length. 2. Apply a classifier to those calculated features. Those features are then used by the classifier of your choice. In the original code they use 2 simple linear classifiers: RidgeClassifierCV and Logistic Regression, but you can use any classifier.  1️⃣ Generate featuresLet's first generate the features. We'll import data from a UCR Time Series dataset.The original method requires the time series to be in a 2d array of shape (samples, len). Remember than only univariate sequences are allow in this original method.<jupyter_code>X_train, y_train, X_valid, y_valid = get_UCR_data('OliveOil')
seq_len = X_train.shape[-1]
X_train = X_train[:, 0]
X_valid = X_valid[:, 0]
labels = np.unique(y_train)
transform = {}
for i, l in enumerate(labels): transform[l] = i
y_train = np.vectorize(transform.get)(y_train)
y_valid = np.vectorize(transform.get)(y_valid)<jupyter_output><empty_output><jupyter_text>Now we normalize the data to mean 0 and std 1 **'per sample'** (recommended by the authors), that is they set each sample to mean 0 and std 1).<jupyter_code>X_train = (X_train - X_train.mean(axis = 1, keepdims = True)) / (X_train.std(axis = 1, keepdims = True) + 1e-8)
X_valid = (X_valid - X_valid.mean(axis = 1, keepdims = True)) / (X_valid.std(axis = 1, keepdims = True) + 1e-8)
X_train.mean(axis = 1, keepdims = True).shape<jupyter_output><empty_output><jupyter_text>To generate the features, we first need to create the 10k random kernels that will be used to process the data.<jupyter_code>kernels = generate_kernels(seq_len, 10000)<jupyter_output><empty_output><jupyter_text>Now we apply those ramdom kernels to the data<jupyter_code>X_train_tfm = apply_kernels(X_train, kernels)
X_valid_tfm = apply_kernels(X_valid, kernels)<jupyter_output><empty_output><jupyter_text>2️⃣ Apply a classifierSo now we have the features, and we are ready to apply a classifier. Let's use a simple, linear RidgeClassifierCV as they propose in the paper. We first instantiate it. Note:alphas: Array of alpha values to try. Regularization strength; must be a positive float. Regularization improves the conditioning of the problem and reduces the variance of the estimates. Larger values specify stronger regularization. Alpha corresponds to C^-1 in other linear models such as LogisticRegression or LinearSVC.<jupyter_code>from sklearn.linear_model import RidgeClassifierCV
classifier = RidgeClassifierCV(alphas=np.logspace(-3, 3, 7), normalize=True)
 classifier.fit(X_train_tfm, y_train)
classifier.score(X_valid_tfm, y_valid)<jupyter_output><empty_output><jupyter_text>☣️ **This is pretty impressive! It matches or exceeds the state-of-the-art performance without any fine tuning in <2 seconds!!!**<jupyter_code>kernels = generate_kernels(seq_len, 10000)
X_train_tfm = apply_kernels(X_train, kernels)
X_valid_tfm = apply_kernels(X_valid, kernels)
from sklearn.linear_model import RidgeClassifierCV
classifier = RidgeClassifierCV(alphas=np.logspace(-3, 3, 7), normalize=True)
classifier.fit(X_train_tfm, y_train)
classifier.score(X_valid_tfm, y_valid)<jupyter_output><empty_output><jupyter_text>⚠️ Bear in mind that this process is not deterministic since there is randomness involved in the kernels. In thiis case, performance may vary between .9 to .933  How to use ROCKET with large and/ or multivariate datasets on GPU? ❓ As stated before, the current ROCKET method doesn't support multivariate time series or GPU. This may be a drawback in some cases. To overcome both limitations I've created a multivariate ROCKET on GPU in Pytorch.  The way to use ROCKET in Pytorch is the following:  1️⃣ Generate featuresFirst you prepare the input data and normalize it per sample. The input to ROCKET Pytorch is a 3d tensor of shape (samples, vars, len), preferrable on gpu.<jupyter_code>X_train, y_train, X_valid, y_valid = get_UCR_data('HandMovementDirection')
_, features, seq_len = X_train.shape
X_train = (X_train - X_train.mean(axis=(1, 2), keepdims=True)) / (
    X_train.std(axis=(1, 2), keepdims=True) + 1e-8)
X_valid = (X_valid - X_valid.mean(axis=(1, 2), keepdims=True)) / (
    X_valid.std(axis=(1, 2), keepdims=True) + 1e-8)
X_train = torch.tensor(X_train, dtype=torch.float32, device=device)
X_valid = torch.tensor(X_valid, dtype=torch.float32, device=device)
labels = np.unique(y_train)
transform = {}
for i, l in enumerate(labels): transform[l] = i
y_train = np.vectorize(transform.get)(y_train)
y_valid = np.vectorize(transform.get)(y_valid)
print(X_train.shape, X_valid.shape, y_train.shape, y_valid.shape)<jupyter_output>torch.Size([160, 10, 400]) torch.Size([74, 10, 400]) (160,) (74,)<jupyter_text>Then you instantiate the Pytorch model.<jupyter_code>n_kernels=10_000
kss=[7, 9, 11]
model = ROCKET(features, seq_len, n_kernels=n_kernels, kss=kss).to(device)<jupyter_output><empty_output><jupyter_text>And we now transform the original data, creating 20k features per sample<jupyter_code>X_train_tfm = model(X_train).cpu().numpy()
X_valid_tfm = model(X_valid).cpu().numpy()<jupyter_output><empty_output><jupyter_text>2️⃣ Apply a classifier  RidgeClassifierCV And now you apply a classifier of your choice. With RidgeClassifierCV in particular, there's no need to normalize the calculated features before passing them to the classifier, as it does it internally (if normalize is set to True as recommended by the authors).<jupyter_code>from sklearn.linear_model import RidgeClassifierCV
ridge = RidgeClassifierCV(alphas=np.logspace(-8, 8, 17), normalize=True)
ridge.fit(X_train_tfm, y_train)
print('alpha: {:.2E}  train: {:.5f}  valid: {:.5f}'.format(ridge.alpha_, 
                                                           ridge.score(X_train_tfm, y_train), 
                                                           ridge.score(X_valid_tfm, y_valid)))<jupyter_output>alpha: 1.00E+01  train: 1.00000  valid: 0.48649<jupyter_text>Logistic Regression In the case of other classifiers (like Logistic Regression), the authors recommend a per-feature normalization.<jupyter_code>eps = 1e-6
Cs = np.logspace(-5, 5, 11)
from sklearn.linear_model import LogisticRegression
best_loss = np.inf
for i, C in enumerate(Cs):
    f_mean = X_train_tfm.mean(axis=0, keepdims=True)
    f_std = X_train_tfm.std(axis=0, keepdims=True) + eps  # epsilon to avoid dividing by 0
    X_train_tfm2 = (X_train_tfm - f_mean) / f_std
    X_valid_tfm2 = (X_valid_tfm - f_mean) / f_std
    classifier = LogisticRegression(penalty='l2', C=C, n_jobs=-1)
    classifier.fit(X_train_tfm2, y_train)
    probas = classifier.predict_proba(X_train_tfm2)
    loss = nn.CrossEntropyLoss()(torch.tensor(probas), torch.tensor(y_train)).item()
    train_score = classifier.score(X_train_tfm2, y_train)
    val_score = classifier.score(X_valid_tfm2, y_valid)
    if loss < best_loss:
        best_eps = eps
        best_C = C
        best_loss = loss
        best_train_score = train_score
        best_val_score = val_score
    print('{:2} eps: {:.2E}  C: {:.2E}  loss: {:.5f}  train_acc: {:.5f}  valid_acc: {:.5f}'.format(
        i, eps, C, loss, train_score, val_score))
print('\nBest result:')
print('eps: {:.2E}  C: {:.2E}  train_loss: {:.5f}  train_acc: {:.5f}  valid_acc: {:.5f}'.format(
        best_eps, best_C, best_loss, best_train_score, best_val_score))<jupyter_output>0 eps: 1.00E-06  C: 1.00E-05  loss: 1.35156  train_acc: 0.82500  valid_acc: 0.40541
 1 eps: 1.00E-06  C: 1.00E-04  loss: 1.15439  train_acc: 1.00000  valid_acc: 0.44595
 2 eps: 1.00E-06  C: 1.00E-03  loss: 0.85366  train_acc: 1.00000  valid_acc: 0.50000
 3 eps: 1.00E-06  C: 1.00E-02  loss: 0.76184  train_acc: 1.00000  valid_acc: 0.52703
 4 eps: 1.00E-06  C: 1.00E-01  loss: 0.74626  train_acc: 1.00000  valid_acc: 0.50000
 5 eps: 1.00E-06  C: 1.00E+00  loss: 0.74401  train_acc: 1.00000  valid_acc: 0.48649
 6 eps: 1.00E-06  C: 1.00E+01  loss: 0.74371  train_acc: 1.00000  valid_acc: 0.48649
 7 eps: 1.00E-06  C: 1.00E+02  loss: 0.74367  train_acc: 1.00000  valid_acc: 0.48649
 8 eps: 1.00E-06  C: 1.00E+03  loss: 0.74367  train_acc: 1.00000  valid_acc: 0.47297
 9 eps: 1.00E-06  C: 1.00E+04  loss: 0.74367  train_acc: 1.00000  valid_acc: 0.47297
10 eps: 1.00E-06  C: 1.00E+05  loss: 0.74367  train_acc: 1.00000  valid_acc: 0.47297

Best result:
eps: 1.00E-06  C: 1.00E+05  train_loss: 0.74367  tr[...]<jupyter_text>This result is amazing!! The previous state of the art (Inceptiontime) was .37837☣️ Note: Epsilon has a large impact on the result. You can actually test several values to find the one that best fits your problem, but bear in mind you can only select C and epsilon based on train data!!!   RandomSearch One way to do this would be to performa a random search using several epsilon and C values<jupyter_code>n_tests = 10
epss = np.logspace(-8, 0, 9)
Cs = np.logspace(-5, 5, 11)

from sklearn.linear_model import LogisticRegression
best_loss = np.inf
for i in range(n_tests):
    eps = np.random.choice(epss)
    C = np.random.choice(Cs)
    f_mean = X_train_tfm.mean(axis=0, keepdims=True)
    f_std = X_train_tfm.std(axis=0, keepdims=True) + eps  # epsilon
    X_train_tfm2 = (X_train_tfm - f_mean) / f_std
    X_valid_tfm2 = (X_valid_tfm - f_mean) / f_std
    classifier = LogisticRegression(penalty='l2', C=C, n_jobs=-1)
    classifier.fit(X_train_tfm2, y_train)
    probas = classifier.predict_proba(X_train_tfm2)
    loss = nn.CrossEntropyLoss()(torch.tensor(probas), torch.tensor(y_train)).item()
    train_score = classifier.score(X_train_tfm2, y_train)
    val_score = classifier.score(X_valid_tfm2, y_valid)
    if loss < best_loss:
        best_eps = eps
        best_C = C
        best_loss = loss
        best_train_score = train_score
        best_val_score = val_score
    print('{:2}  eps: {:.2E}  C: {:.2E}  loss: {:.5f}  train_acc: {:.5f}  valid_acc: {:.5f}'.format(
        i, eps, C, loss, train_score, val_score))
print('\nBest result:')
print('eps: {:.2E}  C: {:.2E}  train_loss: {:.5f}  train_acc: {:.5f}  valid_acc: {:.5f}'.format(
        best_eps, best_C, best_loss, best_train_score, best_val_score))<jupyter_output>0  eps: 1.00E+00  C: 1.00E+02  loss: 0.74368  train_acc: 1.00000  valid_acc: 0.41892
 1  eps: 1.00E-06  C: 1.00E+01  loss: 0.74371  train_acc: 1.00000  valid_acc: 0.48649
 2  eps: 1.00E-05  C: 1.00E-04  loss: 1.15441  train_acc: 1.00000  valid_acc: 0.44595
 3  eps: 1.00E-04  C: 1.00E+00  loss: 0.74401  train_acc: 1.00000  valid_acc: 0.48649
 4  eps: 1.00E-05  C: 1.00E+00  loss: 0.74401  train_acc: 1.00000  valid_acc: 0.48649
 5  eps: 1.00E-02  C: 1.00E+04  loss: 0.74367  train_acc: 1.00000  valid_acc: 0.47297
 6  eps: 1.00E-01  C: 1.00E+01  loss: 0.74373  train_acc: 1.00000  valid_acc: 0.43243
 7  eps: 1.00E-04  C: 1.00E+01  loss: 0.74371  train_acc: 1.00000  valid_acc: 0.48649
 8  eps: 1.00E-03  C: 1.00E+05  loss: 0.74367  train_acc: 1.00000  valid_acc: 0.47297
 9  eps: 1.00E-08  C: 1.00E+05  loss: 0.74367  train_acc: 1.00000  valid_acc: 0.47297

Best result:
eps: 1.00E-08  C: 1.00E+05  train_loss: 0.74367  train_acc: 1.00000  valid_acc: 0.47297<jupyter_text>In general, the original method may be a bit faster than the GPU method, but for larger datasets, there's a great benefit in using the GPU version. In addition to this, I have also run the code on the TSC UCR multivariate datasets (all the ones that don't contain nan values), and the results are also very good, beating the previous state-of-the-art in this category as well by a large margin. For example, ROCKET reduces InceptionTime errors by 26% on average.  And how to use it with other classifiers?  Once you build the 20k features per sample, you can use them to train any classifier of your choice. If you want to use fastai, for example, with a logistic regression, you could should then build a databunch and a logistic regression model.  Integration with fastai These are the steps to use ROCKET with a classifier in fastai v2:1. Normalize the data 'per sample' (if not previously normalized)2. Calculate features3. Normalize calculated features 'per feature' (you'll get 20k means and stds)4. Create the datasets and dataloaders passing the normalized calculated features as 3d arrays.<jupyter_code>dsid = 'OliveOil'
bs = 30
eps = 1e-6

# extract data
X_train, y_train, X_valid, y_valid = get_UCR_data(dsid)

# normalize data 'per sample'
X_train = (X_train - X_train.mean(axis=(1, 2), keepdims=True)) / (
    X_train.std(axis=(1, 2), keepdims=True) + eps)
X_valid = (X_valid - X_valid.mean(axis=(1, 2), keepdims=True)) / (
    X_valid.std(axis=(1, 2), keepdims=True) + eps)

# calculate 20k features
_, features, seq_len = X_train.shape
model = ROCKET(features, seq_len, n_kernels=10000, kss=[7, 9, 11]).to(device)
X_train_tfm = model(torch.tensor(X_train, device=device)).unsqueeze(1).cpu().numpy()
X_valid_tfm = model(torch.tensor(X_valid, device=device)).unsqueeze(1).cpu().numpy()

# normalize 'per feature'
f_mean = X_train_tfm.mean(axis=0, keepdims=True)
f_std = X_train_tfm.std(axis=0, keepdims=True) + eps
X_train_tfm_norm = (X_train_tfm - f_mean) / f_std
X_valid_tfm_norm = (X_valid_tfm - f_mean) / f_std

# To use fastai v2 we concatenate train and valid, and calculate the splits:
X = np.concatenate((X_train_tfm_norm, X_valid_tfm_norm))
y = np.concatenate((y_train, y_valid))
splits = (L(list(np.arange(len(X_train)))), L(list(np.arange(len(X_train), len(X)))))
print(X.shape, y.shape)<jupyter_output>(60, 1, 20000) (60,)<jupyter_text>We now need ot create datasets and dataloaders in fastai v2:<jupyter_code>tfms  = [None, [Categorize()]]
dsets = TSDatasets(X, y, tfms=tfms, splits=splits, inplace=True)
dls   = TSDataLoaders.from_dsets(dsets.train, dsets.valid, bs=30, num_workers=0)<jupyter_output><empty_output><jupyter_text>Let's take a look at the output of ROCKET.<jupyter_code>dls.show_batch()<jupyter_output><empty_output><jupyter_text>In this case we'll use a very simple: a Logistic Regression with 20k input features and 2 classes in this case.<jupyter_code>def init(layer):
    if isinstance(layer, nn.Linear):
        nn.init.constant_(layer.weight.data, 0.)
        nn.init.constant_(layer.bias.data, 0.)
model = nn.Sequential(*[nn.Linear(20_000, dls.c), Squeeze(1)])
model.apply(init)
learn = Learner(dls, model, metrics=accuracy)
learn.save('stage-0')
learn.lr_find()
learn.load('stage-0')
learn.fit_one_cycle(10, max_lr=3e-3, wd=1e2)
learn.recorder.plot_loss()<jupyter_output><empty_output><jupyter_text>This is an excellent result we have achieved with a very simple classifier (Logistic Regression)  Boosted trees: XGBoost<jupyter_code>dsid = 'OliveOil'
eps = 1e-6

# extract data
X_train, y_train, X_valid, y_valid = get_UCR_data(dsid)

# normalize data 'per sample'
X_train = (X_train - X_train.mean(axis=(1, 2), keepdims=True)) / (
    X_train.std(axis=(1, 2), keepdims=True) + eps)
X_valid = (X_valid - X_valid.mean(axis=(1, 2), keepdims=True)) / (
    X_valid.std(axis=(1, 2), keepdims=True) + eps)

# calculate 20k features
_, features, seq_len = X_train.shape
model = ROCKET(features, seq_len, n_kernels=10000, kss=[7, 9, 11]).to(device)
X_train_tfm = model(torch.tensor(X_train, device=device)).cpu().numpy()
X_valid_tfm = model(torch.tensor(X_valid, device=device)).cpu().numpy()

# normalize 'per feature'
f_mean = X_train_tfm.mean(axis=0, keepdims=True)
f_std = X_train_tfm.std(axis=0, keepdims=True) + eps
X_train_tfm_norm = (X_train_tfm - f_mean) / f_std
X_valid_tfm_norm = (X_valid_tfm - f_mean) / f_std
import xgboost as xgb
classifier = xgb.XGBClassifier(max_depth=3,
                               learning_rate=0.1,
                               n_estimators=100,
                               verbosity=1,
                               objective='binary:logistic',
                               booster='gbtree',
                               tree_method='auto',
                               n_jobs=1,
                               gpu_id=default_device().index,
                               gamma=0,
                               min_child_weight=1,
                               max_delta_step=0,
                               subsample=.5,
                               colsample_bytree=1,
                               colsample_bylevel=1,
                               colsample_bynode=1,
                               reg_alpha=0,
                               reg_lambda=1,
                               scale_pos_weight=1,
                               base_score=0.5,
                               random_state=0,
                               missing=None)

classifier.fit(X_train_tfm_norm, y_train)
preds = classifier.predict(X_valid_tfm_norm)
(preds == y_valid).mean()<jupyter_output>[16:23:00] WARNING: /home/conda/feedstock_root/build_artifacts/xgboost_1598185621802/work/src/learner.cc:516: 
Parameters: { scale_pos_weight } might not be used.

  This may not be accurate due to some parameters are only used in language bindings but
  passed down to XGBoost core.  Or some parameters are not used but slip through this
  verification. Please open an issue if you find above cases.
{'fraction_non_alphanumeric': 0.11075837742504409, 'fraction_numerical': 0.05326278659611993, 'mean_word_length': 3.877365446055542, 'pattern_counts': {'":': 0, '<': 76, '<?xml version=': 0, '>': 73, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '139717', 'license': 'Apache-2.0', 'n_tokens_mistral': 7845, 'n_tokens_neox': 7060, 'n_words': 1981, 'path': 'tutorial_nbs/02_ROCKET_a_new_SOTA_classifier.ipynb', 'repo_name': 'williamsdoug/timeseriesAI', 'size': '19877'}