Repository: alan-turing-institute/sktime
Pull request #2964 (2022-07-10 02:15:42)

[ENH] Msm distance cleanup
<!--
Thanks for contributing a pull request! Please ensure you have taken a look
at our contribution guide: https://github.com/alan-turing-institute/sktime/blob/main/CONTRIBUTING.md
-->
This PR cleans up msm redundent functions and improves docstrong

#### Reference Issues/PRs
<!--
Example: Fixes #1234. See also #3456.

Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->


#### What does this implement/fix? Explain your changes.
<!--
A clear and concise description of what you have implemented.
-->
MSM was missing docstrings for variables and additionally there was redundent functions in the file which have been removed.
Additionally code has been simplified using built in python functions.

#### Does your contribution introduce a new dependency? If yes, which one?

<!--
If your contribution does add a new hard dependency, we may suggest to initially develop your contribution in a separate companion package in https://github.com/sktime/ to keep external dependencies of the core sktime package to a minimum.
-->

#### What should a reviewer concentrate their feedback on?

<!-- This section is particularly useful if you have a pull request that is still in development. You can guide the reviews to focus on the parts that are ready for their comments. We suggest using bullets (indicated by * or -) and filled checkboxes [x] here -->

#### Any other comments?
<!--
Please be aware that we are a loose team of volunteers so patience is necessary; assistance handling other issues is very welcome. We value all user contributions, no matter how minor they are. If we are slow to review, either the pull request needs some benchmarking, tinkering, convincing, etc. or more likely the reviewers are simply busy. In either case, we ask for your understanding during the review process.
-->

#### PR checklist
<!--
Please go through the checklist below. Please feel free to remove points if they are not applicable.
-->

