reproduce_experiments_additional.py

# Angus Dempster, Francois Petitjean, Geoff Webb

# Dempster A, Petitjean F, Webb GI (2019) ROCKET: Exceptionally fast and
# accurate time series classification using random convolutional kernels.
# arXiv:1910.13051

import argparse
import numpy as np
import pandas as pd
import time

from numba import njit, prange

from sklearn.linear_model import RidgeClassifierCV

from rocket_functions import generate_kernels, apply_kernels, apply_kernel

# == notes =====================================================================

# - This script is intended to allow for reproduction of the experiments on the
#   additional 2018 datasets in the UCR archive, using the txt versions of those
#   datasets from timeseriesclassification.com (Univariate2018_arff.zip).
# - This code has significant overlap with *reproduce_experiments_bakeoff.py*
#   but, for convenience, is provided as a separate script.
# - The differences from *reproduce_experiments_bakeoff.py* relate to:
#   - normalising input time series;
#   - handling missing values (missing values are interpolated); and
#   - hadling variable length time series (time series are rescaled or used
#     "as is", using a variation of *apply_kernels(...)*, as determined by
#     10-fold cross-validation).
# - The required arguments for this script are:
#   - -i or --input_path, the parent directory for the datasets; and
#   - -o or --output_path, to save "results_additional.csv".
# - Optional arguments allow you to set the number of runs, -n or --num_runs,
#   and the number of kernels, -k or --num_kernels.
# - If input_path is ".../Univariate_arff/", then each dataset should be
#   located at "{input_path}/{dataset_name}/{dataset_name}_TRAIN.txt", etc.

# == parse arguments ===========================================================

parser = argparse.ArgumentParser()

parser.add_argument("-i", "--input_path", required = True)
parser.add_argument("-o", "--output_path", required = True)
parser.add_argument("-n", "--num_runs", type = int, default = 10)
parser.add_argument("-k", "--num_kernels", type = int, default = 10_000)

arguments = parser.parse_args()

# == additional dataset names ==================================================

dataset_names_additional = \
(
    "ACSF1",
    "AllGestureWiimoteX",
    "AllGestureWiimoteY",
    "AllGestureWiimoteZ",
    "BME",
    "Chinatown",
    "Crop",
    "DodgerLoopDay",
    "DodgerLoopGame",
    "DodgerLoopWeekend",
    "EOGHorizontalSignal",
    "EOGVerticalSignal",
    "EthanolLevel",
    "FreezerRegularTrain",
    "FreezerSmallTrain",
    "Fungi",
    "GestureMidAirD1",
    "GestureMidAirD2",
    "GestureMidAirD3",
    "GesturePebbleZ1",
    "GesturePebbleZ2",
    "GunPointAgeSpan",
    "GunPointMaleVersusFemale",
    "GunPointOldVersusYoung",
    "HouseTwenty",
    "InsectEPGRegularTrain",
    "InsectEPGSmallTrain",
    "MelbournePedestrian",
    "MixedShapesRegularTrain",
    "MixedShapesSmallTrain",
    "PLAID",
    "PickupGestureWiimoteZ",
    "PigAirwayPressure",
    "PigArtPressure",
    "PigCVP",
    "PowerCons",
    "Rock",
    "SemgHandGenderCh2",
    "SemgHandMovementCh2",
    "SemgHandSubjectCh2",
    "ShakeGestureWiimoteZ",
    "SmoothSubspace",
    "UMD"
)

# == apply kernels, variable input lengths =====================================

# if variable length time series are used "as is", the effective length of some
# kernels, including dilation (if set with reference to, e.g., the longest
# time series in a dataset, the default behaviour), may be larger than some
# input time series; this is irrelevant if padding is applied; even where
# padding is not applied, this should only affect a minority of kernels
# (dilation for most kernels is relatively small); the default behaviour of
# *apply_kernels_jagged(...)* is to "skip" incompatible kernels (i.e., where the
# effective size of the kernel including dilation is larger than the input time
# series including padding)

@njit(parallel = True, fastmath = True)
def apply_kernels_jagged(X, kernels, input_lengths):

    weights, lengths, biases, dilations, paddings = kernels

    num_examples = len(X)
    num_kernels = len(weights)

    # initialise output
    _X = np.zeros((num_examples, num_kernels * 2)) # 2 features per kernel

    for i in prange(num_examples):

        for j in range(num_kernels):

            # skip incompatible kernels (effective length is "too big" without padding)
            if (input_lengths[i] + (2 * paddings[j])) > ((lengths[j] - 1) * dilations[j]):

                _X[i, (j * 2):((j * 2) + 2)] = \
                apply_kernel(X[i][:input_lengths[i]], weights[j][:lengths[j]], lengths[j], biases[j], dilations[j], paddings[j])

    return _X

