fix(kernelcpd): bug fix in pelt

ruptures/detection/_detection/ekcpd_pelt_computation.c

@@ -3,6 +3,12 @@
 #include <math.h>
 
 #include "kernels.h"
+static inline int max_int(int a, int b)
+{
+    if (a > b)
+        return a;
+    return b;
+}
 
 /**
  * @brief Efficient kernel change point detection
@@ -37,13 +43,33 @@ void ekcpd_pelt_compute(double *signal, int n_samples, int n_dims, double beta,
         D[t] = 0.0;
         S[t] = 0.0;
         M_V[t] = 0.0;
-        M_path[t] = 0.0;
+        M_path[t] = 0;
         M_pruning[t] = 0.0;
     }
 
+    // for t<2*min_size, there cannot be any change point.
+    for (t = 1; t < 2 * min_size; t++)
+    {
+        diag_element = kernel_value_by_name(&(signal[(t - 1) * n_dims]), &(signal[(t - 1) * n_dims]), n_dims, kernelDescObj);
+        D[t] = D[t - 1] + diag_element;
+
+        // Compute S[t-1] = S_{t-1, t}, S[t-2] = S_{t-2, t}, ..., S[0] = S_{0, t}
+        // S_{t-1, t} can be computed with S_{t-1, t-1}.
+        // S_{t-1, t-1} was stored in S[t-1]
+        // S_{t-1, t} will be stored in S[t-1] as well
+        c_r = 0.0;
+        for (s = t - 1; s >= 0; s--)
+        {
+            c_r += kernel_value_by_name(&(signal[s * n_dims]), &(signal[(t - 1) * n_dims]), n_dims, kernelDescObj);
+            S[s] += 2 * c_r - diag_element;
+        }
+        c_cost = D[t] - D[s] - S[0] / t;
+        M_V[t] = c_cost + beta;
+    }
+
     // Computation loop
     // Handle y_{0..t} = {y_0, ..., y_{t-1}}
-    for (t = 1; t < (n_samples + 1); t++)
+    for (t = 2 * min_size; t < (n_samples + 1); t++)
     {
         diag_element = kernel_value_by_name(&(signal[(t - 1) * n_dims]), &(signal[(t - 1) * n_dims]), n_dims, kernelDescObj);
         D[t] = D[t - 1] + diag_element;
@@ -70,7 +96,7 @@ void ekcpd_pelt_compute(double *signal, int n_samples, int n_dims, double beta,
         M_V[t] = c_cost_sum;
         M_path[t] = s;
         // search for minimum (penalized) sum of cost
-        for (s = min_size + 1; s < t - min_size + 1; s++)
+        for (s = max_int(s_min, min_size) + 1; s < t - min_size + 1; s++)
         {
             // Compute cost on y_{s..t}
             // D_{s..t} = D_{0..t} - D{0..s} <--> D_{s..t} = D[t] - D[s]
@@ -89,7 +115,14 @@ void ekcpd_pelt_compute(double *signal, int n_samples, int n_dims, double beta,
         // Pruning
         while ((M_pruning[s_min] >= M_V[t]) && (s_min < t - min_size + 1))
         {
-            s_min++;
+            if (s_min == 0)
+            {
+                s_min += min_size;
+            }
+            else
+            {
+                s_min++;
+            }
         }
     }