fix: indexing in kdtw (#2826)

Author: SebastianSchmidl

aeon/clustering/_kernel_k_means.py

@@ -26,20 +26,20 @@ def _kdtw_lk(x, y, local_kernel):
     diagonal_weights = np.zeros(max(x_timepoints, y_timepoints))
 
     min_timepoints = min(x_timepoints, y_timepoints)
-    diagonal_weights[1] = 1.0
+    diagonal_weights[0] = 1.0
     for i in range(1, min_timepoints):
         diagonal_weights[i] = local_kernel[i - 1, i - 1]
 
     cost_matrix[0, 0] = 1
     cumulative_dp_diag[0, 0] = 1
 
     for i in range(1, x_timepoints):
-        cost_matrix[i, 1] = cost_matrix[i - 1, 1] * local_kernel[i - 1, 2]
-        cumulative_dp_diag[i, 1] = cumulative_dp_diag[i - 1, 1] * diagonal_weights[i]
+        cost_matrix[i, 0] = cost_matrix[i - 1, 0] * local_kernel[i - 1, 0]
+        cumulative_dp_diag[i, 0] = cumulative_dp_diag[i - 1, 0] * diagonal_weights[i]
 
     for j in range(1, y_timepoints):
-        cost_matrix[1, j] = cost_matrix[1, j - 1] * local_kernel[2, j - 1]
-        cumulative_dp_diag[1, j] = cumulative_dp_diag[1, j - 1] * diagonal_weights[j]
+        cost_matrix[0, j] = cost_matrix[0, j - 1] * local_kernel[0, j - 1]
+        cumulative_dp_diag[0, j] = cumulative_dp_diag[0, j - 1] * diagonal_weights[j]
 
     for i in range(1, x_timepoints):
         for j in range(1, y_timepoints):

aeon/clustering/tests/test_kernel_k_means.py

@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from aeon.clustering._kernel_k_means import TimeSeriesKernelKMeans
+from aeon.clustering._kernel_k_means import TimeSeriesKernelKMeans, _kdtw
 from aeon.datasets import load_basic_motions
 from aeon.utils.validation._dependencies import _check_estimator_deps
 
@@ -19,18 +19,18 @@
 
 expected_results_kdtw = [0, 2, 0, 0, 0]
 
+max_train = 5
+
+X_train, y_train = load_basic_motions(split="train")
+X_test, y_test = load_basic_motions(split="test")
+
 
 @pytest.mark.skipif(
     not _check_estimator_deps(TimeSeriesKernelKMeans, severity="none"),
     reason="skip test if required soft dependencies not available",
 )
-def test_kernel_k_means():
-    """Test implementation of kernel k means."""
-    max_train = 5
-
-    X_train, y_train = load_basic_motions(split="train")
-    X_test, y_test = load_basic_motions(split="test")
-
+def test_kernel_k_means_gak():
+    """Test implementation of kernel k means with GAK kernel."""
     kernel_kmeans = TimeSeriesKernelKMeans(random_state=1, n_clusters=3)
     kernel_kmeans.fit(X_train[0:max_train])
     test_shape_result = kernel_kmeans.predict(X_test[0:max_train])
@@ -44,6 +44,13 @@ def test_kernel_k_means():
     for val in proba:
         assert np.count_nonzero(val == 1.0) == 1
 
+
+@pytest.mark.skipif(
+    not _check_estimator_deps(TimeSeriesKernelKMeans, severity="none"),
+    reason="skip test if required soft dependencies not available",
+)
+def test_kernel_k_means_kdtw():
+    """Test implementation of kernel k means with KDTW kernel."""
     kernel_kmeans_kdtw = TimeSeriesKernelKMeans(
         kernel="kdtw",
         random_state=1,
@@ -61,3 +68,17 @@ def test_kernel_k_means():
 
     for val in kdtw_proba:
         assert np.count_nonzero(val == 1.0) == 1
+
+
+def test_kdtw_kernel_univariate():
+    """Test kdtw kernel for univariate time series."""
+    # expected value created with the original (Matlab) code from:
+    # https://people.irisa.fr/Pierre-Francois.Marteau/REDK/KDTW/KDTW.html
+    x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9], dtype=np.float64).reshape(-1, 1)
+    y = np.array([5, 6, 7, 8, 9, 1, 2], dtype=np.float64).reshape(-1, 1)
+    sigma = 0.125
+    epsilon = 1e-20
+    expected_distance = 1.2814e-102
+
+    distance = _kdtw(x, y, sigma=sigma, epsilon=epsilon)
+    np.testing.assert_allclose(expected_distance, distance, rtol=1e-4, atol=1e-106)