# == additional convenience function ===========================================

def run_additional(training_data, test_data, num_runs = 10, num_kernels = 10_000):

    # assumes variable length time series are padded with nan
    get_input_lengths = lambda X : X.shape[1] - (~np.isnan(np.flip(X, 1))).argmax(1)

    def rescale(X, reference_length):
        _X = np.zeros([len(X), reference_length])
        input_lengths = get_input_lengths(X)
        for i in range(len(X)):
            _X[i] = np.interp(np.linspace(0, 1, reference_length), np.linspace(0, 1, input_lengths[i]), X[i][:input_lengths[i]])
        return _X

    def interpolate_nan(X):
        _X = X.copy()
        good = ~np.isnan(X)
        for i in np.where(np.any(~good, 1))[0]:
            _X[i] = np.interp(np.arange(len(X[i])), np.where(good[i])[0], X[i][good[i]])
        return _X

    results = np.zeros(num_runs)
    timings = np.zeros([4, num_runs]) # training transform, test transform, training, test

    Y_training, X_training = training_data[:, 0].astype(np.int), training_data[:, 1:]
    Y_test, X_test = test_data[:, 0].astype(np.int), test_data[:, 1:]

    variable_lengths = False

    # handle three cases: (1) same lengths, no missing values; (2) same lengths,
    # missing values; and (3) variable lengths, no missing values

    if np.any(np.isnan(X_training)):

        input_lengths_training = get_input_lengths(X_training)
        input_lengths_training_max = input_lengths_training.max()
        input_lengths_test = get_input_lengths(X_test)

        # missing values (same lengths)
        if np.all(input_lengths_training == input_lengths_training_max):

            X_training = interpolate_nan(X_training)
            X_test = interpolate_nan(X_test)

        # variable lengths (no missing values)
        else:

            variable_lengths = True
            num_folds = 10
            cross_validation_results = np.zeros([2, num_folds])

    # normalise time series
    X_training = (X_training - np.nanmean(X_training, axis = 1, keepdims = True)) / (np.nanstd(X_training, axis = 1, keepdims = True) + 1e-8)
    X_test = (X_test - np.nanmean(X_test, axis = 1, keepdims = True)) / (np.nanstd(X_test, axis = 1, keepdims = True) + 1e-8)

    for i in range(num_runs):

        # -- variable lengths --------------------------------------------------

        if variable_lengths:

            kernels = generate_kernels(input_lengths_training_max, num_kernels)

            time_a = time.perf_counter()
            X_training_transform_rescale = apply_kernels(rescale(X_training, input_lengths_training_max), kernels)
            X_training_transform_jagged = apply_kernels_jagged(X_training, kernels, input_lengths_training)
            time_b = time.perf_counter()
            timings[0, i] = time_b - time_a

            # indices for cross-validation folds
            I = np.random.permutation(len(X_training))
            I = np.array_split(I, num_folds)

            time_a = time.perf_counter()

            # j = 0 -> rescale
            # j = 1 -> "as is" ("jagged")
            for j in range(2):

                for k in range(num_folds):

                    VA, *TR = np.roll(I, k, axis = 0)
                    TR = np.concatenate(TR)

                    classifier = RidgeClassifierCV(alphas = 10 ** np.linspace(-3, 3, 10), normalize = True)

                    if j == 0: # rescale

                        classifier.fit(X_training_transform_rescale[TR], Y_training[TR])
                        cross_validation_results[j][k] = classifier.score(X_training_transform_rescale[VA], Y_training[VA])

                    elif j == 1: # jagged

                        classifier.fit(X_training_transform_jagged[TR], Y_training[TR])
                        cross_validation_results[j][k] = classifier.score(X_training_transform_jagged[VA], Y_training[VA])

            best = cross_validation_results.sum(1).argmax()
            time_b = time.perf_counter()
            timings[2, i] = time_b - time_a

            classifier = RidgeClassifierCV(alphas = 10 ** np.linspace(-3, 3, 10), normalize = True)

            if best == 0: # rescale

                time_a = time.perf_counter()
                X_test_transform_rescale = apply_kernels(rescale(X_test, input_lengths_training_max), kernels)
                time_b = time.perf_counter()
                timings[1, i] = time_b - time_a

                time_a = time.perf_counter()
                classifier.fit(X_training_transform_rescale, Y_training)
                time_b = time.perf_counter()
                timings[2, i] += time_b - time_a

                time_a = time.perf_counter()
                results[i] = classifier.score(X_test_transform_rescale, Y_test)
                time_b = time.perf_counter()
                timings[3, i] = time_b - time_a

            elif best == 1: # jagged

                time_a = time.perf_counter()
                X_test_transform_jagged = apply_kernels_jagged(X_test, kernels, input_lengths_test)
                time_b = time.perf_counter()
                timings[1, i] = time_b - time_a

                time_a = time.perf_counter()
                classifier.fit(X_training_transform_jagged, Y_training)
                time_b = time.perf_counter()
                timings[2, i] += time_b - time_a

                time_a = time.perf_counter()
                results[i] = classifier.score(X_test_transform_jagged, Y_test)
                time_b = time.perf_counter()
                timings[3, i] = time_b - time_a

        # -- same lengths ------------------------------------------------------

        else:

            kernels = generate_kernels(X_training.shape[1], num_kernels)

            # -- transform training --------------------------------------------

            time_a = time.perf_counter()
            X_training_transform = apply_kernels(X_training, kernels)
            time_b = time.perf_counter()
            timings[0, i] = time_b - time_a

            # -- transform test ------------------------------------------------

            time_a = time.perf_counter()
            X_test_transform = apply_kernels(X_test, kernels)
            time_b = time.perf_counter()
            timings[1, i] = time_b - time_a

            # -- training ------------------------------------------------------

            time_a = time.perf_counter()
            classifier = RidgeClassifierCV(alphas = 10 ** np.linspace(-3, 3, 10), normalize = True)
            classifier.fit(X_training_transform, Y_training)
            time_b = time.perf_counter()
            timings[2, i] = time_b - time_a

            # -- test ----------------------------------------------------------

            time_a = time.perf_counter()
            results[i] = classifier.score(X_test_transform, Y_test)
            time_b = time.perf_counter()
            timings[3, i] = time_b - time_a

    return results, timings

