kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/DynamicTimeWarping.kt

@@ -5,70 +5,35 @@
 
 package space.kscience.kmath.series
 
+import space.kscience.kmath.PerformancePitfall
 import space.kscience.kmath.nd.*
 import space.kscience.kmath.nd.DoubleBufferND
 import space.kscience.kmath.nd.ShapeND
-import space.kscience.kmath.nd.ndAlgebra
-import space.kscience.kmath.operations.DoubleField
-import space.kscience.kmath.structures.DoubleBuffer
 import kotlin.math.abs
 
-public const val LEFT_OFFSET : Int = -1
-public const val BOTTOM_OFFSET : Int = 1
-public const val DIAGONAL_OFFSET : Int = 0
-
-
-
-// TODO: Change container for alignMatrix to kmath special ND structure
-public data class DynamicTimeWarpingData(val totalCost : Double = 0.0,
-                                     val alignMatrix : IntBufferND = IntRingOpsND.structureND(ShapeND(0, 0)) { (i, j) -> 0}) {
-    override fun equals(other: Any?): Boolean {
-        if (this === other) return true
-        if (other == null || this::class != other::class) return false
-
-        other as DynamicTimeWarpingData
-
-        if (totalCost != other.totalCost) return false
-
-        return true
-    }
-}
 
 /**
- * costMatrix calculates special matrix of costs alignment for two series.
+ *  Offset constants which will be used later. Added them for avoiding "magical numbers" problem.
- * Formula: costMatrix[i, j] = euqlidNorm(series1(i), series2(j)) + min(costMatrix[i - 1, j],
- *                                                                      costMatrix[i, j - 1],
- *                                                                      costMatrix[i - 1, j - 1]).
- * There is special cases for i = 0 or j = 0.
  */
+internal const val LEFT_OFFSET : Int = -1
+internal const val BOTTOM_OFFSET : Int = 1
+internal const val DIAGONAL_OFFSET : Int = 0
 
-public fun costMatrix(series1 : DoubleBuffer, series2 : DoubleBuffer) : DoubleBufferND {
+
-    val dtwMatrix = DoubleField.ndAlgebra.structureND(ShapeND(series1.size, series2.size)) {
+/**
-            (row, col) -> abs(series1[row] - series2[col])
+ * Public class to store result of method. Class contains total penalty cost for series alignment.
-    }
+ * Also, this class contains align matrix (which point of the first series matches to point of the other series).
-    for ( (row, col) in dtwMatrix.indices) {
+ */
-        dtwMatrix[row, col] += when {
+public data class DynamicTimeWarpingData(val totalCost : Double = 0.0,
-            row == 0 && col == 0 -> 0.0
+                                         val alignMatrix : DoubleBufferND = DoubleFieldOpsND.structureND(ShapeND(0, 0)) { (_, _) -> 0.0})
-            row == 0 -> dtwMatrix[row, col - 1]
-            col == 0 -> dtwMatrix[row - 1, col]
-            else -> minOf(
-                dtwMatrix[row, col - 1],
-                dtwMatrix[row - 1, col],
-                dtwMatrix[row - 1, col - 1]
-            )
-        }
-    }
-    return dtwMatrix
-}
 
 /**
  * PathIndices class for better code perceptibility.
  * Special fun moveOption represent offset for indices. Arguments of this function
  * is flags for bottom, diagonal or left offsets respectively.
  */
-
+internal data class PathIndices (var id_x: Int, var id_y: Int) {
-public data class PathIndices (var id_x: Int, var id_y: Int) {
+    fun moveOption (direction: Int) {
-    public fun moveOption (direction: Int) {
         when(direction) {
             BOTTOM_OFFSET -> id_x--
             DIAGONAL_OFFSET -> {
@@ -85,16 +50,36 @@ public data class PathIndices (var id_x: Int, var id_y: Int) {
  * Final DTW method realization. Returns alignment matrix
  * for two series comparing and penalty for this alignment.
  */
-
+@OptIn(PerformancePitfall::class)
-public fun dynamicTimeWarping(series1 : DoubleBuffer, series2 : DoubleBuffer) : DynamicTimeWarpingData {
+public fun DoubleFieldOpsND.dynamicTimeWarping(series1 : Series<Double>, series2 : Series<Double>) : DynamicTimeWarpingData {
     var cost = 0.0
     var pathLength = 0
-    val costMatrix = costMatrix(series1, series2)
+    // Special matrix of costs alignment for two series.
-    val alignMatrix: IntBufferND = IntRingOpsND.structureND(ShapeND(series1.size, series2.size)) {(row, col) -> 0}
+    val costMatrix = structureND(ShapeND(series1.size, series2.size)) {
+            (row, col) -> abs(series1[row] - series2[col])
+    }
+    // Formula: costMatrix[i, j] = euqlidNorm(series1(i), series2(j)) +
+    // min(costMatrix[i - 1, j], costMatrix[i, j - 1], costMatrix[i - 1, j - 1]).
+    for ( (row, col) in costMatrix.indices) {
+        costMatrix[row, col] += when {
+            // There is special cases for i = 0 or j = 0.
+            row == 0 && col == 0 -> 0.0
+            row == 0 -> costMatrix[row, col - 1]
+            col == 0 -> costMatrix[row - 1, col]
+            else -> minOf(
+                costMatrix[row, col - 1],
+                costMatrix[row - 1, col],
+                costMatrix[row - 1, col - 1]
+            )
+        }
+    }
+    // alignMatrix contains non-zero values at position where two points from series matches
+    // Values are penalty for concatenation of current points.
+    val alignMatrix = structureND(ShapeND(series1.size, series2.size)) {(_, _) -> 0.0}
     val indexes = PathIndices(alignMatrix.indices.shape.first() - 1,  alignMatrix.indices.shape.last() - 1)
 
     with(indexes) {
-        alignMatrix[id_x, id_y] = 1
+        alignMatrix[id_x, id_y] = costMatrix[id_x, id_y]
         cost += costMatrix[id_x, id_y]
         pathLength++
         while (id_x != 0 || id_y != 0) {
@@ -109,7 +94,7 @@ public fun dynamicTimeWarping(series1 : DoubleBuffer, series2 : DoubleBuffer) :
                     moveOption(DIAGONAL_OFFSET)
                 }
             }
-            alignMatrix[id_x, id_y] = 1
+            alignMatrix[id_x, id_y] = costMatrix[id_x, id_y]
             cost += costMatrix[id_x, id_y]
             pathLength++
         }