##### For all contributions
- [x] I've added myself to the [list of contributors](https://github.com/alan-turing-institute/sktime/blob/main/.all-contributorsrc).
- [x] Optionally, I've updated sktime's [CODEOWNERS](https://github.com/alan-turing-institute/sktime/blob/main/CODEOWNERS) to receive notifications about future changes to these files.
- [x] I've added unit tests and made sure they pass locally.
- [x] The PR title starts with either [ENH], [MNT], [DOC], or [BUG] indicating whether the PR topic is related to enhancement, maintenance, documentation, or bug.

##### For new estimators
- [x] I've added the estimator to the online documentation.
- [x] I've updated the existing example notebooks or provided a new one to showcase how my estimator works.


<!--
Thanks for contributing!
-->



Patch:
diff --git a/sktime/distances/_distance.py b/sktime/distances/_distance.py
index e0b647335..f4724dae0 100644
--- a/sktime/distances/_distance.py
+++ b/sktime/distances/_distance.py
@@ -694,8 +694,11 @@ def dtw_distance(
 def msm_distance(
     x: np.ndarray,
     y: np.ndarray,
-    c: float = 0.0,
-    **kwargs: Any,
+    c: float = 1.0,
+    window: float = None,
+    itakura_max_slope: float = None,
+    bounding_matrix: np.ndarray = None,
+    **kwargs: dict,
 ) -> float:
     """Compute the move-split-merge distance.
 
@@ -713,8 +716,21 @@ def msm_distance(
         First time series.
     y: np.ndarray (1d or 2d array)
         Second time series.
-    kwargs: Any
-        Extra kwargs.
+    c: float, default = 1.0
+        Cost for split or merge operation.
+    window: Float, defaults = None
+        Float that is the radius of the sakoe chiba window (if using Sakoe-Chiba
+        lower bounding). Must be between 0 and 1.
+    itakura_max_slope: float, defaults = None
+        Gradient of the slope for itakura parallelogram (if using Itakura
+        Parallelogram lower bounding). Must be between 0 and 1.
+    bounding_matrix: np.ndarray (2d array of shape (m1,m2)), defaults = None
+        Custom bounding matrix to use. If defined then other lower_bounding params
+        are ignored. The matrix should be structure so that indexes considered in
+        bound should be the value 0. and indexes outside the bounding matrix should
+        be infinity.
+    kwargs: any
+        extra kwargs.
 
     Returns
     -------
@@ -731,7 +747,6 @@ def msm_distance(
         NumbaDistance.
         If a resolved metric is not no_python compiled.
         If the metric type cannot be determined
-
     References
     ----------
     .. [1]A.  Stefan,  V.  Athitsos,  and  G.  Das.   The  Move-Split-Merge  metric
@@ -740,6 +755,9 @@ def msm_distance(
     """
     format_kwargs = {
         "c": c,
+        "window": window,
+        "itakura_max_slope": itakura_max_slope,
+        "bounding_matrix": bounding_matrix,
     }
     format_kwargs = {**format_kwargs, **kwargs}
 
@@ -1601,8 +1619,11 @@ def msm_alignment_path(
     x: np.ndarray,
     y: np.ndarray,
     return_cost_matrix: bool = False,
-    c: float = 0.0,
-    **kwargs: Any,
+    c: float = 1.0,
+    window: float = None,
+    itakura_max_slope: float = None,
+    bounding_matrix: np.ndarray = None,
+    **kwargs: dict,
 ) -> AlignmentPathReturn:
     """Compute the move-split-merge alignment path.
 
@@ -1622,8 +1643,21 @@ def msm_alignment_path(
         Second time series.
     return_cost_matrix: bool, defaults = False
         Boolean that when true will also return the cost matrix.
-    kwargs: Any
-        Extra kwargs.
+    c: float, default = 1.0
+        Cost for split or merge operation.
+    window: float, defaults = None
+        Float that is the radius of the sakoe chiba window (if using Sakoe-Chiba
+        lower bounding). Must be between 0 and 1.
+    itakura_max_slope: float, defaults = None
+        Gradient of the slope for itakura parallelogram (if using Itakura
+        Parallelogram lower bounding). Must be between 0 and 1.
+    bounding_matrix: np.ndarray (2d array of shape (m1,m2)), defaults = None
+        Custom bounding matrix to use. If defined then other lower_bounding params
+        are ignored. The matrix should be structure so that indexes considered in
+        bound should be the value 0. and indexes outside the bounding matrix should
+        be infinity.
+    kwargs: any
+        extra kwargs.
 
     Returns
     -------
@@ -1645,7 +1679,6 @@ def msm_alignment_path(
         NumbaDistance.
         If a resolved metric is not no_python compiled.
         If the metric type cannot be determined
-
     References
     ----------
     .. [1]A.  Stefan,  V.  Athitsos,  and  G.  Das.   The  Move-Split-Merge  metric
@@ -1654,6 +1687,9 @@ def msm_alignment_path(
     """
     format_kwargs = {
         "c": c,
+        "window": window,
+        "itakura_max_slope": itakura_max_slope,
+        "bounding_matrix": bounding_matrix,
     }
     format_kwargs = {**format_kwargs, **kwargs}
 
diff --git a/sktime/distances/_msm.py b/sktime/distances/_msm.py
index b68e3f63c..d15126268 100644
--- a/sktime/distances/_msm.py
+++ b/sktime/distances/_msm.py
@@ -49,7 +49,7 @@ class _MsmDistance(NumbaDistance):
         x: np.ndarray,
         y: np.ndarray,
         return_cost_matrix: bool = False,
-        c: float = 0.0,
+        c: float = 1.0,
         window: float = None,
         itakura_max_slope: float = None,
         bounding_matrix: np.ndarray = None,
@@ -58,19 +58,31 @@ class _MsmDistance(NumbaDistance):
         """Create a no_python compiled MSM distance path callable.
 
         Series should be shape (1, m), where m is the series length. Series can be
-        different
-        lengths.
+        different lengths.
 
         Parameters
         ----------
-        x: np.ndarray (2d array of shape (1,m1)).
+        x: np.ndarray (2d array of shape (d,m1)).
             First time series.
-        y: np.ndarray (2d array of shape (1,m2)).
+        y: np.ndarray (2d array of shape (d,m2)).
             Second time series.
         return_cost_matrix: bool, defaults = False
             Boolean that when true will also return the cost matrix.
-        c: float
-            parameter used in MSM (update later!)
+        c: float, default = 1.0
+            Cost for split or merge operation.
+        window: float, defaults = None
+            Float that is the radius of the sakoe chiba window (if using Sakoe-Chiba
+            lower bounding). Must be between 0 and 1.
+        itakura_max_slope: float, defaults = None
+            Gradient of the slope for itakura parallelogram (if using Itakura
+            Parallelogram lower bounding). Must be between 0 and 1.
+        bounding_matrix: np.ndarray (2d array of shape (m1,m2)), defaults = None
+            Custom bounding matrix to use. If defined then other lower_bounding params
+            are ignored. The matrix should be structure so that indexes considered in
+            bound should be the value 0. and indexes outside the bounding matrix should
+            be infinity.
+        kwargs: any
+            extra kwargs.
 
         Returns
         -------
@@ -93,6 +105,7 @@ class _MsmDistance(NumbaDistance):
                 f"univariate series, passed seris shape {x.shape[0]} and"
                 f"shape {y.shape[0]}"
             )
+
         _bounding_matrix = resolve_bounding_matrix(
             x, y, window, itakura_max_slope, bounding_matrix
         )
@@ -104,7 +117,7 @@ class _MsmDistance(NumbaDistance):
                 _x: np.ndarray,
                 _y: np.ndarray,
             ) -> Tuple[List, float, np.ndarray]:
-                cost_matrix = _cost_matrix(_x, _y, c, _bounding_matrix)
+                cost_matrix = _cost_matrix(_x, _y, _bounding_matrix, c)
                 path = compute_min_return_path(cost_matrix, _bounding_matrix)
                 return path, cost_matrix[-1, -1], cost_matrix
 
@@ -115,7 +128,7 @@ class _MsmDistance(NumbaDistance):
                 _x: np.ndarray,
                 _y: np.ndarray,
             ) -> Tuple[List, float]:
-                cost_matrix = _cost_matrix(_x, _y, c, _bounding_matrix)
+                cost_matrix = _cost_matrix(_x, _y, _bounding_matrix, c)
                 path = compute_min_return_path(cost_matrix, _bounding_matrix)
                 return path, cost_matrix[-1, -1]
 
@@ -125,7 +138,7 @@ class _MsmDistance(NumbaDistance):
         self,
         x: np.ndarray,
         y: np.ndarray,
-        c: float = 0.0,
+        c: float = 1.0,
         window: float = None,
         itakura_max_slope: float = None,
         bounding_matrix: np.ndarray = None,
@@ -133,18 +146,30 @@ class _MsmDistance(NumbaDistance):
     ) -> DistanceCallable:
         """Create a no_python compiled MSM distance callable.
 
-        Series should be shape (1, m), where m is the series length. Series can be
-        different
-        lengths.
+        Series should be shape (d, m) where m is the series length. Series can be
+        different lengths.
 
         Parameters
         ----------
-        x: np.ndarray (2d array of shape (1,m1)).
+        x: np.ndarray (2d array of shape (d,m1)).
             First time series.
-        y: np.ndarray (2d array of shape (1,m2)).
+        y: np.ndarray (2d array of shape (d,m2)).
             Second time series.
-        c: float
-            parameter used in MSM (update later!)
+        c: float, default = 1.0
+            Cost for split or merge operation.
+        window: float, defaults = None
+            Float that is the radius of the sakoe chiba window (if using Sakoe-Chiba
+            lower bounding). Must be between 0 and 1.
+        itakura_max_slope: float, defaults = None
+            Gradient of the slope for itakura parallelogram (if using Itakura
+            Parallelogram lower bounding). Must be between 0 and 1.
+        bounding_matrix: np.ndarray (2d array of shape (m1,m2)), defaults = None
+            Custom bounding matrix to use. If defined then other lower_bounding params
+            are ignored. The matrix should be structure so that indexes considered in
+            bound should be the value 0. and indexes outside the bounding matrix should
+            be infinity.
+        kwargs: any
+            extra kwargs.
 
         Returns
         -------
@@ -176,69 +201,44 @@ class _MsmDistance(NumbaDistance):
             _x: np.ndarray,
             _y: np.ndarray,
         ) -> float:
-            cost_matrix = _cost_matrix(_x, _y, c, _bounding_matrix)
+            cost_matrix = _cost_matrix(_x, _y, _bounding_matrix, c)
             return cost_matrix[-1, -1]
 
         return numba_msm_distance
 
 
-@njit(cache=True, fastmath=True)
-def _dimension_sum(x: np.ndarray, j: int):
-    total = 0
-    for i in range(x.shape[0]):
-        total += x[i][j]
-
-    return total
-
-
 @njit(cache=True)
-def _calc_cost_cell(
-    new_point: np.ndarray,
-    x: np.ndarray,
-    y: np.ndarray,
-    c: float,
-) -> float:
-    """Cost calculation function for MSM."""
-    new_point_sum = _dimension_sum(new_point)
-    x_sum = _dimension_sum(x)
-    y_sum = _dimension_sum(y)
-
-    if ((x_sum <= new_point_sum) and (new_point_sum <= y_sum)) or (
-        (y_sum <= new_point_sum) and new_point_sum <= x_sum
-    ):
-        return c
-    else:
-        a = np.abs(new_point_sum - x_sum)
-        b = np.abs(new_point_sum - y_sum)
-
-        if a < b:
-            return c + a
-        else:
-            return c + b
-        # return c + np.min([np.abs(new_point - x), np.abs(new_point - y)])
+def _cost_function(x: float, y: float, z: float, c: float) -> float:
+    """Compute the cost function for the MSM algorithm.
 
+    Parameters
+    ----------
+    x: float
+        First value.
+    y: float
+        Second value.
+    z: float
+        Third value.
+    c: float
+        Parameter used in MSM (update later!)
 
-@njit(cache=True)
-def _cost_function(x: float, y: float, z: float, c: float) -> float:
-    if (y <= x and x <= z) or (y >= x and x >= z):
+    Returns
+    -------
+    float
+        Cost function value.
+    """
+    if (y <= x <= z) or (y >= x >= z):
         return c
-    a = x - y
-    if a < 0:
-        a = -a
-    b = x - z
-    if b < 0:
-        b = -b
-    if a > b:
-        return c + b
-    return c + a
+
+    return c + min(abs(x - y), abs(x - z))
 
 
 @njit(cache=True)
 def _cost_matrix(
     x: np.ndarray,
     y: np.ndarray,
-    c: float,
     bounding_matrix: np.ndarray,
+    c: float,
 ) -> float:
     """MSM distance compiled to no_python.
 
@@ -276,16 +276,9 @@ def _cost_matrix(
     for i in range(1, x_size):
         for j in range(1, y_size):
             if np.isfinite(bounding_matrix[i, j]):
-                d1 = cost[i - 1][j - 1] + np.abs(x[0][i] - y[0][j])
-                d2 = cost[i - 1][j] + _cost_function(x[0][i], x[0][i - 1], y[0][j], c)
-                d3 = cost[i][j - 1] + _cost_function(y[0][j], x[0][i], y[0][j - 1], c)
-
-            temp = d1
-            if d2 < temp:
-                temp = d2
-            if d3 < temp:
-                temp = d3
-
-            cost[i][j] = temp
+                d1 = cost[i - 1, j - 1] + np.abs(x[0][i] - y[0][j])
+                d2 = cost[i - 1, j] + _cost_function(x[0][i], x[0][i - 1], y[0][j], c)
+                d3 = cost[i, j - 1] + _cost_function(y[0][j], x[0][i], y[0][j - 1], c)
+                cost[i][j] = min(d1, d2, d3)
 
     return cost[0:, 0:]


Test patch:
diff --git a/sktime/classification/distance_based/tests/test_time_series_neighbors.py b/sktime/classification/distance_based/tests/test_time_series_neighbors.py
index 59f2d3404..513b5cb99 100644
--- a/sktime/classification/distance_based/tests/test_time_series_neighbors.py
+++ b/sktime/classification/distance_based/tests/test_time_series_neighbors.py
@@ -28,6 +28,17 @@ expected_correct = {
     "edr": 20,
 }
 
+# expected correct on test set using window params.
+expected_correct_window = {
+    "euclidean": 19,
+    "dtw": 21,
+    "wdtw": 21,
+    "msm": 10,
+    "erp": 19,
+    "edr": 20,
+    "lcss": 12,
+}
+
 
 @pytest.mark.parametrize("distance_key", distance_functions)
 def test_knn_on_unit_test(distance_key):
@@ -61,4 +72,4 @@ def test_knn_bounding_matrix(distance_key):
     for j in range(0, len(pred)):
         if pred[j] == y_test[j]:
             correct = correct + 1
-    assert correct == expected_correct[distance_key]
+    assert correct == expected_correct_window[distance_key]
diff --git a/sktime/distances/tests/_expected_results.py b/sktime/distances/tests/_expected_results.py
index 54161e854..3943ae416 100644
--- a/sktime/distances/tests/_expected_results.py
+++ b/sktime/distances/tests/_expected_results.py
@@ -13,7 +13,7 @@ _expected_distance_results = {
     "ddtw": [0.0, 2.0884818837222006, 34.837800040564005],
     "wdtw": [12.343758137512241, 0.985380547171357, 21.265839226825413],
     "wddtw": [0.0, 1.0442409418611003, 17.418900020282003],
-    "msm": [5.0, 3.0922226886554434, None],
+    "msm": [5.0, 6.828557434224288, None],
     "twe": [5.0, 11.548698748091073, 39.87793560457224],
 }