# == run through the additional datasets =======================================

results_additional = pd.DataFrame(index = dataset_names_additional,
                                  columns = ["accuracy_mean",
                                             "accuracy_standard_deviation",
                                             "time_training_seconds",
                                             "time_test_seconds"],
                               data = 0)
results_additional.index.name = "dataset"

compiled = False

print(f"RUNNING".center(80, "="))

for dataset_name in dataset_names_additional:

    print(f"{dataset_name}".center(80, "-"))

    # -- read data -------------------------------------------------------------

    print(f"Loading data".ljust(80 - 5, "."), end = "", flush = True)

    if dataset_name != "PLAID":

        training_data = np.loadtxt(f"{arguments.input_path}/{dataset_name}/{dataset_name}_TRAIN.txt")
        test_data = np.loadtxt(f"{arguments.input_path}/{dataset_name}/{dataset_name}_TEST.txt")

    else:

        training_data = np.loadtxt(f"{arguments.input_path}/{dataset_name}/{dataset_name}_TRAIN.txt", delimiter = ",")
        test_data = np.loadtxt(f"{arguments.input_path}/{dataset_name}/{dataset_name}_TEST.txt", delimiter = ",")

    print("Done.")

    # -- precompile ------------------------------------------------------------

    if not compiled:

        print(f"Compiling ROCKET functions (once only)".ljust(80 - 5, "."), end = "", flush = True)

        _ = generate_kernels(100, 10)
        apply_kernels(np.zeros_like(training_data)[:, 1:], _)
        apply_kernels_jagged(np.zeros_like(training_data)[:, 1:], _, np.array([training_data.shape[1]] * len(training_data)))
        compiled = True

        print("Done.")

    # -- run -------------------------------------------------------------------

    print(f"Performing runs".ljust(80 - 5, "."), end = "", flush = True)

    results, timings = run_additional(training_data, test_data,
                                      num_runs = arguments.num_runs,
                                      num_kernels = arguments.num_kernels)
    timings_mean = timings.mean(1)

    print("Done.")

    # -- store results ---------------------------------------------------------

    results_additional.loc[dataset_name, "accuracy_mean"] = results.mean()
    results_additional.loc[dataset_name, "accuracy_standard_deviation"] = results.std()
    results_additional.loc[dataset_name, "time_training_seconds"] = timings_mean[[0, 2]].sum()
    results_additional.loc[dataset_name, "time_test_seconds"] = timings_mean[[1, 3]].sum()

print(f"FINISHED".center(80, "="))

results_additional.to_csv(f"{arguments.output_path}/results_additional.csv")