Repository: alan-turing-institute/sktime
Pull request #3352 (2022-08-26 16:14:37)

[BUG] fixes incorrect warning condition in `InvertTransform`
fixes incorrect warning condition in `InvertTransform`, which was flipped (superfluous `not`)


Patch:
diff --git a/examples/AA_datatypes_and_datasets.ipynb b/examples/AA_datatypes_and_datasets.ipynb
index 4353e6f5c..a7bdceb36 100644
--- a/examples/AA_datatypes_and_datasets.ipynb
+++ b/examples/AA_datatypes_and_datasets.ipynb
@@ -6,7 +6,7 @@
    "source": [
     "**Abstract:** this notebook give an introduction to `sktime` in-memory data containers and data sets, with associated functionality such as in-memory format validation, conversion, and data set loading.\n",
     "\n",
-    "**Set-up instructions:** on binder, this nootebook should run out-of-the-box.\n",
+    "**Set-up instructions:** on binder, this notebook should run out-of-the-box.\n",
     "\n",
     "To run this notebook as intended, ensure that `sktime` with basic dependency requirements is installed in your python environment.\n",
     "\n",
@@ -295,7 +295,7 @@
     "\n",
     "In the `\"numpy3D\"` mtype, time series panels are represented by an in-memory container `obj: np.ndarray` as follows.\n",
     "\n",
-    "* structure convention: `obj` must be 3D, i.e., `obj.shape` must have length 2.\n",
+    "* structure convention: `obj` must be 3D, i.e., `obj.shape` must have length 3.\n",
     "* instances: instances correspond to axis 0 elements of `obj`.\n",
     "* instance index: the instance index is implicit and by-convention. The `i`-th element of axis 0 (for an integer `i`) is interpreted as indicative of observing instance `i`. \n",
     "* variables: variables correspond to axis 1 elements of `obj`.\n",
diff --git a/examples/hidalgo_segmentation.ipynb b/examples/hidalgo_segmentation.ipynb
new file mode 100644
index 000000000..c7614d712
--- /dev/null
+++ b/examples/hidalgo_segmentation.ipynb
@@ -0,0 +1,66 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 29,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
+       "       1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
+       "       1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,\n",
+       "       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
+       "       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=int64)"
+      ]
+     },
+     "execution_count": 29,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "from sktime.transformations.series.hidalgo import Hidalgo\n",
+    "\n",
+    "np.random.seed(123)\n",
+    "X = np.random.rand(100, 3)\n",
+    "X[:60, 1:] += 10\n",
+    "X[60:, 1:] = 0\n",
+    "model = Hidalgo(K=2, burn_in=0.8, n_iter=1000, seed=10)\n",
+    "fitted_model = model.fit(X)\n",
+    "Z = fitted_model.transform(X)\n",
+    "Z"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3.9.13 ('sktime-dev')",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.12"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "3672c9971f48dc8b5371767a83a5c009b8758f695819d57d51811323dc5c559e"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/pyproject.toml b/pyproject.toml
index 31fe452b4..6b9509db6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -72,9 +72,11 @@ all_extras = [
     "tensorflow",
     "tsfresh>=0.17.0; python_version < '3.10'",
     "tslearn>=0.5.2",
+    "xarray",
 ]
 
 dev = [
+    "backoff",
     "httpx",
     "pre-commit",
     "pytest",
diff --git a/sktime/annotation/ggs.py b/sktime/annotation/ggs.py
new file mode 100644
index 000000000..928bfd14d
--- /dev/null
+++ b/sktime/annotation/ggs.py
@@ -0,0 +1,468 @@
+# -*- coding: utf-8 -*-
+"""
+Greedy Gaussian Segmentation (GGS).
+
+Based on the work from [1]_.
+
+- source code adapted based on: https://github.com/cvxgrp/GGS
+- paper available at: https://stanford.edu/~boyd/papers/pdf/ggs.pdf
+
+References
+----------
+.. [1] Hallac, D., Nystrup, P. & Boyd, S.
+   "Greedy Gaussian segmentation of multivariate time series.",
+    Adv Data Anal Classif 13, 727–751 (2019).
+    https://doi.org/10.1007/s11634-018-0335-0
+"""
+
+import logging
+import math
+from typing import Dict, List, Tuple
+
+import numpy as np
+import numpy.typing as npt
+from attrs import asdict, define, field
+from sklearn.utils.validation import check_random_state
+
+from sktime.base import BaseEstimator
+
+logger = logging.getLogger(__name__)
+
+
+@define
+class GGS:
+    """
+    Greedy Gaussian Segmentation.
+
+    Parameters
+    ----------
+    k_max: int, default=10
+        Maximum number of change points to find
+    lamb: : float, default=1.0
+        Regularization parameter lambda
+    max_shuffles: int, default=250
+        maximum number of shuffles
+    verbose: bool, default=False
+        If ``True`` verbose output is enabled.
+    random_state: int or np.random.RandomState, default=None
+        Either random seed or an instance of ``np.random.RandomState``
+
+    Notes
+    -----
+    Based on the work from [1]_.
+
+    - source code adapted based on: https://github.com/cvxgrp/GGS
+    - paper available at: https://stanford.edu/~boyd/papers/pdf/ggs.pdf
+
+    References
+    ----------
+    .. [1] Hallac, D., Nystrup, P. & Boyd, S.,
+       "Greedy Gaussian segmentation of multivariate time series.",
+       Adv Data Anal Classif 13, 727–751 (2019).
+       https://doi.org/10.1007/s11634-018-0335-0
+    """
+
+    k_max: int = 10
+    lamb: float = 1.0
+    max_shuffles: int = 250
+    verbose: bool = False
+    random_state: int = None
+
+    change_points_: npt.ArrayLike = field(init=False, default=[])
+    _intermediate_change_points: List[List[int]] = field(init=False, default=[])
+    _intermediate_ll: List[float] = field(init=False, default=[])
+
+    def initialize_intermediates(self):
+        """Fit."""
+        self._intermediate_change_points = []
+        self._intermediate_ll = []
+
+    def log_likelihood(self, cov: float, nrows: int, ncols: int) -> float:
+        """
+        Compute the GGS log likelihood.
+
+        Parameters
+        ----------
+        cov: float
+            covariance
+        nrows: int
+            number of observations
+        ncols: int
+            number of columns
+
+        Returns
+        -------
+        log_likelihood
+        """
+        (_, logdet) = np.linalg.slogdet(
+            cov + float(self.lamb) * np.identity(ncols) / nrows
+        )
+
+        return nrows * logdet - float(self.lamb) * np.trace(
+            np.linalg.inv(cov + float(self.lamb) * np.identity(ncols) / nrows)
+        )
+
+    def cumulative_log_likelihood(
+        self, data: npt.ArrayLike, change_points: List[int]
+    ) -> float:
+        """
+        Calculate cumulative GGS log likelihood for all segments.
+
+        Args
+        ----
+        data: array_like
+            time series data
+        change_points: list of ints
+            list of indexes of change points
+
+        Returns
+        -------
+        log_likelihood: cumulative log likelihood
+        """
+        log_likelihood = 0
+        for start, stop in zip(change_points[:-1], change_points[1:]):
+            segment = data[start:stop, :]
+            nrows, ncols = segment.shape
+            cov = np.cov(segment.T, bias=True)
+            log_likelihood -= self.log_likelihood(cov, nrows, ncols)
+        return log_likelihood
+
+    def add_new_change_point(self, data: npt.ArrayLike) -> Tuple[int, float]:
+        """
+        Add change point.
+
+        Parameters
+        ----------
+        data: array_like
+            time series data
+
+        Returns
+        -------
+        index: change point index
+        gll: gained log likelihood
+        """
+        # Initialize parameters
+        m, n = data.shape
+        orig_mean = np.mean(data, axis=0)
+        orig_cov = np.cov(data.T, bias=True)
+        orig_ll = self.log_likelihood(orig_cov, m, n)
+        totSum = m * (orig_cov + np.outer(orig_mean, orig_mean))
+        mu_left = data[0, :] / n
+        mu_right = (m * orig_mean - data[0, :]) / (m - 1)
+        runSum = np.outer(data[0, :], data[0, :])
+        # Loop through all samples
+        # find point where breaking the segment would have the largest LL increase
+        minLL = orig_ll
+        new_index = 0
+        for i in range(2, m - 1):
+            # Update parameters
+            runSum = runSum + np.outer(data[i - 1, :], data[i - 1, :])
+            mu_left = ((i - 1) * mu_left + data[i - 1, :]) / (i)
+            mu_right = ((m - i + 1) * mu_right - data[i - 1, :]) / (m - i)
+            sigLeft = runSum / (i) - np.outer(mu_left, mu_left)
+            sigRight = (totSum - runSum) / (m - i) - np.outer(mu_right, mu_right)
+
+            # Compute Cholesky, LogDet, and Trace
+            Lleft = np.linalg.cholesky(sigLeft + float(self.lamb) * np.identity(n) / i)
+            Lright = np.linalg.cholesky(
+                sigRight + float(self.lamb) * np.identity(n) / (m - i)
+            )
+            llLeft = 2 * sum(map(math.log, np.diag(Lleft)))
+            llRight = 2 * sum(map(math.log, np.diag(Lright)))
+            (trLeft, trRight) = (0, 0)
+            if self.lamb > 0:
+                trLeft = math.pow(np.linalg.norm(np.linalg.inv(Lleft)), 2)
+                trRight = math.pow(np.linalg.norm(np.linalg.inv(Lright)), 2)
+            LL = (
+                i * llLeft
+                - float(self.lamb) * trLeft
+                + (m - i) * llRight
+                - float(self.lamb) * trRight
+            )
+            # Keep track of the best point so far
+            if LL < minLL:
+                minLL = LL
+                new_index = i
+        # Return break, increase in LL
+        return new_index, minLL - orig_ll
+
+    def adjust_change_points(
+        self, data: npt.ArrayLike, change_points: List[int], new_index: List[int]
+    ) -> List[int]:
+        """
+        Adjust change points.
+
+        Parameters
+        ----------
+        data: array_like
+            time series data
+        change_points: list of ints
+            change points indexes
+        new_index: list of ints
+            new change points
+
+        Returns
+        -------
+        change_points: list of ints
+            change point indexes
+        """
+        rng = check_random_state(self.random_state)
+        bp = change_points[:]
+
+        # Just one breakpoint, no need to adjust anything
+        if len(bp) == 3:
+            return bp
+        # Keep track of what change_points have changed,
+        # so that we don't have to adjust ones which we know are constant
+        last_pass = {}
+        this_pass = {b: 0 for b in bp}
+        for i in new_index:
+            this_pass[i] = 1
+        for _ in range(self.max_shuffles):
+            last_pass = dict(this_pass)
+            this_pass = {b: 0 for b in bp}
+            switch_any = False
+            ordering = list(range(1, len(bp) - 1))
+            rng.shuffle(ordering)
+            for i in ordering:
+                # Check if we need to adjust it
+                if (
+                    last_pass[bp[i - 1]] == 1
+                    or last_pass[bp[i + 1]] == 1
+                    or this_pass[bp[i - 1]] == 1
+                    or this_pass[bp[i + 1]] == 1
+                ):
+                    tempData = data[bp[i - 1] : bp[i + 1], :]
+                    ind, val = self.add_new_change_point(tempData)
+                    if bp[i] != ind + bp[i - 1] and val != 0:
+                        last_pass[ind + bp[i - 1]] = last_pass[bp[i]]
+                        del last_pass[bp[i]]
+                        del this_pass[bp[i]]
+                        this_pass[ind + bp[i - 1]] = 1
+                        if self.verbose:
+                            logger.info(
+                                f"Moving {bp[i]} to {ind + bp[i - 1]}"
+                                f"length = {tempData.shape[0]}, {ind}"
+                            )
+                        bp[i] = ind + bp[i - 1]
+                        switch_any = True
+            if not switch_any:
+                return bp
+        return bp
+
+    def initialize_change_points(self, data: npt.ArrayLike) -> List[int]:
+        """Initialize change points."""
+        return [0, data.shape[0] + 1]
+
+    def find_change_points(self, data: npt.ArrayLike) -> List[int]:
+        """
+        Find ``k`` change points on the data at a specific lambda.
+
+        Parameters
+        ----------
+        data: array_like
+            time series data
+
+        Returns
+        -------
+        The K change points, along with all intermediate change points (for k < K)
+        and their corresponding covariance-regularized maximum likelihoods.
+        """
+        change_points = self.initialize_change_points(data)
+        self._intermediate_change_points = [change_points[:]]
+        self._intermediate_ll = [self.cumulative_log_likelihood(data, change_points)]
+
+        # Start GGS Algorithm
+        for _ in range(self.k_max):
+            new_index = -1
+            new_value = +1
+            # For each segment, find change point and increase in LL
+            for start, stop in zip(change_points[:-1], change_points[1:]):
+                segment = data[start:stop, :]
+                ind, val = self.add_new_change_point(segment)
+                if val < new_value:
+                    new_index = ind + start
+                    new_value = val
+
+            # Check if our algorithm is finished
+            if new_value == 0:
+                logger.info("Adding change points!")
+                return change_points
+
+            # Add new change point
+            change_points.append(new_index)
+            change_points.sort()
+            if self.verbose:
+                logger.info(f"Change point occurs at: {new_index}, LL: {new_value}")
+
+            # Adjust current locations of the change points
+            change_points = self.adjust_change_points(data, change_points, [new_index])[
+                :
+            ]
+
+            # Calculate likelihood
+            ll = self.cumulative_log_likelihood(data, change_points)
+            self._intermediate_change_points.append(change_points[:])
+            self._intermediate_ll.append(ll)
+
+        return change_points
+
+    def predict(self, X: npt.ArrayLike):
+        """Predict."""
+        self.change_points_ = self.find_change_points(X)
+
+        labels = np.zeros(X.shape[0], dtype=np.int32)
+        for i, (start, stop) in enumerate(
+            zip(self.change_points_[:-1], self.change_points_[1:])
+        ):
+            labels[start:stop] = i
+        return labels
+
+
+class GreedyGaussianSegmentation(BaseEstimator):
+    """Greedy Gaussian Segmentation Estimator.
+
+    Parameters
+    ----------
+    k_max: int, default=10
+        Maximum number of change points to find
+    lamb: : float, default=1.0
+        Regularization parameter lambda
+    max_shuffles: int, default=250
+        maximum number of shuffles
+    verbose: bool, default=False
+        If ``True`` verbose output is enabled.
+    random_state: int or np.random.RandomState, default=None
+        Either random seed or an instance of ``np.random.RandomState``
+
+    Notes
+    -----
+    Based on the work from [1]_.
+
+    - source code adapted based on: https://github.com/cvxgrp/GGS
+    - paper available at: https://stanford.edu/~boyd/papers/pdf/ggs.pdf
+
+    References
+    ----------
+    .. [1] Hallac, D., Nystrup, P. & Boyd, S.,
+       "Greedy Gaussian segmentation of multivariate time series.",
+       Adv Data Anal Classif 13, 727–751 (2019).
+       https://doi.org/10.1007/s11634-018-0335-0
+    """
+
+    def __init__(
+        self,
+        k_max: int = 10,
+        lamb: float = 1.0,
+        max_shuffles: int = 250,
+        verbose: bool = False,
+        random_state: int = None,
+    ):
+        # this is ugly and necessary only because of dum `test_constructor`
+        self.k_max = k_max
+        self.lamb = lamb
+        self.max_shuffles = max_shuffles
+        self.verbose = verbose
+        self.random_state = random_state
+
+        self._adaptee_class = GGS
+        self._adaptee = self._adaptee_class(
+            k_max=k_max,
+            lamb=lamb,
+            max_shuffles=max_shuffles,
+            verbose=verbose,
+            random_state=random_state,
+        )
+
+    def fit(self, X: npt.ArrayLike, y: npt.ArrayLike = None):
+        """Fit method for compatibility with sklearn-type estimator interface.
+
+        It sets the internal state of the estimator and returns the initialized
+        instance.
+
+        Parameters
+        ----------
+        X: array_like
+            2D `array_like` representing time series with sequence index along
+            the first dimension and value series as columns.
+        y: array_like
+            Placeholder for compatibility with sklearn-api, not used, default=None.
+        """
+        self._adaptee.initialize_intermediates()
+        return self
+
+    def predict(self, X: npt.ArrayLike, y: npt.ArrayLike = None) -> npt.ArrayLike:
+        """Perform segmentation.
+
+        Parameters
+        ----------
+        X: array_like
+            2D `array_like` representing time series with sequence index along
+            the first dimension and value series as columns.
+        y: array_like
+            Placeholder for compatibility with sklearn-api, not used, default=None.
+
+        Returns
+        -------
+        y_pred : array_like
+            1D array with predicted segmentation of the same size as the first
+            dimension of X. The numerical values represent distinct segments
+            labels for each of the data points.
+        """
+        return self._adaptee.predict(X)
+
+    def fit_predict(self, X: npt.ArrayLike, y: npt.ArrayLike = None) -> npt.ArrayLike:
+        """Perform fit and predict.
+
+        Parameters
+        ----------
+        X: array_like
+            2D `array_like` representing time series with sequence index along
+            the first dimension and value series as columns.
+        y: array_like
+            Placeholder for compatibility with sklearn-api, not used, default=None.
+
+        Returns
+        -------
+        y_pred : array_like
+            1D array with predicted segmentation of the same size as the first
+            dimension of X. The numerical values represent distinct segments
+            labels for each of the data points.
+        """
+        return self.fit(X, y).predict(X, y)
+
+    def get_params(self, deep: bool = True) -> Dict:
+        """Return initialization parameters.
+
+        Parameters
+        ----------
+        deep: bool
+            Dummy argument for compatibility with sklearn-api, not used.
+
+        Returns
+        -------
+        params: dict
+            Dictionary with the estimator's initialization parameters, with
+            keys being argument names and values being argument values.
+        """
+        return asdict(self._adaptee, filter=lambda attr, value: attr.init is True)
+
+    def set_params(self, **parameters):
+        """Set the parameters of this object.
+
+        Parameters
+        ----------
+        parameters : dict
+            Initialization parameters for th estimator.
+
+        Returns
+        -------
+        self : reference to self (after parameters have been set)
+        """
+        for key, value in parameters.items():
+            setattr(self._adaptee, key, value)
+        return self
+
+    def __repr__(self) -> str:
+        """Return a string representation of the estimator."""
+        return self._adaptee.__repr__()
diff --git a/sktime/classification/deep_learning/cnn.py b/sktime/classification/deep_learning/cnn.py
index 7fb3678a0..5c1aa0a6a 100644
--- a/sktime/classification/deep_learning/cnn.py
+++ b/sktime/classification/deep_learning/cnn.py
@@ -62,8 +62,8 @@ class CNNClassifier(BaseDeepClassifier):
     >>> from sktime.datasets import load_unit_test
     >>> X_train, y_train = load_unit_test(split="train", return_X_y=True)
     >>> X_test, y_test = load_unit_test(split="test", return_X_y=True)
-    >>> cnn = CNNClassifier()
-    >>> cnn.fit(X_train, y_train)
+    >>> cnn = CNNClassifier()  # doctest: +SKIP
+    >>> cnn.fit(X_train, y_train)  # doctest: +SKIP
     CNNClassifier(...)
     """
 
diff --git a/sktime/classification/dictionary_based/_muse.py b/sktime/classification/dictionary_based/_muse.py
index 9a721016a..e47e66d38 100644
--- a/sktime/classification/dictionary_based/_muse.py
+++ b/sktime/classification/dictionary_based/_muse.py
@@ -12,6 +12,7 @@ import math
 import warnings
 
 import numpy as np
+from joblib import Parallel, delayed
 from numba import njit
 from sklearn.feature_extraction import DictVectorizer
 from sklearn.feature_selection import chi2
@@ -106,7 +107,7 @@ class MUSE(BaseClassifier):
     _tags = {
         "capability:multivariate": True,
         "capability:multithreading": True,
-        "X_inner_mtype": "nested_univ",  # MUSE requires nested datafrane
+        "X_inner_mtype": "nested_univ",  # MUSE requires nested dataframe
         "classifier_type": "dictionary",
     }
 
@@ -181,7 +182,6 @@ class MUSE(BaseClassifier):
 
         self.n_dims = len(self.col_names)
         self.highest_dim_bit = (math.ceil(math.log2(self.n_dims))) + 1
-        self.highest_bits = np.zeros(self.n_dims)
 
         if self.n_dims == 1:
             warnings.warn(
@@ -189,17 +189,19 @@ class MUSE(BaseClassifier):
                 + " multivariate series. It is recommended WEASEL is used instead."
             )
 
-        self.SFA_transformers = [[] for _ in range(self.n_dims)]
+        def _parallel_fit(ind, column):
+            all_words = [
+                [] for x in range(X.shape[0])
+            ]  # no dict needed, array is enough
+            relevant_features_count = 0
 
-        # the words of all dimensions and all time series
-        all_words = [dict() for _ in range(X.shape[0])]
-
-        # On each dimension, perform SFA
-        for ind, column in enumerate(self.col_names):
+            # On each dimension, perform SFA
             X_dim = X[[column]]
             X_dim = from_nested_to_3d_numpy(X_dim)
             series_length = X_dim.shape[-1]  # TODO compute minimum over all ts ?
 
+            SFA_transformers = []
+
             # increment window size in steps of 'win_inc'
             win_inc = self._compute_window_inc(series_length)
 
@@ -213,13 +215,13 @@ class MUSE(BaseClassifier):
                     f"the constructor, but the classifier may not work at "
                     f"all with very short series"
                 )
-            self.window_sizes.append(
+
+            window_sizes = np.array(
                 list(range(self.min_window, self.max_window, win_inc))
             )
+            highest_bits = math.ceil(math.log2(self.max_window)) + 1
 
-            self.highest_bits[ind] = math.ceil(math.log2(self.max_window)) + 1
-
-            for window_size in self.window_sizes[ind]:
+            for window_size in window_sizes:
                 transformer = SFA(
                     word_length=rng.choice(self.word_lengths),
                     alphabet_size=self.alphabet_size,
@@ -236,12 +238,12 @@ class MUSE(BaseClassifier):
 
                 sfa_words = transformer.fit_transform(X_dim, y)
 
-                self.SFA_transformers[ind].append(transformer)
+                SFA_transformers.append(transformer)
                 bag = sfa_words[0]
+                apply_chi_squared = self.p_threshold < 1
 
                 # chi-squared test to keep only relevant features
                 relevant_features = {}
-                apply_chi_squared = self.p_threshold < 1
                 if apply_chi_squared:
                     vectorizer = DictVectorizer(sparse=True, dtype=np.int32, sort=False)
                     bag_vec = vectorizer.fit_transform(bag)
@@ -251,11 +253,12 @@ class MUSE(BaseClassifier):
                     relevant_features = set(
                         np.array(vectorizer.feature_names_)[relevant_features_idx]
                     )
+                    relevant_features_count += len(relevant_features_idx)
 
                 # merging bag-of-patterns of different window_sizes
                 # to single bag-of-patterns with prefix indicating
                 # the used window-length
-                highest = np.int32(self.highest_bits[ind])
+                highest = np.int32(highest_bits)
                 for j in range(len(bag)):
                     for (key, value) in bag[j].items():
                         # chi-squared test
@@ -265,7 +268,43 @@ class MUSE(BaseClassifier):
                             word = MUSE._shift_left(
                                 key, highest, ind, self.highest_dim_bit, window_size
                             )
-                            all_words[j][word] = value
+                            all_words[j].append((word, value))
+
+            return (
+                all_words,
+                relevant_features_count,
+                SFA_transformers,
+                window_sizes,
+                highest_bits,
+            )
+
+        parallel_res = Parallel(n_jobs=self._threads_to_use)(
+            delayed(_parallel_fit)(ind, column)
+            for (ind, column) in enumerate(self.col_names.values)
+        )
+
+        relevant_features_count = 0
+        all_words = [dict() for x in range(len(X))]
+
+        self.SFA_transformers = []
+        self.window_sizes = []
+        self.highest_bits = []
+
+        for (
+            sfa_words,
+            rel_features_count,
+            transformers,
+            window_sizes,
+            highest_bits,
+        ) in parallel_res:
+            relevant_features_count += rel_features_count
+
+            self.window_sizes.append(window_sizes)
+            self.SFA_transformers.append(transformers)
+            self.highest_bits.append(highest_bits)
+
+            for idx, bag in enumerate(sfa_words):
+                all_words[idx].update(bag)
 
         self.clf = make_pipeline(
             DictVectorizer(sparse=True, sort=False),
@@ -325,15 +364,14 @@ class MUSE(BaseClassifier):
         if self.use_first_order_differences:
             X = self._add_first_order_differences(X)
 
-        bag_all_words = [dict() for _ in range(len(X))]
+        def _parallel_transform_words(X, ind, column):
+            bag_all_words = [[] for _ in range(len(X))]
 
-        # On each dimension, perform SFA
-        for ind, column in enumerate(self.col_names):
+            # On each dimension, perform SFA
             X_dim = X[[column]]
             X_dim = from_nested_to_3d_numpy(X_dim)
 
             for i, window_size in enumerate(self.window_sizes[ind]):
-
                 # SFA transform
                 sfa_words = self.SFA_transformers[ind][i].transform(X_dim)
                 bag = sfa_words[0]
@@ -349,9 +387,18 @@ class MUSE(BaseClassifier):
                         word = MUSE._shift_left(
                             key, highest, ind, self.highest_dim_bit, window_size
                         )
-                        bag_all_words[j][word] = value
+                        bag_all_words[j].append((word, value))
+            return bag_all_words
 
-        return bag_all_words
+        parallel_res = Parallel(n_jobs=self._threads_to_use)(
+            delayed(_parallel_transform_words)(X, ind, column)
+            for ind, column in enumerate(self.col_names)
+        )
+        all_words = [dict() for x in range(len(X))]
+        for sfa_words in parallel_res:
+            for idx, bag in enumerate(sfa_words):
+                all_words[idx].update(bag)
+        return all_words
 
     def _add_first_order_differences(self, X):
         X_copy = X.copy()
diff --git a/sktime/classification/dictionary_based/_weasel.py b/sktime/classification/dictionary_based/_weasel.py
index 007065032..bc24324df 100644
--- a/sktime/classification/dictionary_based/_weasel.py
+++ b/sktime/classification/dictionary_based/_weasel.py
@@ -193,7 +193,7 @@ class WEASEL(BaseClassifier):
             window_size,
         ):
             rng = check_random_state(window_size)
-            all_words = [dict() for x in range(len(X))]
+            all_words = [[] for x in range(len(X))]  # no dict needed, array is enough
             relevant_features_count = 0
 
             # for window_size in self.window_sizes:
@@ -209,6 +209,7 @@ class WEASEL(BaseClassifier):
                 remove_repeat_words=False,
                 lower_bounding=False,
                 save_words=False,
+                n_jobs=self._threads_to_use,
             )
 
             sfa_words = transformer.fit_transform(X, y)
@@ -218,6 +219,7 @@ class WEASEL(BaseClassifier):
             apply_chi_squared = self.p_threshold < 1
 
             # chi-squared test to keep only relevant features
+            relevant_features = {}
             if apply_chi_squared:
                 vectorizer = DictVectorizer(sparse=True, dtype=np.int32, sort=False)
                 bag_vec = vectorizer.fit_transform(bag)
@@ -229,21 +231,19 @@ class WEASEL(BaseClassifier):
                 )
                 relevant_features_count += len(relevant_features_idx)
 
-                # merging bag-of-patterns of different window_sizes
-                # to single bag-of-patterns with prefix indicating
-                # the used window-length
-                for j in range(len(bag)):
-                    for (key, value) in bag[j].items():
-                        # chi-squared test
-                        if (not apply_chi_squared) or (key in relevant_features):
-                            # append the prefixes to the words to
-                            # distinguish between window-sizes
-                            word = WEASEL._shift_left(
-                                key, self.highest_bit, window_size
-                            )
-                            all_words[j][word] = value
-
-                return all_words, transformer, relevant_features_count
+            # merging bag-of-patterns of different window_sizes
+            # to single bag-of-patterns with prefix indicating
+            # the used window-length
+            for j in range(len(bag)):
+                for (key, value) in bag[j].items():
+                    # chi-squared test
+                    if (not apply_chi_squared) or (key in relevant_features):
+                        # append the prefixes to the words to
+                        # distinguish between window-sizes
+                        word = WEASEL._shift_left(key, self.highest_bit, window_size)
+                        all_words[j].append((word, value))
+
+            return all_words, transformer, relevant_features_count
 
         parallel_res = Parallel(n_jobs=self._threads_to_use)(
             delayed(_parallel_fit)(window_size) for window_size in self.window_sizes
@@ -253,13 +253,11 @@ class WEASEL(BaseClassifier):
         all_words = [dict() for x in range(len(X))]
 
         for sfa_words, transformer, rel_features_count in parallel_res:
-            transformer.n_jobs = self._threads_to_use
             self.SFA_transformers.append(transformer)
             relevant_features_count += rel_features_count
 
             for idx, bag in enumerate(sfa_words):
-                for word, count in bag.items():
-                    all_words[idx][word] = count
+                all_words[idx].update(bag)
 
         self.clf = make_pipeline(
             DictVectorizer(sparse=True, sort=False),
@@ -313,8 +311,9 @@ class WEASEL(BaseClassifier):
         return self.clf.predict_proba(bag)
 
     def _transform_words(self, X):
-        bag_all_words = [dict() for _ in range(len(X))]
-        for transformer in self.SFA_transformers:
+        def _parallel_transform_words(X, transformer):
+            bag_all_words = [[] for x in range(len(X))]
+
             # SFA transform
             sfa_words = transformer.transform(X)
             bag = sfa_words[0]
@@ -324,14 +323,25 @@ class WEASEL(BaseClassifier):
             # the used window-length
             for j in range(len(bag)):
                 for (key, value) in bag[j].items():
+                    # if value > 1:
                     # append the prefices to the words to distinguish
                     # between window-sizes
                     word = WEASEL._shift_left(
                         key, self.highest_bit, transformer.window_size
                     )
-                    bag_all_words[j][word] = value
+                    bag_all_words[j].append((word, value))
+
+            return bag_all_words
 
-        return bag_all_words
+        parallel_res = Parallel(n_jobs=self._threads_to_use)(
+            delayed(_parallel_transform_words)(X, transformer)
+            for transformer in self.SFA_transformers
+        )
+        all_words = [dict() for x in range(len(X))]
+        for sfa_words in parallel_res:
+            for idx, bag in enumerate(sfa_words):
+                all_words[idx].update(bag)
+        return all_words
 
     def _compute_window_inc(self):
         win_inc = self.window_inc
diff --git a/sktime/datasets/_single_problem_loaders.py b/sktime/datasets/_single_problem_loaders.py
index e1d014cd3..85b3d07d3 100644
--- a/sktime/datasets/_single_problem_loaders.py
+++ b/sktime/datasets/_single_problem_loaders.py
@@ -38,6 +38,7 @@ __all__ = [
 ]
 
 import os
+from urllib.error import HTTPError
 
 import numpy as np
 import pandas as pd
@@ -1053,31 +1054,54 @@ def load_solar(
 
     Examples
     --------
-    >>> from sktime.datasets import load_solar
-    >>> y = load_solar()
+    >>> from sktime.datasets import load_solar  # doctest: +SKIP
+    >>> y = load_solar()  # doctest: +SKIP
     """
-    url = "https://api0.solar.sheffield.ac.uk/pvlive/api/"
-    url = url + api_version + "/gsp/0?"
-    url = url + "start=" + start + "T00:00:00&"
-    url = url + "end=" + end + "T00:00:00&"
-    url = url + "extra_fields=capacity_mwp&"
-    url = url + "data_format=csv"
-
-    df = (
-        pd.read_csv(
-            url, index_col=["gsp_id", "datetime_gmt"], parse_dates=["datetime_gmt"]
+    from sktime.utils.validation._dependencies import _check_soft_dependencies
+
+    _check_soft_dependencies("backoff")
+
+    import backoff
+
+    @backoff.on_exception(wait_gen=backoff.expo, exception=HTTPError, max_tries=5)
+    def _load_solar(
+        start="2021-05-01",
+        end="2021-09-01",
+        normalise=True,
+        return_full_df=False,
+        api_version="v4",
+    ):
+        """Private loader, for decoration with backoff."""
+        url = "https://api0.solar.sheffield.ac.uk/pvlive/api/"
+        url = url + api_version + "/gsp/0?"
+        url = url + "start=" + start + "T00:00:00&"
+        url = url + "end=" + end + "T00:00:00&"
+        url = url + "extra_fields=capacity_mwp&"
+        url = url + "data_format=csv"
+
+        df = (
+            pd.read_csv(
+                url, index_col=["gsp_id", "datetime_gmt"], parse_dates=["datetime_gmt"]
+            )
+            .droplevel(0)
+            .sort_index()
         )
-        .droplevel(0)
-        .sort_index()
-    )
-    df = df.asfreq("30T")
-    df["generation_pu"] = df["generation_mw"] / df["capacity_mwp"]
-
-    if return_full_df:
+        df = df.asfreq("30T")
         df["generation_pu"] = df["generation_mw"] / df["capacity_mwp"]
-        return df
-    else:
-        if normalise:
-            return df["generation_pu"].rename("solar_gen")
+
+        if return_full_df:
+            df["generation_pu"] = df["generation_mw"] / df["capacity_mwp"]
+            return df
         else:
-            return df["generation_mw"].rename("solar_gen")
+            if normalise:
+                return df["generation_pu"].rename("solar_gen")
+            else:
+                return df["generation_mw"].rename("solar_gen")
+
+    return _load_solar(
+        start=start,
+        end=end,
+        normalise=normalise,
+        return_full_df=return_full_df,
+        api_version=api_version,
+    )
diff --git a/sktime/datatypes/_series/_check.py b/sktime/datatypes/_series/_check.py
index 017f697e6..9c98b26cd 100644
--- a/sktime/datatypes/_series/_check.py
+++ b/sktime/datatypes/_series/_check.py
@@ -40,6 +40,7 @@ __all__ = ["check_dict"]
 import numpy as np
 import pandas as pd
 
+from sktime.utils.validation._dependencies import _check_soft_dependencies
 from sktime.utils.validation.series import is_in_valid_index_types
 
 VALID_INDEX_TYPES = (pd.RangeIndex, pd.PeriodIndex, pd.DatetimeIndex)
@@ -235,3 +236,71 @@ def _index_equally_spaced(index):
     all_equal = np.all(diffs == diffs[0])
 
     return all_equal
+
+
+if _check_soft_dependencies("xarray", severity="none"):
+    import xarray as xr
+
+    def check_xrdataarray_series(obj, return_metadata=False, var_name="obj"):
+        metadata = {}
+
+        def ret(valid, msg, metadata, return_metadata):
+            if return_metadata:
+                return valid, msg, metadata
+            return valid
+
+        if not isinstance(obj, xr.DataArray):
+            msg = f"{var_name} must be a xarray.DataArray, found {type(obj)}"
+            return ret(False, msg, None, return_metadata)
+
+        # we now know obj is a xr.DataArray
+        if len(obj.dims) > 2:  # Without multi indexing only two dimensions are possible
+            msg = f"{var_name} must have two or less dimension, found {type(obj.dims)}"
+            return ret(False, msg, None, return_metadata)
+
+        # The first dimension is the index of the time series in sktimelen
+        index = obj.indexes[obj.dims[0]]
+
+        metadata["is_empty"] = len(index) < 1 or len(obj.values) < 1
+        # The second dimension is the set of columns
+        metadata["is_univariate"] = len(obj.dims) == 1 or len(obj[obj.dims[1]]) < 2
+
+        # check that columns are unique
+        msg = f"{var_name} must have " f"unique column indices, but found {obj.dims}"
+        assert len(obj.dims) == len(set(obj.dims)), msg
+
+        # check whether the time index is of valid type
+        if not is_in_valid_index_types(index):
+            msg = (
+                f"{type(index)} is not supported for {var_name}, use "
+                f"one of {VALID_INDEX_TYPES} or integer index instead."
+            )
+            return ret(False, msg, None, return_metadata)
+
+        # check that the dtype is not object
+        if "object" == obj.dtype:
+            msg = f"{var_name} should not have column of 'object' dtype"
+            return ret(False, msg, None, return_metadata)
+
+        # Check time index is ordered in time
+        if not index.is_monotonic:
+            msg = (
+                f"The (time) index of {var_name} must be sorted "
+                f"monotonically increasing, but found: {index}"
+            )
+            return ret(False, msg, None, return_metadata)
+
+        if FREQ_SET_CHECK and isinstance(index, pd.DatetimeIndex):
+            if index.freq is None:
+                msg = f"{var_name} has DatetimeIndex, but no freq attribute set."
+                return ret(False, msg, None, return_metadata)
+
+        # check whether index is equally spaced or if there are any nans
+        #   compute only if needed
+        if return_metadata:
+            metadata["is_equally_spaced"] = _index_equally_spaced(index)
+            metadata["has_nans"] = obj.isnull().values.any()
+
+        return ret(True, None, metadata, return_metadata)
+
+    check_dict[("xr.DataArray", "Series")] = check_xrdataarray_series
diff --git a/sktime/datatypes/_series/_convert.py b/sktime/datatypes/_series/_convert.py
index fea20e113..1ee68eac9 100644
--- a/sktime/datatypes/_series/_convert.py
+++ b/sktime/datatypes/_series/_convert.py
@@ -37,6 +37,8 @@ import pandas as pd
 ##############################################################
 # methods to convert one machine type to another machine type
 ##############################################################
+from sktime.datatypes._registry import MTYPE_LIST_SERIES
+from sktime.utils.validation._dependencies import _check_soft_dependencies
 
 convert_dict = dict()
 
@@ -47,7 +49,7 @@ def convert_identity(obj, store=None):
 
 
 # assign identity function to type conversion to self
-for tp in ["pd.Series", "pd.DataFrame", "np.ndarray"]:
+for tp in ["pd.Series", "pd.DataFrame", "np.ndarray", "xr.DataArray"]:
     convert_dict[(tp, tp, "Series")] = convert_identity
 
 
@@ -172,3 +174,70 @@ def convert_np_to_UvS_as_Series(obj: np.ndarray, store=None) -> pd.Series:
 
 
 convert_dict[("np.ndarray", "pd.Series", "Series")] = convert_np_to_UvS_as_Series
+
+
+# obtain other conversions from/to numpyflat via concatenation to numpy3D
+def _concat(fun1, fun2):
+    def concat_fun(obj, store=None):
+        obj1 = fun1(obj, store=store)
+        obj2 = fun2(obj1, store=store)
+        return obj2
+
+    return concat_fun
+
+
+def _extend_conversions(mtype, anchor_mtype, convert_dict):
+    keys = convert_dict.keys()
+    scitype = list(keys)[0][2]
+
+    for tp in set(MTYPE_LIST_SERIES).difference([mtype, anchor_mtype]):
+        if (anchor_mtype, tp, scitype) in convert_dict.keys():
+            convert_dict[(mtype, tp, scitype)] = _concat(
+                convert_dict[(mtype, anchor_mtype, scitype)],
+                convert_dict[(anchor_mtype, tp, scitype)],
+            )
+        if (tp, anchor_mtype, scitype) in convert_dict.keys():
+            convert_dict[(tp, mtype, scitype)] = _concat(
+                convert_dict[(tp, anchor_mtype, scitype)],
+                convert_dict[(anchor_mtype, mtype, scitype)],
+            )
+
+
+if _check_soft_dependencies("xarray", severity="none"):
+    import xarray as xr
+
+    def convert_xrdataarray_to_Mvs_as_Series(
+        obj: xr.DataArray, store=None
+    ) -> pd.DataFrame:
+        if not isinstance(obj, xr.DataArray):
+            raise TypeError("input must be a xr.DataArray")
+
+        if isinstance(store, dict):
+            store["coords"] = list(obj.coords.keys())
+
+        index = obj.indexes[obj.dims[0]]
+        columns = obj.indexes[obj.dims[1]] if len(obj.dims) == 2 else None
+        return pd.DataFrame(obj.values, index=index, columns=columns)
+
+    convert_dict[
+        ("xr.DataArray", "pd.DataFrame", "Series")
+    ] = convert_xrdataarray_to_Mvs_as_Series
+
+    def convert_Mvs_to_xrdatarray_as_Series(
+        obj: pd.DataFrame, store=None
+    ) -> xr.DataArray:
+        if not isinstance(obj, pd.DataFrame):
+            raise TypeError("input must be a xr.DataArray")
+
+        result = xr.DataArray(obj.values, coords=[obj.index, obj.columns])
+        if isinstance(store, dict) and "coords" in store:
+            result = result.rename(
+                dict(zip(list(result.coords.keys()), store["coords"]))
+            )
+        return result
+
+    convert_dict[
+        ("pd.DataFrame", "xr.DataArray", "Series")
+    ] = convert_Mvs_to_xrdatarray_as_Series
+
+    _extend_conversions("xr.DataArray", "pd.DataFrame", convert_dict)
diff --git a/sktime/datatypes/_series/_examples.py b/sktime/datatypes/_series/_examples.py
index 42296f846..7b635ee59 100644
--- a/sktime/datatypes/_series/_examples.py
+++ b/sktime/datatypes/_series/_examples.py
@@ -32,6 +32,8 @@ overall, conversions from non-lossy representations to any other ones
 import numpy as np
 import pandas as pd
 
+from sktime.utils.validation._dependencies import _check_soft_dependencies
+
 example_dict = dict()
 example_dict_lossy = dict()
 example_dict_metadata = dict()
@@ -54,6 +56,18 @@ arr = np.array([[1], [4], [0.5], [-3]])
 example_dict[("np.ndarray", "Series", 0)] = arr
 example_dict_lossy[("np.ndarray", "Series", 0)] = True
 
+if _check_soft_dependencies("xarray", severity="none"):
+    import xarray as xr
+
+    da = xr.DataArray(
+        [[1], [4], [0.5], [-3]],
+        coords=[[0, 1, 2, 3], ["a"]],
+    )
+
+    example_dict[("xr.DataArray", "Series", 0)] = da
+    example_dict_lossy[("xr.DataArray", "Series", 0)] = False
+
+
 example_dict_metadata[("Series", 0)] = {
     "is_univariate": True,
     "is_equally_spaced": True,
@@ -76,6 +90,16 @@ arr = np.array([[1, 3], [4, 7], [0.5, 2], [-3, -3 / 7]])
 
 example_dict[("np.ndarray", "Series", 1)] = arr
 example_dict_lossy[("np.ndarray", "Series", 1)] = True
+if _check_soft_dependencies("xarray", severity="none"):
+    import xarray as xr
+
+    da = xr.DataArray(
+        [[1, 3], [4, 7], [0.5, 2], [-3, -3 / 7]],
+        coords=[[0, 1, 2, 3], ["a", "b"]],
+    )
+
+    example_dict[("xr.DataArray", "Series", 1)] = da
+    example_dict_lossy[("xr.DataArray", "Series", 1)] = False
 
 example_dict_metadata[("Series", 1)] = {
     "is_univariate": False,
@@ -100,6 +124,18 @@ arr = np.array([[1, 3], [4, 7], [0.5, 2], [3, 3 / 7]])
 example_dict[("np.ndarray", "Series", 2)] = arr
 example_dict_lossy[("np.ndarray", "Series", 2)] = True
 
+if _check_soft_dependencies("xarray", severity="none"):
+    import xarray as xr
+
+    da = xr.DataArray(
+        [[1, 3], [4, 7], [0.5, 2], [3, 3 / 7]],
+        coords=[[0, 1, 2, 3], ["a", "b"]],
+    )
+
+    example_dict[("xr.DataArray", "Series", 2)] = da
+    example_dict_lossy[("xr.DataArray", "Series", 2)] = False
+
+
 example_dict_metadata[("Series", 2)] = {
     "is_univariate": False,
     "is_equally_spaced": True,
diff --git a/sktime/datatypes/_series/_registry.py b/sktime/datatypes/_series/_registry.py
index 73f76a489..a04007d19 100644
--- a/sktime/datatypes/_series/_registry.py
+++ b/sktime/datatypes/_series/_registry.py
@@ -20,6 +20,11 @@ MTYPE_REGISTER_SERIES = [
         "Series",
         "2D numpy.ndarray with rows=samples, cols=variables, index=integers",
     ),
+    (
+        "xr.DataArray",
+        "Series",
+        "xr.DataArray representation of a uni- or multivariate series",
+    ),
 ]
 
 MTYPE_LIST_SERIES = pd.DataFrame(MTYPE_REGISTER_SERIES)[0].values
diff --git a/sktime/dists_kernels/scipy_dist.py b/sktime/dists_kernels/scipy_dist.py
index 5ff618db5..a49b15bf9 100644
--- a/sktime/dists_kernels/scipy_dist.py
+++ b/sktime/dists_kernels/scipy_dist.py
@@ -7,8 +7,8 @@ Interface module to scipy.spatial's pairwise distance function cdist
 
 __author__ = ["fkiraly"]
 
+import numpy as np
 import pandas as pd
-
 from scipy.spatial.distance import cdist
 
 from sktime.dists_kernels._base import BasePairwiseTransformer
@@ -23,15 +23,15 @@ class ScipyDist(BasePairwiseTransformer):
 
     Parameters
     ----------
-    metric: string or function, as in cdist; default = 'euclidean'
+    metric : string or function, as in cdist; default = 'euclidean'
         if string, one of: 'braycurtis', 'canberra', 'chebyshev', 'cityblock',
             'correlation', 'cosine', 'dice', 'euclidean', 'hamming', 'jaccard',
             'jensenshannon', 'kulsinski', 'mahalanobis', 'matching', 'minkowski',
             'rogerstanimoto', 'russellrao', 'seuclidean', 'sokalmichener',
             'sokalsneath', 'sqeuclidean', 'yule'
         if function, should have signature 1D-np.array x 1D-np.array -> float
-    p: if metric='minkowski', the "p" in "p-norm", otherwise irrelevant
-    colalign: string, one of 'intersect' (default), 'force-align', 'none'
+    p:  if metric='minkowski', the "p" in "p-norm", otherwise irrelevant
+    colalign : string, one of 'intersect' (default), 'force-align', 'none'
         controls column alignment if X, X2 passed in fit are pd.DataFrame
         columns between X and X2 are aligned via column names
         if 'intersect', distance is computed on columns occurring both in X and X2,
@@ -40,17 +40,34 @@ class ScipyDist(BasePairwiseTransformer):
             column ordering in X2 is copied from X
         if 'none', X and X2 are passed through unmodified (no columns are aligned)
             note: this will potentially align "non-matching" columns
+    var_weights : 1D np.array of float or None, default=None
+        weight/scaling vector applied to variables in X/X2
+        before being passed to cdist, i-th col of X/X2 is multiplied by var_weights[i]
+        if None, equivalent to all-ones vector
+    metric_kwargs : dict, optional, default=None
+        any kwargs passed to the metric in addition, i.e., to the function cdist
+        common kwargs: "w" : array-like, same length as X.columns, weights for metric
+        refer to scipy.spatial.distance.dist for a documentation of other extra kwargs
     """
 
     _tags = {
         "symmetric": True,  # all the distances are symmetric
     }
 
-    def __init__(self, metric="euclidean", p=2, colalign="intersect"):
+    def __init__(
+        self,
+        metric="euclidean",
+        p=2,
+        colalign="intersect",
+        var_weights=None,
+        metric_kwargs=None,
+    ):
 
         self.metric = metric
         self.p = p
         self.colalign = colalign
+        self.var_weights = var_weights
+        self.metric_kwargs = metric_kwargs
 
         super(ScipyDist, self).__init__()
 
@@ -78,6 +95,10 @@ class ScipyDist(BasePairwiseTransformer):
         """
         p = self.p
         metric = self.metric
+        var_weights = self.var_weights
+        metric_kwargs = self.metric_kwargs
+        if metric_kwargs is None:
+            metric_kwargs = {}
 
         if isinstance(X, pd.DataFrame):
             X = X.select_dtypes("number").to_numpy(dtype="float")
@@ -85,9 +106,45 @@ class ScipyDist(BasePairwiseTransformer):
         if isinstance(X2, pd.DataFrame):
             X2 = X2.select_dtypes("number").to_numpy(dtype="float")
 
-        if metric == "minkowski":
-            distmat = cdist(XA=X, XB=X2, metric=metric, p=p)
+        if np.ndim(var_weights) == 1:
+            if len(var_weights) == len(X.columns) == len(X2.columns):
+                X = var_weights * X
+                X2 = var_weights * X2
+            else:
+                raise ValueError(
+                    "weights vector length must be equal to X and X2 number of columns"
+                )
+
+        if metric == "minkowski" and "p" not in metric_kwargs.keys():
+            distmat = cdist(XA=X, XB=X2, metric=metric, p=p, **metric_kwargs)
         else:
-            distmat = cdist(XA=X, XB=X2, metric=metric)
+            distmat = cdist(XA=X, XB=X2, metric=metric, **metric_kwargs)
 
         return distmat
+
+    @classmethod
+    def get_test_params(cls, parameter_set="default"):
+        """Return testing parameter settings for the estimator.
+
+        Parameters
+        ----------
+        parameter_set : str, default="default"
+            Name of the set of test parameters to return, for use in tests. If no
+            special parameters are defined for a value, will return `"default"` set.
+            There are currently no reserved values for distance/kernel transformers.
+
+        Returns
+        -------
+        params : dict or list of dict, default = {}
+            Parameters to create testing instances of the class
+            Each dict are parameters to construct an "interesting" test instance, i.e.,
+            `MyClass(**params)` or `MyClass(**params[i])` creates a valid test instance.
+            `create_test_instance` uses the first (or only) dictionary in `params`
+        """
+        # default settings
+        params1 = {}
+
+        # using kwargs
+        params2 = {"metric": "minkowski", "p": 3}
+
+        return [params1, params2]
diff --git a/sktime/forecasting/model_selection/_split.py b/sktime/forecasting/model_selection/_split.py
index 0b57a79aa..c796443ee 100644
--- a/sktime/forecasting/model_selection/_split.py
+++ b/sktime/forecasting/model_selection/_split.py
@@ -1106,7 +1106,7 @@ class SlidingWindowSplitter(BaseWindowSplitter):
         Step length between windows
     initial_window : int or timedelta or pd.DateOffset, optional (default=None)
         Window length of first window
-    start_with_window : bool, optional (default=False)
+    start_with_window : bool, optional (default=True)
         - If True, starts with full window.
         - If False, starts with empty window.
 
diff --git a/sktime/transformations/base.py b/sktime/transformations/base.py
index 7f3ae85af..d003db519 100644
--- a/sktime/transformations/base.py
+++ b/sktime/transformations/base.py
@@ -145,9 +145,10 @@ class BaseTransformer(BaseEstimator):
         "pd_multiindex_hier",
     ]
 
-    def __init__(self):
+    def __init__(self, _output_convert="auto"):
 
         self._converter_store_X = dict()  # storage dictionary for in/output conversion
+        self._output_convert = _output_convert
 
         super(BaseTransformer, self).__init__()
         _check_estimator_deps(self)
@@ -469,7 +470,10 @@ class BaseTransformer(BaseEstimator):
             Xt = self._vectorize("transform", X=X_inner, y=y_inner)
 
         # convert to output mtype
-        X_out = self._convert_output(Xt, metadata=metadata)
+        if not hasattr(self, "_output_convert") or self._output_convert == "auto":
+            X_out = self._convert_output(Xt, metadata=metadata)
+        else:
+            X_out = Xt
 
         return X_out
 
@@ -582,7 +586,10 @@ class BaseTransformer(BaseEstimator):
             Xt = self._vectorize("inverse_transform", X=X_inner, y=y_inner)
 
         # convert to output mtype
-        X_out = self._convert_output(Xt, metadata=metadata, inverse=True)
+        if self._output_convert == "auto":
+            X_out = self._convert_output(Xt, metadata=metadata, inverse=True)
+        else:
+            X_out = Xt
 
         return X_out
 
diff --git a/sktime/transformations/compose.py b/sktime/transformations/compose.py
index ff561bb67..cca2d9ea5 100644
--- a/sktime/transformations/compose.py
+++ b/sktime/transformations/compose.py
@@ -1146,7 +1146,7 @@ class InvertTransform(_DelegatedTransformer):
                 "If the tag was correctly set, this transformer will likely crash"
             )
         inner_output = transformer.get_tag("scitype:transform-output")
-        if not transformer.get_tag("scitype:transform-output") != "Series":
+        if transformer.get_tag("scitype:transform-output") != "Series":
             warn(
                 f"transformer output is not Series but {inner_output}, "
                 "according to scitype:transform-output tag. "
diff --git a/sktime/transformations/panel/dictionary_based/_paa.py b/sktime/transformations/panel/dictionary_based/_paa.py
index 4a0718d5f..9c1ef2d98 100644
--- a/sktime/transformations/panel/dictionary_based/_paa.py
+++ b/sktime/transformations/panel/dictionary_based/_paa.py
@@ -1,14 +1,16 @@
 # -*- coding: utf-8 -*-
+"""Piecewise Aggregate Approximation Transformer (PAA)."""
 import pandas as pd
-from sktime.transformations.base import _PanelToPanelTransformer
+
 from sktime.datatypes._panel._convert import from_nested_to_2d_array
-from sktime.utils.validation.panel import check_X
+from sktime.transformations.base import BaseTransformer
 
-__author__ = "Matthew Middlehurst"
+__author__ = "MatthewMiddlehurst"
 
 
-class PAA(_PanelToPanelTransformer):
-    """
+class PAA(BaseTransformer):
+    """Piecewise Aggregate Approximation Transformer (PAA).
+
     (PAA) Piecewise Aggregate Approximation Transformer, as described in
     Eamonn Keogh, Kaushik Chakrabarti, Michael Pazzani, and Sharad Mehrotra.
     Dimensionality reduction for fast similarity search in large time series
@@ -25,18 +27,30 @@ class PAA(_PanelToPanelTransformer):
     Parameters
     ----------
     num_intervals   : int, dimension of the transformed data (default 8)
-
     """
 
+    _tags = {
+        "scitype:transform-input": "Series",
+        # what is the scitype of X: Series, or Panel
+        "scitype:transform-output": "Series",
+        # what scitype is returned: Primitives, Series, Panel
+        "scitype:instancewise": True,  # is this an instance-wise transform?
+        "X_inner_mtype": "nested_univ",  # which mtypes do _fit/_predict support for X?
+        "y_inner_mtype": "None",  # which mtypes do _fit/_predict support for X?
+    }
+
     def __init__(self, num_intervals=8):
         self.num_intervals = num_intervals
         super(PAA, self).__init__()
 
     def set_num_intervals(self, n):
+        """Set self.num_intervals to n."""
         self.num_intervals = n
 
-    def transform(self, X, y=None):
-        """
+    # todo: looks like this just loops over series instances
+    # so should be refactored to work on Series directly
+    def _transform(self, X, y=None):
+        """Transform data.
 
         Parameters
         ----------
@@ -48,10 +62,6 @@ class PAA(_PanelToPanelTransformer):
         dims: Pandas data frame with first dimension in column zero,
               second in column one etc.
         """
-        # Check the data
-        self.check_is_fitted()
-        X = check_X(X, enforce_univariate=False, coerce_to_pandas=True)
-
         # Get information about the dataframe
         num_atts = len(X.iloc[0, 0])
         col_names = X.columns
@@ -115,7 +125,8 @@ class PAA(_PanelToPanelTransformer):
         return dims
 
     def _check_parameters(self, num_atts):
-        """
+        """Check parameters of PAA.
+
         Function for checking the values of parameters inserted into PAA.
         For example, the number of subsequences cannot be larger than the
         time series length.
diff --git a/sktime/transformations/panel/dictionary_based/_sax.py b/sktime/transformations/panel/dictionary_based/_sax.py
index 39c2901e1..6782c7fc4 100644
--- a/sktime/transformations/panel/dictionary_based/_sax.py
+++ b/sktime/transformations/panel/dictionary_based/_sax.py
@@ -1,25 +1,26 @@
 # -*- coding: utf-8 -*-
+"""Symbolic Aggregate approXimation (SAX) transformer."""
+
 import sys
 
 import numpy as np
 import pandas as pd
 import scipy.stats
 
-from sktime.transformations.base import _PanelToPanelTransformer
+from sktime.transformations.base import BaseTransformer
 from sktime.transformations.panel.dictionary_based import PAA
 
 #    TO DO: verify this returned pandas is consistent with sktime
 #    definition. Timestamps?
-from sktime.utils.validation.panel import check_X
 
 # from numba import types
 # from numba.experimental import jitclass
 
-__author__ = "Matthew Middlehurst"
+__author__ = "MatthewMiddlehurst"
 
 
-class SAX(_PanelToPanelTransformer):
-    """SAX (Symbolic Aggregate approXimation) transformer.
+class SAX(BaseTransformer):
+    """Symbolic Aggregate approXimation (SAX) transformer.
 
     as described in
     Jessica Lin, Eamonn Keogh, Li Wei and Stefano Lonardi,
@@ -38,29 +39,38 @@ class SAX(_PanelToPanelTransformer):
 
     Parameters
     ----------
-        word_length:         int, length of word to shorten window to (using
-        PAA) (default 8)
-        alphabet_size:       int, number of values to discretise each value
-        to (default to 4)
-        window_size:         int, size of window for sliding. Input series
-        length for whole series transform (default to 12)
-        remove_repeat_words: boolean, whether to use numerosity reduction (
-        default False)
-        save_words:          boolean, whether to use numerosity reduction (
-        default False)
-
-        return_pandas_data_series:          boolean, default = True
-            set to true to return Pandas Series as a result of transform.
-            setting to true reduces speed significantly but is required for
-            automatic test.
+    word_length:         int, length of word to shorten window to (using
+    PAA) (default 8)
+    alphabet_size:       int, number of values to discretise each value
+    to (default to 4)
+    window_size:         int, size of window for sliding. Input series
+    length for whole series transform (default to 12)
+    remove_repeat_words: boolean, whether to use numerosity reduction (
+    default False)
+    save_words:          boolean, whether to use numerosity reduction (
+    default False)
+
+    return_pandas_data_series:          boolean, default = True
+        set to true to return Pandas Series as a result of transform.
+        setting to true reduces speed significantly but is required for
+        automatic test.
 
     Attributes
     ----------
-        words:      history = []
-
+    words:      history = []
     """
 
-    _tags = {"univariate-only": True, "fit_is_empty": True}
+    _tags = {
+        "univariate-only": True,
+        "fit_is_empty": True,
+        "scitype:transform-input": "Series",
+        # what is the scitype of X: Series, or Panel
+        "scitype:transform-output": "Series",
+        # what scitype is returned: Primitives, Series, Panel
+        "scitype:instancewise": True,  # is this an instance-wise transform?
+        "X_inner_mtype": "numpy3D",  # which mtypes do _fit/_predict support for X?
+        "y_inner_mtype": "None",  # which mtypes do _fit/_predict require for y?
+    }
 
     def __init__(
         self,
@@ -79,9 +89,11 @@ class SAX(_PanelToPanelTransformer):
         self.return_pandas_data_series = return_pandas_data_series
         self.words = []
 
-        super(SAX, self).__init__()
+        super(SAX, self).__init__(_output_convert="off")
 
-    def transform(self, X, y=None):
+    # todo: looks like this just loops over series instances
+    # so should be refactored to work on Series directly
+    def _transform(self, X, y=None):
         """Transform data.
 
         Parameters
@@ -93,8 +105,6 @@ class SAX(_PanelToPanelTransformer):
         -------
         dims: Pandas data frame with first dimension in column zero
         """
-        self.check_is_fitted()
-        X = check_X(X, enforce_univariate=True, coerce_to_numpy=True)
         X = X.squeeze(1)
 
         if self.alphabet_size < 2 or self.alphabet_size > 4:
diff --git a/sktime/transformations/panel/dictionary_based/_sfa.py b/sktime/transformations/panel/dictionary_based/_sfa.py
index a045157e9..e3932f21b 100644
--- a/sktime/transformations/panel/dictionary_based/_sfa.py
+++ b/sktime/transformations/panel/dictionary_based/_sfa.py
@@ -4,7 +4,7 @@
 Configurable SFA transform for discretising time series into words.
 """
 
-__author__ = ["Matthew Middlehurst", "Patrick Schäfer"]
+__author__ = ["MatthewMiddlehurst", "patrickzib"]
 __all__ = ["SFA"]
 
 import math
@@ -20,7 +20,7 @@ from sklearn.feature_selection import f_classif
 from sklearn.preprocessing import KBinsDiscretizer
 from sklearn.tree import DecisionTreeClassifier
 
-from sktime.transformations.base import _PanelToPanelTransformer
+from sktime.transformations.base import BaseTransformer
 from sktime.utils.validation.panel import check_X
 
 # The binning methods to use: equi-depth, equi-width, information gain or kmeans
@@ -29,7 +29,7 @@ binning_methods = {"equi-depth", "equi-width", "information-gain", "kmeans"}
 # TODO remove imag-part from dc-component component
 
 
-class SFA(_PanelToPanelTransformer):
+class SFA(BaseTransformer):
     """Symbolic Fourier Approximation (SFA) Transformer.
 
     Overview: for each series:
@@ -43,51 +43,51 @@ class SFA(_PanelToPanelTransformer):
 
     Parameters
     ----------
-        word_length:         int, default = 8
-            length of word to shorten window to (using PAA)
+    word_length:         int, default = 8
+        length of word to shorten window to (using PAA)
 
-        alphabet_size:       int, default = 4
-            number of values to discretise each value to
+    alphabet_size:       int, default = 4
+        number of values to discretise each value to
 
-        window_size:         int, default = 12
-            size of window for sliding. Input series
-            length for whole series transform
+    window_size:         int, default = 12
+        size of window for sliding. Input series
+        length for whole series transform
 
-        norm:                boolean, default = False
-            mean normalise words by dropping first fourier coefficient
+    norm:                boolean, default = False
+        mean normalise words by dropping first fourier coefficient
 
-        binning_method:      {"equi-depth", "equi-width", "information-gain", "kmeans"},
-                             default="equi-depth"
-            the binning method used to derive the breakpoints.
+    binning_method:      {"equi-depth", "equi-width", "information-gain", "kmeans"},
+                            default="equi-depth"
+        the binning method used to derive the breakpoints.
 
-        anova:               boolean, default = False
-            If True, the Fourier coefficient selection is done via a one-way
-            ANOVA test. If False, the first Fourier coefficients are selected.
-            Only applicable if labels are given
+    anova:               boolean, default = False
+        If True, the Fourier coefficient selection is done via a one-way
+        ANOVA test. If False, the first Fourier coefficients are selected.
+        Only applicable if labels are given
 
-        bigrams:             boolean, default = False
-            whether to create bigrams of SFA words
+    bigrams:             boolean, default = False
+        whether to create bigrams of SFA words
 
-        skip_grams:          boolean, default = False
-            whether to create skip-grams of SFA words
+    skip_grams:          boolean, default = False
+        whether to create skip-grams of SFA words
 
-        remove_repeat_words: boolean, default = False
-            whether to use numerosity reduction (default False)
+    remove_repeat_words: boolean, default = False
+        whether to use numerosity reduction (default False)
 
-        levels:              int, default = 1
-            Number of spatial pyramid levels
+    levels:              int, default = 1
+        Number of spatial pyramid levels
 
-        save_words:          boolean, default = False
-            whether to save the words generated for each series (default False)
+    save_words:          boolean, default = False
+        whether to save the words generated for each series (default False)
 
-        return_pandas_data_series:          boolean, default = False
-            set to true to return Pandas Series as a result of transform.
-            setting to true reduces speed significantly but is required for
-            automatic test.
+    return_pandas_data_series:          boolean, default = False
+        set to true to return Pandas Series as a result of transform.
+        setting to true reduces speed significantly but is required for
+        automatic test.
 
-        n_jobs:              int, optional, default = 1
-            The number of jobs to run in parallel for both `transform`.
-            ``-1`` means using all processors.
+    n_jobs:              int, optional, default = 1
+        The number of jobs to run in parallel for both `transform`.
+        ``-1`` means using all processors.
 
     Attributes
     ----------
@@ -96,7 +96,6 @@ class SFA(_PanelToPanelTransformer):
     num_insts = 0
     num_atts = 0
 
-
     References
     ----------
     .. [1] Schäfer, Patrick, and Mikael Högqvist. "SFA: a symbolic fourier approximation
@@ -104,7 +103,17 @@ class SFA(_PanelToPanelTransformer):
     15th international conference on extending database technology. 2012.
     """
 
-    _tags = {"univariate-only": True}
+    _tags = {
+        "univariate-only": True,
+        "scitype:transform-input": "Series",
+        # what is the scitype of X: Series, or Panel
+        "scitype:transform-output": "Series",
+        # what scitype is returned: Primitives, Series, Panel
+        "scitype:instancewise": False,  # is this an instance-wise transform?
+        "X_inner_mtype": "numpy3D",  # which mtypes do _fit/_predict support for X?
+        "y_inner_mtype": "pd_Series_Table",  # which mtypes does y require?
+        "requires_y": True,  # does y need to be passed in fit?
+    }
 
     def __init__(
         self,
@@ -160,7 +169,8 @@ class SFA(_PanelToPanelTransformer):
         self.skip_grams = skip_grams
 
         self.return_pandas_data_series = return_pandas_data_series
-        self.use_fallback_dft = (
+        self.use_fallback_dft = use_fallback_dft
+        self._use_fallback_dft = (
             use_fallback_dft if word_length < window_size - offset else True
         )
         self.typed_dict = typed_dict
@@ -178,6 +188,9 @@ class SFA(_PanelToPanelTransformer):
 
         super(SFA, self).__init__()
 
+        if not return_pandas_data_series:
+            self._output_convert = "off"
+
     def fit(self, X, y=None):
         """Calculate word breakpoints using MCB or IGB.
 
@@ -241,7 +254,7 @@ class SFA(_PanelToPanelTransformer):
         self._is_fitted = True
         return self
 
-    def transform(self, X, y=None):
+    def _transform(self, X, y=None):
         """Transform data into SFA words.
 
         Parameters
@@ -253,8 +266,6 @@ class SFA(_PanelToPanelTransformer):
         -------
         List of dictionaries containing SFA words
         """
-        self.check_is_fitted()
-        X = check_X(X, enforce_univariate=True, coerce_to_numpy=True)
         X = X.squeeze(1)
 
         with warnings.catch_warnings():
@@ -516,7 +527,7 @@ class SFA(_PanelToPanelTransformer):
                     self.inverse_sqrt_win_size,
                     self.lower_bounding,
                 )
-                if self.use_fallback_dft
+                if self._use_fallback_dft
                 else self._fast_fourier_transform(row)
             )
 
@@ -623,7 +634,7 @@ class SFA(_PanelToPanelTransformer):
         transformed = np.zeros((end, length))
 
         # first run with dft
-        if self.use_fallback_dft:
+        if self._use_fallback_dft:
             mft_data = self._discrete_fourier_transform(
                 series[0 : self.window_size],
                 self.dft_length,
diff --git a/sktime/transformations/panel/signature_based/_augmentations.py b/sktime/transformations/panel/signature_based/_augmentations.py
index d17052155..8a5f63c23 100644
--- a/sktime/transformations/panel/signature_based/_augmentations.py
+++ b/sktime/transformations/panel/signature_based/_augmentations.py
@@ -199,6 +199,7 @@ class _CumulativeSum(BaseTransformer):
 
     def __init__(self, append_zero=False):
         self.append_zero = append_zero
+        super(_CumulativeSum, self).__init__()
 
     def _transform(self, X, y=None):
         if self.append_zero:
diff --git a/sktime/transformations/series/boxcox.py b/sktime/transformations/series/boxcox.py
index daad5c38a..f0790f3cb 100644
--- a/sktime/transformations/series/boxcox.py
+++ b/sktime/transformations/series/boxcox.py
@@ -265,6 +265,7 @@ class LogTransformer(BaseTransformer):
         y : ignored argument for interface compatibility
             Additional data, e.g., labels for transformation
 
+
         Returns
         -------
         Xt : 2D np.ndarray
@@ -285,6 +286,7 @@ class LogTransformer(BaseTransformer):
         y : ignored argument for interface compatibility
             Additional data, e.g., labels for transformation
 
+
         Returns
         -------
         Xt : 2D np.ndarray
diff --git a/sktime/transformations/series/hidalgo.py b/sktime/transformations/series/hidalgo.py
new file mode 100644
index 000000000..85db4305b
--- /dev/null
+++ b/sktime/transformations/series/hidalgo.py
@@ -0,0 +1,762 @@
+# -*- coding: utf-8 -*-
+
+"""Hidalgo (Heterogeneous Intrinsic Dimensionality Algorithm) Segmentation."""
+
+__author__ = ["KatieBuc"]
+__all__ = ["Hidalgo"]
+
+
+from functools import reduce
+from typing import Union
+
+import numpy as np
+from sklearn.neighbors import NearestNeighbors
+from sklearn.utils.validation import check_random_state
+
+from sktime.transformations.base import BaseTransformer
+
+
+class Hidalgo(BaseTransformer):
+    """Heteregeneous Intrinsic Dimensionality Algorithm (Hidalgo) model.
+
+    Hidalgo is a robust approach in discriminating regions with
+    different local intrinsic dimensionality (topological feature
+    measuring complexity). Hidalgo offers unsupervised segmentation
+    of high-dimensional data.
+
+    Parameters
+    ----------
+    metric : str, or callable, optional, default="euclidean"
+        directly passed to sklearn KNearestNeighbors,
+        must be str or callable that can be passed to KNearestNeighbors
+        distance used in the nearest neighbors part of the algorithm
+    K : int, optional, default=2
+        number of manifolds used in algorithm
+    zeta : float, optional, defualt=0.8
+        "local homogeneity level" used in the algorithm, see equation (4)
+    q : int, optional, default=3
+        number of points for local Z interaction, "local homogeneity range"
+        see equation (4)
+    n_iter : int, optional, default=1000
+        number of Gibbs sampling iterations
+    n_replicas : int, optional, default=1
+        number of random starts to run Gibbs sampling
+    burn_in : float, optional, default=0.9
+        percentage of Gibbs sampling iterations discarded, "burn-in fraction"
+    fixed_Z : bool, optional, default=False
+        estimate parameters with fixed z (joint posterior approximation via Gibbs)
+        z = (z_1, ..., z_K) is a latent variable introduced, where z_i = k
+        indicates point i belongs to manifold K
+    use_Potts : bool, optional, default=True
+        if using local interaction between z, see equation (4)
+    estimate_zeta : bool, optional, default=False
+        update zeta in the sampling
+    sampling_rate: int, optional, default=10
+        rate at which to save samples for each n_iter
+    a : np.ArrayLike, optional, default=None
+        prior parameters of d, the dimensionality of manifold k
+    b : np.ArrayLike, optional, default=None
+        prior parameters of d, the dimensionality of manifold k
+    c : np.ArrayLike, optional, default=None
+        prior parameters of p, the probability that point belongs to manifold k
+    f : np.ArrayLike, optional, default=None
+        parameters of zeta
+    seed : int, optional, default = 1
+        generate random numbers with seed
+
+    References
+    ----------
+    Allegra, Michele, et al. "Data segmentation based on the local
+    intrinsic dimension." Scientific reports 10.1 (2020): 1-12.
+    https://www.nature.com/articles/s41598-020-72222-0
+
+
+    Examples
+    --------
+    >>> from sktime.transformations.series.hidalgo import Hidalgo
+    >>> import numpy as np
+    >>> np.random.seed(123)
+    >>> X = np.random.rand(10,3)
+    >>> X[:6, 1:] += 10
+    >>> X[6:, 1:] = 0
+    >>> model = Hidalgo(K=2, burn_in=0.8, n_iter=100, seed=10)
+    >>> fitted_model = model.fit(X)
+    >>> Z = fitted_model.transform(X)
+    >>> Z.tolist()
+    [1, 1, 1, 1, 1, 1, 0, 0, 0, 0]
+    """
+
+    _tags = {
+        "scitype:transform-input": "Series",
+        # what is the scitype of X: Series, or Panel
+        "scitype:transform-output": "Series",
+        # what scitype is returned: Primitives, Series, Panel
+        "transform-returns-same-time-index": True,
+        "univariate-only": False,
+        "fit_is_empty": False,
+    }
+
+    def __init__(
+        self,
+        metric="euclidean",
+        K=1,
+        zeta=0.8,
+        q=3,
+        n_iter=1000,
+        n_replicas=1,
+        burn_in=0.9,
+        fixed_Z=False,
+        use_Potts=True,
+        estimate_zeta=False,
+        sampling_rate=10,
+        a=None,
+        b=None,
+        c=None,
+        f=None,
+        seed=1,
+    ):
+
+        self.metric = metric
+        self.K = K
+        self.zeta = zeta
+        self.q = q
+        self.n_iter = n_iter
+        self.burn_in = burn_in
+        self.n_replicas = n_replicas
+        self.fixed_Z = fixed_Z
+        self.use_Potts = use_Potts
+        self.estimate_zeta = estimate_zeta
+        self.sampling_rate = sampling_rate
+        self.a = a
+        self.b = b
+        self.c = c
+        self.f = f
+        self.seed = seed
+
+        super(Hidalgo, self).__init__()
+
+    def _get_neighbourhood_params(self, X):
+        """
+        Neighbourhood information from input data X.
+
+        Parameters
+        ----------
+        X : 2D np.ndarray of shape (N, dim), where dim > 1
+            data to fit the algorithm to
+
+        Returns
+        -------
+        N : int
+            number of rows of X
+        mu : 1D np.ndarray of length N
+            paramerer in Pereto distribtion estimated by r2/r1
+        Iin : 1D np.ndarray of length N * q
+            encodes the q neighbour index values for point index i in 0:N-1
+            e.g. popint i=0 has neighbours 2, 4, 7 and point i=1
+            has neighbours 3, 9, 4
+            Iin = np.array([2, 4, 7, 3, 9, 4,...])
+        Iout : 1D np.ndarray of length N * q
+            array of indices for which i is neighbour for i in 0:N-1
+            e.g. point i=0 is also neighbour of points 2, 4 and
+            point i=1 is a neighbour of point 3 only
+            Iout = np.array([2, 4, 3,...])
+        Iout_count : 1D np.ndarray of length N
+            count of how many neighbours point i has for i in 0:N-1
+            e.g. Iout_count = np.array([2, 1, 1,...])
+        Iout_track : 1D np.ndarray of length N
+            cumulative sum of Iout_count at i-1 for i in 1:N-1
+            e.g. Iout_track = np.array([0, 2, 3, 4,...])
+        """
+        q = self.q
+        metric = self.metric
+
+        N, _ = np.shape(X)
+
+        nbrs = NearestNeighbors(
+            n_neighbors=q + 1, algorithm="ball_tree", metric=metric
+        ).fit(X)
+        distances, Iin = nbrs.kneighbors(X)
+        mu = np.divide(distances[:, 2], distances[:, 1])
+
+        nbrmat = np.zeros((N, N))
+        for n in range(q):
+            nbrmat[Iin[:, 0], Iin[:, n + 1]] = 1
+
+        Iout_count = np.sum(nbrmat, axis=0).astype(int)
+        Iout = np.where(nbrmat.T)[1].astype(int)
+        Iout_track = np.cumsum(Iout_count)
+        Iout_track = np.append(0, Iout_track[:-1]).astype(int)
+        Iin = Iin[:, 1:]
+        Iin = np.reshape(Iin, (N * q,)).astype(int)
+
+        return N, mu, Iin, Iout, Iout_count, Iout_track
+
+    def _update_zeta_prior(self, Z, N, Iin):
+        """Update prior parameters for zeta."""
+        q = self.q
+        f = self.f
+
+        if f is None:
+            f = np.ones(2)
+
+        N_in = sum([Z[Iin[q * i + j]] == Z[i] for j in range(q) for i in range(N)])
+
+        f1 = np.empty(shape=2)
+        f1[0] = f[0] + N_in
+        f1[1] = f[1] + N * q - N_in
+
+        return N_in, f1
+
+    def _initialise_params(self, N, mu, Iin, _rng):
+        """
+        Initialise parameters used in algorithm.
+
+        Outputs
+        ----------
+        V : 1D np.ndarray of length K
+            sum(log(mu_i)) for k in 0:K-1, when mu_i belongs to manifold k
+        NN : 1D np.ndarray of length K
+            count for k in 0:K-1, when data at index i belongs to manifold k
+        a1 : 1D np.ndarray of length K
+            prior parameters of d
+        b1 : 1D np.ndarray of length K
+            prior parameters of d
+        c1 : 1D np.ndarray of length K
+            prior parameters of p
+        Z : 1D np.ndarray of length N
+            segmentation based on manifold k
+        f1 : 1D np.ndarray of length 2
+            parameters of zeta
+        N_in : int
+            parameters of zeta
+        """
+        K = self.K
+        a = self.a
+        b = self.b
+        c = self.c
+        fixed_Z = self.fixed_Z
+
+        if a is None:
+            a = np.ones(K)
+        if b is None:
+            b = np.ones(K)
+        if c is None:
+            c = np.ones(K)
+
+        if not fixed_Z:
+            random_z = _rng.randint(0, K, N)
+            Z = np.array(random_z, dtype=int)
+        else:
+            Z = np.zeros(N, dtype=int)
+
+        V = [sum(np.log(mu[[Z[i] == k for i in range(N)]])) for k in range(K)]
+        NN = [sum([Z[i] == k for i in range(N)]) for k in range(K)]
+
+        a1 = a + NN
+        b1 = b + V
+        c1 = c + NN
+
+        N_in, f1 = self._update_zeta_prior(Z, N, Iin)
+
+        return (V, NN, a1, b1, c1, Z, f1, N_in)
+
+    def _gibbs_sampling(
+        self,
+        N,
+        mu,
+        Iin,
+        Iout,
+        Iout_count,
+        Iout_track,
+        V,
+        NN,
+        a1,
+        b1,
+        c1,
+        Z,
+        f1,
+        N_in,
+        _rng,
+    ):
+        """
+        Gibbs sampling method to find joint posterior distribution of target variables.
+
+        Notes
+        -----
+        Target parameters are d, p, Z
+        zeta must be computed for the probability distribution of the q-Neighbourhood
+        matrix
+
+        Parameters
+        ----------
+        V : 1D np.ndarray of length K
+            sum(log(mu_i)) for k in 0:K-1, when mu_i belongs to manifold k
+        NN : 1D np.ndarray of length K
+            count for k in 0:K-1, when data at index i belongs to manifold k
+        a1 : 1D np.ndarray of length K
+            prior parameters of d
+        b1 : 1D np.ndarray of length K
+            prior parameters of d
+        c1 : 1D np.ndarray of length K
+            prior parameters of p
+        Z : 1D np.ndarray of length N
+            segmentation based on manifold k
+        f1 : 1D np.ndarray of length 2
+            parameters of zeta
+        N_in : int
+            parameters of zeta
+
+        Returns
+        -------
+        sampling : 2D np.ndarray of shape (n_iter, Npar), where Npar = N + 2 * K + 2 + 1
+            posterior samples of d, p, Z and likelihood samples, respectively.
+
+        """
+        zeta = self.zeta
+        q = self.q
+        K = self.K
+        n_iter = self.n_iter
+        fixed_Z = self.fixed_Z
+        use_Potts = self.use_Potts
+        estimate_zeta = self.estimate_zeta
+
+        sampling = np.empty(shape=0)
+        pp = (K - 1) / K
+        p = np.ones(shape=K) / K
+
+        def sample_d(K, a1, b1, _rng):
+            """Sample d, the dimension of manifold k."""
+            d = np.empty(shape=K)
+            for k in range(K):
+                stop = False
+
+                while stop is False:
+                    r1 = _rng.random() * 200  # random sample for d[k]
+                    r2 = _rng.random()  # random number for accepting
+
+                    rmax = (a1[k] - 1) / b1[k]
+
+                    if a1[k] - 1 > 0:
+                        assert rmax > 0
+                        frac = np.exp(
+                            -b1[k] * (r1 - rmax)
+                            - (a1[k] - 1) * (np.log(rmax) - np.log(r1))
+                        )
+                    else:
+                        frac = np.exp(-b1[k] * r1)
+
+                    if frac > r2:
+                        stop = True
+                        d[k] = r1
+
+            return d
+
+        def sample_p(K, p, pp, c1, _rng):
+            """Sample p, the prior probability that a point belongs to manifold k."""
+            for k in range(K - 1):
+                stop = False
+
+                while stop is False:
+                    r1 = _rng.random()  # random sample for p[k]
+                    r2 = _rng.random()  # random number for accepting
+
+                    rmax = (c1[k] - 1) / (c1[k] - 1 + c1[K - 1] - 1)
+                    frac = ((r1 / rmax) ** (c1[k] - 1)) * (
+                        ((1 - r1) / (1 - rmax)) ** (c1[K - 1] - 1)
+                    )
+
+                    if frac > r2:
+                        stop = True
+                        r1 = r1 * (1.0 - pp + p[k])
+                        p[K - 1] += p[k] - r1
+                        pp -= p[k] - r1
+                        p[k] = r1
+
+            return (p, pp)
+
+        def sample_zeta(N, K, zeta, use_Potts, estimate_zeta, q, NN, f1, it, _rng):
+            """Sample zeta, sample probability neighbour of point from same manifold."""
+            stop = False
+            maxval = -100000
+
+            if use_Potts and estimate_zeta:
+                for zeta_candidates in range(10):
+                    zeta1 = 0.5 + 0.05 * zeta_candidates
+                    ZZ = [partition_function(N, NN[k], zeta1, q) for k in range(K)]
+                    h = [NN[k] * np.log(ZZ[k]) for k in range(K)]
+                    val = (
+                        (f1[0] - 1) * np.log(zeta1)
+                        + (f1[1] - 1) * np.log(1 - zeta1)
+                        - h
+                    )
+
+                    if val > maxval:
+                        maxval = val
+
+                while stop is False:
+                    r1 = _rng.random()  # random sample for zeta
+                    r2 = _rng.random()  # random number for accepting
+
+                    ZZ = [partition_function(N, NN[k], r1, q) for k in range(K)]
+                    h = [NN[k] * np.log(ZZ[k]) for k in range(K)]
+                    val = (f1[0] - 1) * np.log(r1) + (f1[1] - 1) * np.log(1 - r1) - h
+                    frac = np.exp(val - maxval)
+
+                    if frac > r2:
+                        stop = True
+                        if it > 0:
+                            zeta = r1
+
+            return zeta
+
+        def sample_Z(
+            N,
+            mu,
+            Iin,
+            Iout,
+            Iout_count,
+            Iout_track,
+            Z,
+            NN,
+            a1,
+            c1,
+            V,
+            b1,
+            zeta,
+            fixed_Z,
+            q,
+            _rng,
+        ):
+            """Sample z, latent variable indicating point i belongs to manifold k."""
+            if (abs(zeta - 1) < 1e-5) or fixed_Z:
+                return Z, NN, a1, c1, V, b1
+
+            for i in range(N):
+                stop = False
+                gg = np.empty(shape=K)
+                gmax = 0
+
+                for k1 in range(K):
+                    g = 0
+                    if use_Potts:
+
+                        n_in = sum([Z[Iin[q * i + j]] == k1 for j in range(q)])
+
+                        m_in = sum(
+                            [
+                                Iout[Iout_track[i] + j] > -1
+                                and Z[Iout[Iout_track[i] + j]] == k1
+                                for j in range(Iout_count[i])
+                            ]
+                        )
+
+                        g = (n_in + m_in) * np.log(zeta / (1 - zeta)) - np.log(
+                            partition_function(N, NN[k1], zeta, q)
+                        )
+                        var = partition_function(
+                            N, NN[k1] - 1, zeta, q
+                        ) / partition_function(N, NN[k1], zeta, q)
+                        assert var > 0
+                        g = g + np.log(var) * (NN[k1] - 1)
+
+                    if g > gmax:
+                        gmax = g
+                    gg[k1] = g
+
+                gg = [np.exp(gg[k1] - gmax) for k1 in range(K)]
+
+                prob = p * d * mu[i] ** (-(d + 1)) * gg
+                prob /= prob.sum()
+
+                while stop is False:
+                    r1 = int(np.floor(_rng.random() * K))  # random sample for Z
+                    r2 = _rng.random()  # random number for accepting
+
+                    if prob[r1] > r2:
+                        stop = True
+                        # minus values
+                        NN[Z[i]] -= 1
+                        a1[Z[i]] -= 1
+                        c1[Z[i]] -= 1
+                        V[Z[i]] -= np.log(mu[i])
+                        b1[Z[i]] -= np.log(mu[i])
+                        # change Z, add values
+                        Z[i] = r1
+                        NN[Z[i]] += 1
+                        a1[Z[i]] += 1
+                        c1[Z[i]] += 1
+                        V[Z[i]] += np.log(mu[i])
+                        b1[Z[i]] += np.log(mu[i])
+
+            return Z, NN, a1, c1, V, b1
+
+        def sample_likelihood(N, mu, p, d, Z, N_in, zeta, NN):
+            """Sample likelihood values of mu, and local inhomogeneity penalisation."""
+            lik0 = 0
+            for i in range(N):
+                lik0 = (
+                    lik0
+                    + np.log(p[Z[i]])
+                    + np.log(d[Z[i]])
+                    - (d[Z[i]] + 1) * np.log(mu[i])
+                )
+
+            lik1 = lik0 + np.log(zeta / (1 - zeta)) * N_in
+
+            for k1 in range(K):
+                lik1 = lik1 - (NN[k1] * np.log(partition_function(N, NN[k1], zeta, q)))
+
+            return lik0, lik1
+
+        for it in range(n_iter):
+
+            d = sample_d(K, a1, b1, _rng)
+            sampling = np.append(sampling, d)
+
+            (p, pp) = sample_p(K, p, pp, c1, _rng)
+            sampling = np.append(sampling, p[: K - 1])
+            sampling = np.append(sampling, (1 - pp))
+
+            zeta = sample_zeta(
+                N, K, zeta, use_Potts, estimate_zeta, q, NN, f1, it, _rng
+            )
+            sampling = np.append(sampling, zeta)
+
+            Z, NN, a1, c1, V, b1 = sample_Z(
+                N,
+                mu,
+                Iin,
+                Iout,
+                Iout_count,
+                Iout_track,
+                Z,
+                NN,
+                a1,
+                c1,
+                V,
+                b1,
+                zeta,
+                fixed_Z,
+                q,
+                _rng,
+            )
+            sampling = np.append(sampling, Z)
+
+            N_in, f1 = self._update_zeta_prior(Z, N, Iin)
+
+            lik = sample_likelihood(N, mu, p, d, Z, N_in, zeta, NN)
+            sampling = np.append(sampling, lik)
+
+        return sampling
+
+    def _fit(self, X, y=None):
+        """Run the Hidalgo algorithm.
+
+        Find parameter esimates as distributions in sampling.
+        Iterate through n_replicas random starts and get posterior
+        samples with best max likelihood.
+
+        Notes
+        -----
+        Writes to self
+        _d : 1D np.ndarray of length K
+            posterior mean of d, from posterior sample in gibbs_sampling
+        _derr : 1D np.ndarray of length K
+            posterior std of d, from posterior sample in gibbs_sampling
+        _p : 1D np.ndarray of length K
+            posterior mean of p, from posterior sample in gibbs_sampling
+        _perr : 1D np.ndarray of length K
+            posterior std of p, from posterior sample in gibbs_sampling
+        _lik : float
+            mean of likelihood, from sample in gibbs_sampling
+        _likerr : float
+            std of likelihood, from sample in gibbs_sampling
+        _Pi : 2D np.ndarray of shape (K, N)
+            probability of posterior of z_i = k, point i can be safely
+            assigned to manifold k if Pi > 0.8
+        _Z : 1D np.ndarray of length N
+            base-zero integer values corresponsing to segment (manifold k)
+
+        Parameters
+        ----------
+        X : 2D np.ndarray of shape (N, dim), where dim > 1
+            data to fit the algorithm to
+
+        Returns
+        -------
+        self
+        """
+        K = self.K
+        n_replicas = self.n_replicas
+        n_iter = self.n_iter
+        sampling_rate = self.sampling_rate
+        burn_in = self.burn_in
+        seed = self.seed
+
+        _rng = check_random_state(seed)
+
+        N, mu, Iin, Iout, Iout_count, Iout_track = self._get_neighbourhood_params(X)
+        V, NN, a1, b1, c1, Z, f1, N_in = self._initialise_params(N, mu, Iin, _rng)
+
+        Npar = N + 2 * K + 2 + 1
+        bestsampling = np.zeros(shape=0)
+        maxlik = -1e10
+
+        for _ in range(n_replicas):
+
+            sampling = self._gibbs_sampling(
+                N,
+                mu,
+                Iin,
+                Iout,
+                Iout_count,
+                Iout_track,
+                V,
+                NN,
+                a1,
+                b1,
+                c1,
+                Z,
+                f1,
+                N_in,
+                _rng,
+            )
+            sampling = np.reshape(sampling, (n_iter, Npar))
+
+            idx = [
+                it
+                for it in range(n_iter)
+                if it % sampling_rate == 0 and it >= n_iter * burn_in
+            ]
+            sampling = sampling[
+                idx,
+            ]
+
+            lik = np.mean(sampling[:, -1], axis=0)
+
+            if lik > maxlik:
+                bestsampling = sampling
+                maxlik = lik
+
+        self._d = np.mean(bestsampling[:, :K], axis=0)
+        self._derr = np.std(bestsampling[:, :K], axis=0)
+        self._p = np.mean(bestsampling[:, K : 2 * K], axis=0)
+        self._perr = np.std(bestsampling[:, K : 2 * K], axis=0)
+        self._lik = np.mean(bestsampling[:, -1], axis=0)
+        self._likerr = np.std(bestsampling[:, -1], axis=0)
+
+        Pi = np.zeros((K, N))
+
+        for k in range(K):
+            Pi[k, :] = np.sum(bestsampling[:, (2 * K) + 1 : 2 * K + N + 1] == k, axis=0)
+
+        Pi = Pi / len(idx)
+        self._Pi = Pi
+
+        Z = np.argmax(Pi, axis=0)
+        pZ = np.max(Pi, axis=0)
+        Z[np.where(pZ < 0.8)] = -1
+        self._Z = Z
+
+        return self
+
+    def _transform(self, X, y=None):
+        """Transform X and return a transformed version.
+
+        private _transform containing core logic, called from transform
+
+        Parameters
+        ----------
+        X : Series of mtype X_inner_mtype
+            if X_inner_mtype is list, _transform must support all types in it
+            Data to be transformed
+        y : Series of mtype y_inner_mtype, default=None
+            Not used in this unsupervised implementation
+
+        Returns
+        -------
+        transformed version of X
+        """
+        return self._Z
+
+    @classmethod
+    def get_test_params(cls, parameter_set="default"):
+        """Return testing parameter settings for the estimator.
+
+        Parameters
+        ----------
+        parameter_set : str, default="default"
+            Name of the set of test parameters to return, for use in tests. If no
+            special parameters are defined for a value, will return `"default"` set.
+            Reserved values for classifiers:
+                "results_comparison" - used for identity testing in some classifiers
+                    should contain parameter settings comparable to "TSC bakeoff"
+
+        Returns
+        -------
+        params : dict or list of dict, default = {}
+            Parameters to create testing instances of the class
+            Each dict are parameters to construct an "interesting" test instance, i.e.,
+            `MyClass(**params)` or `MyClass(**params[i])` creates a valid test instance.
+            `create_test_instance` uses the first (or only) dictionary in `params`
+        """
+        return {
+            "metric": "euclidean",
+            "K": 1,
+            "zeta": 0.8,
+            "q": 3,
+            "n_iter": 10,
+            "n_replicas": 1,
+            "burn_in": 0.5,
+            "fixed_Z": False,
+            "use_Potts": True,
+            "estimate_zeta": False,
+            "sampling_rate": 2,
+            "a": None,
+            "b": None,
+            "c": None,
+            "f": None,
+            "seed": 1,
+        }
+
+
+def binom(N: Union[int, float], q: Union[int, float]):
+    """Calculate the binomial coefficient.
+
+    Parameters
+    ----------
+    N : int, float
+        number of fixed elements from qhich q is chosen
+    q : int, float
+        number of subset q elements chosen from N
+    """
+    if q == 0:
+        return 1.0
+    if N < 0:
+        return 0.0
+    return reduce(lambda x, y: x * y, [(N - q1) / (q1 + 1) for q1 in range(q)])
+
+
+def partition_function(N, N1, zeta, q):
+    """Partition function for Z.
+
+    Parameters
+    ----------
+    N : int, float
+        number of rows of input data X
+    N1 : int, float
+        parameter value from NN[k] for k=0:K-1
+    zeta : float
+        parameter value zeta
+    q : int, float
+        parameter value q
+    """
+    return sum(
+        [
+            binom(N1 - 1, q1)
+            * binom(N - N1, q - q1)
+            * zeta ** (q1)
+            * (1 - zeta) ** (q - q1)
+            for q1 in range(q + 1)
+        ]
+    )


Test patch:
diff --git a/sktime/annotation/tests/test_ggs.py b/sktime/annotation/tests/test_ggs.py
new file mode 100644
index 000000000..6f09abe50
--- /dev/null
+++ b/sktime/annotation/tests/test_ggs.py
@@ -0,0 +1,34 @@
+# -*- coding: utf-8 -*-
+"""Tests for GGS module."""
+
+import numpy as np
+import pytest
+
+from sktime.annotation.ggs import GGS, GreedyGaussianSegmentation
+
+
+@pytest.fixture
+def univariate_mean_shift():
+    """Generate simple mean shift time series."""
+    x = np.concatenate(tuple(np.ones(5) * i**2 for i in range(4)))
+    return x[:, np.newaxis]
+
+
+def test_GGS_find_change_points(univariate_mean_shift):
+    """Test the GGS core estimator."""
+    ggs = GGS(k_max=10, lamb=1.0)
+    pred = ggs.find_change_points(univariate_mean_shift)
+    assert isinstance(pred, list)
+    assert len(pred) == 5
+
+
+def test_GreedyGaussianSegmentation(univariate_mean_shift):
+    """Test the GreedyGaussianSegmentation."""
+    ggs = GreedyGaussianSegmentation(k_max=5, lamb=0.5)
+    assert ggs.get_params() == {
+        "k_max": 5,
+        "lamb": 0.5,
+        "verbose": False,
+        "max_shuffles": 250,
+        "random_state": None,
+    }
diff --git a/sktime/datasets/tests/test_data_io.py b/sktime/datasets/tests/test_data_io.py
index 7e74c0546..bfbc969c0 100644
--- a/sktime/datasets/tests/test_data_io.py
+++ b/sktime/datasets/tests/test_data_io.py
@@ -31,6 +31,7 @@ from sktime.datasets._data_io import (
     _load_provided_dataset,
 )
 from sktime.datatypes import MTYPE_LIST_PANEL, check_is_mtype
+from sktime.utils.validation._dependencies import _check_soft_dependencies
 
 # Disabling test for these mtypes since they don't support certain functionality yet
 _TO_DISABLE = ["pd-long", "pd-wide", "numpyflat"]
@@ -1406,6 +1407,10 @@ def test_convert_tsf_to_multiindex(freq):
     )
 
 
+@pytest.mark.skipif(
+    not _check_soft_dependencies("backoff", severity="none"),
+    reason="load_solar requires backoff in the environment",
+)
 @pytest.mark.parametrize("return_df", [False, True])
 def test_load_solar(return_df):
     """Test function for loading solar data through the Sheffiled Solar API."""
diff --git a/sktime/tests/_config.py b/sktime/tests/_config.py
index facb82162..2d9c41628 100644
--- a/sktime/tests/_config.py
+++ b/sktime/tests/_config.py
@@ -92,6 +92,10 @@ EXCLUDED_TESTS = {
         "test_update_predict_single",  # see 2997, sporadic failure, unknown cause
         "test__y_when_refitting",  # see 3176
     ],
+    # GGS inherits from BaseEstimator which breaks this test
+    "GreedyGaussianSegmentation": ["test_inheritance", "test_create_test_instance"],
+    "SAX": "test_fit_transform_output",  # SAX returns strange output format
+    # this needs to be fixed, was not tested previously due to legacy exception
 }
 
 # We use estimator tags in addition to class hierarchies to further distinguish
diff --git a/sktime/transformations/series/tests/test_hidalgo.py b/sktime/transformations/series/tests/test_hidalgo.py
new file mode 100644
index 000000000..e08780398
--- /dev/null
+++ b/sktime/transformations/series/tests/test_hidalgo.py
@@ -0,0 +1,288 @@
+# -*- coding: utf-8 -*-
+"""Test for hidalgo segmentation."""
+import numpy as np
+from sklearn.utils.validation import check_random_state
+
+from sktime.transformations.series.hidalgo import Hidalgo
+
+# generate dataset
+np.random.seed(10002)
+X = np.zeros((10, 6))
+
+# half the points from one generating regime
+for j in range(1):
+    X[:5, j] = np.random.normal(0, 3, 5)
+
+# the other half from another
+for j in range(3):
+    X[5:, j] = np.random.normal(2, 1, 5)
+
+
+def test_X():
+    """Test if innput data is of expected dimension and type."""
+    assert isinstance(X, np.ndarray), "X should be a numpy array"
+    assert len(np.shape(X)) == 2, "X should be a two-dimensional numpy array"
+
+
+def test_get_neighbourhood_params():
+    """Test for neighbourhood parameter generation."""
+    model = Hidalgo(K=2, n_iter=10, burn_in=0.5, sampling_rate=2, seed=1)
+    (
+        N_actual,
+        mu_actual,
+        Iin_actual,
+        Iout_actual,
+        Iout_count_actual,
+        Iout_track_actual,
+    ) = model._get_neighbourhood_params(X)
+
+    N_expected = 10
+    mu_expected = np.array(
+        [
+            1.464722,
+            5.40974013,
+            2.464722,
+            4.78175011,
+            1.17460641,
+            1.95890874,
+            1.05863706,
+            1.06313386,
+            1.85040636,
+            1.28065266,
+        ]
+    )
+    Iin_expected = np.array(
+        [
+            2,
+            4,
+            7,
+            3,
+            9,
+            4,
+            0,
+            4,
+            7,
+            1,
+            4,
+            9,
+            0,
+            9,
+            3,
+            8,
+            6,
+            9,
+            8,
+            5,
+            9,
+            5,
+            8,
+            9,
+            5,
+            6,
+            7,
+            5,
+            7,
+            4,
+        ]
+    )
+    Iout_expected = np.array(
+        [
+            2,
+            4,
+            3,
+            0,
+            1,
+            4,
+            0,
+            1,
+            2,
+            3,
+            9,
+            6,
+            7,
+            8,
+            9,
+            5,
+            8,
+            0,
+            2,
+            8,
+            9,
+            5,
+            6,
+            7,
+            1,
+            3,
+            4,
+            5,
+            6,
+            7,
+        ]
+    )
+    Iout_count_expected = np.array([2, 1, 1, 2, 5, 4, 2, 4, 3, 6])
+    Iout_track_expected = np.array([0, 2, 3, 4, 6, 11, 15, 17, 21, 24])
+
+    assert np.allclose(N_actual, N_expected)
+    assert np.allclose(mu_actual, mu_expected)
+    assert np.allclose(Iin_actual, Iin_expected)
+    assert np.allclose(Iout_actual, Iout_expected)
+    assert np.allclose(Iout_count_actual, Iout_count_expected)
+    assert np.allclose(Iout_track_actual, Iout_track_expected)
+
+
+def test_initialise_params():
+    """Test for initialise parameters."""
+    model = Hidalgo(K=2, n_iter=10, burn_in=0.5, sampling_rate=2, seed=1)
+    N = 10
+    mu = np.array(
+        [
+            1.464722,
+            5.40974013,
+            2.464722,
+            4.78175011,
+            1.17460641,
+            1.95890874,
+            1.05863706,
+            1.06313386,
+            1.85040636,
+            1.28065266,
+        ]
+    )
+    Iin = np.array(
+        [
+            2,
+            4,
+            7,
+            3,
+            9,
+            4,
+            0,
+            4,
+            7,
+            1,
+            4,
+            9,
+            0,
+            9,
+            3,
+            8,
+            6,
+            9,
+            8,
+            5,
+            9,
+            5,
+            8,
+            9,
+            5,
+            6,
+            7,
+            5,
+            7,
+            4,
+        ]
+    )
+    _rng = check_random_state(model.seed)
+
+    (
+        V_actual,
+        NN_actual,
+        a1_actual,
+        b1_actual,
+        c1_actual,
+        Z_actual,
+        f1_actual,
+        N_in_actual,
+    ) = model._initialise_params(N, mu, Iin, _rng)
+
+    V_expected = [2.7142554736925324, 3.6367957651001124]
+    NN_expected = [3, 7]
+    a1_expected = np.array([4.0, 8.0])
+    b1_expected = np.array([3.71425547, 4.63679577])
+    c1_expected = np.array([4.0, 8.0])
+    Z_expected = np.array([1, 1, 0, 0, 1, 1, 1, 1, 1, 0])
+    f1_expected = np.array([15.0, 17.0])
+    N_in_expected = 14
+
+    assert np.allclose(V_actual, V_expected)
+    assert np.allclose(NN_actual, NN_expected)
+    assert np.allclose(a1_actual, a1_expected)
+    assert np.allclose(b1_actual, b1_expected)
+    assert np.allclose(c1_actual, c1_expected)
+    assert np.allclose(Z_actual, Z_expected)
+    assert np.allclose(f1_actual, f1_expected)
+    assert N_in_actual == N_in_expected
+
+
+def test_gibbs():
+    """Test _transform method including gibbs sampling iterations."""
+    expected = [
+        [
+            0.72269469,
+            3.20993546,
+            0.59599049,
+            0.40400951,
+            0.8,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            1.0,
+            1.0,
+            1.0,
+            1.0,
+            1.0,
+            -17.9650666,
+            -4.50125762,
+        ],
+        [
+            0.85025968,
+            2.32993642,
+            0.80166206,
+            0.19833794,
+            0.8,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            0.0,
+            1.0,
+            1.0,
+            1.0,
+            1.0,
+            0.0,
+            -19.21897121,
+            -7.78773922,
+        ],
+    ]
+
+    model = Hidalgo(K=2, n_iter=10, burn_in=0.5, sampling_rate=2, seed=1)
+    _rng = check_random_state(model.seed)
+
+    N, mu, Iin, Iout, Iout_count, Iout_track = model._get_neighbourhood_params(X)
+    V, NN, a1, b1, c1, Z, f1, N_in = model._initialise_params(N, mu, Iin, _rng)
+
+    sampling = model._gibbs_sampling(
+        N,
+        mu,
+        Iin,
+        Iout,
+        Iout_count,
+        Iout_track,
+        V,
+        NN,
+        a1,
+        b1,
+        c1,
+        Z,
+        f1,
+        N_in,
+        _rng,
+    )
+    sampling = np.reshape(sampling, (10, 17))
+    actual = sampling[
+        [6, 8],
+    ]
+
+    assert np.allclose(actual, expected)
diff --git a/sktime/utils/_testing/deep_equals.py b/sktime/utils/_testing/deep_equals.py
index 8a7b7d4cb..bd6fa1e5c 100644
--- a/sktime/utils/_testing/deep_equals.py
+++ b/sktime/utils/_testing/deep_equals.py
@@ -118,9 +118,10 @@ def deep_equals(x, y, return_msg=False):
         )
     elif type(x).__name__ == "ForecastingHorizon":
         return ret(*_fh_equals(x, y, return_msg=True))
-    elif x != y:
+    elif isinstance(x != y, bool) and x != y:
+        return ret(False, f" !=, {x} != {y}")
+    elif np.any(x != y):
         return ret(False, f" !=, {x} != {y}")
-
     return ret(True, "")