Tests and peak finding examples

examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt

@@ -5,20 +5,24 @@
 
 package space.kscience.kmath.series
 
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
 import space.kscience.kmath.structures.*
 import space.kscience.kmath.operations.invoke
 import space.kscience.plotly.*
 import space.kscience.plotly.models.Scatter
 import kotlin.math.sin
 
-fun main(): Unit = (Double.seriesAlgebra()) {
+private val customAlgebra = (Double.algebra.bufferAlgebra) { SeriesAlgebra(this) { it.toDouble() } }
+
+fun main(): Unit = (customAlgebra) {
     val signal = DoubleArray(800) {
         sin(it.toDouble() / 10.0) + 3.5 * sin(it.toDouble() / 60.0)
     }.asBuffer().moveTo(0)
-
     val emd = empiricalModeDecomposition(
         sConditionThreshold = 1,
         maxSiftIterations = 15,
+        siftingDelta = 1e-2,
         nModes = 4
     ).decompose(signal)
     println("EMD: ${emd.modes.size} modes extracted, terminated because ${emd.terminatedBecause}")
@@ -26,7 +30,7 @@ fun main(): Unit = (Double.seriesAlgebra()) {
     fun Plot.series(name: String, buffer: Buffer<Double>, block: Scatter.() -> Unit = {}) {
         this.scatter {
             this.name = name
-            this.x.numbers = buffer.offsetIndices
+            this.x.numbers = buffer.labels
             this.y.doubles = buffer.toDoubleArray()
             block()
         }

examples/src/main/kotlin/space/kscience/kmath/series/peakFinding.kt

@@ -0,0 +1,85 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.series
+
+import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.*
+import space.kscience.plotly.*
+import space.kscience.plotly.models.Scatter
+import space.kscience.plotly.models.ScatterMode
+import kotlin.math.sin
+
+private val customAlgebra = (Double.algebra.bufferAlgebra) { SeriesAlgebra(this) { it * 50.0 / 599.0 } }
+
+fun main(): Unit = (customAlgebra) {
+    /*
+    val signal = DoubleArray(600) {
+        val x = it * 50.0 / 599
+        (3.0 * sin(x) + 0.5 * cos(7.0 * x)).coerceIn(-3.0 .. 3.0)
+    }.asBuffer().moveTo(0)
+    val peaks = signal.peaks()
+    val troughs = signal.troughs()
+    println(peaks)
+    println(troughs)
+
+    fun Plot.series(name: String, buffer: Buffer<Double>, block: Scatter.() -> Unit = {}) {
+        scatter {
+            this.name = name
+            this.x.numbers = buffer.labels
+            this.y.doubles = buffer.toDoubleArray()
+            block()
+        }
+    }
+    Plotly.plot {
+        series("Signal", signal)
+        scatter {
+            name = "Peaks"
+            mode = ScatterMode.markers
+            x.doubles = peaks.map { signal.labels[it] }.toDoubleArray()
+            y.doubles = peaks.map { signal[it] }.toDoubleArray()
+        }
+        scatter {
+            name = "Troughs"
+            mode = ScatterMode.markers
+            x.doubles = troughs.map {  signal.labels[it] }.toDoubleArray()
+            y.doubles = troughs.map { signal[it] }.toDoubleArray()
+        }
+    }.makeFile(resourceLocation = ResourceLocation.REMOTE)
+     */
+    val nSamples = 600
+    val signal = DoubleArray(nSamples) {
+        val x = it * 12.0 / (nSamples - 1)
+        (3.5 * sin(x)).coerceIn(-3.0 .. 3.0)
+    }.asBuffer().moveTo(0)
+    val peaks = signal.peaks(PlateauEdgePolicy.KEEP_ALL_EDGES)
+    val troughs = signal.troughs(PlateauEdgePolicy.KEEP_ALL_EDGES)
+    println(peaks)
+    println(troughs)
+
+    fun Plot.series(name: String, buffer: Buffer<Double>, block: Scatter.() -> Unit = {}) {
+        scatter {
+            this.name = name
+            this.x.numbers = buffer.labels
+            this.y.doubles = buffer.toDoubleArray()
+            block()
+        }
+    }
+    Plotly.plot {
+        series("Signal", signal)
+        scatter {
+            name = "Peaks"
+            mode = ScatterMode.markers
+            x.doubles = peaks.map { signal.labels[it] }.toDoubleArray()
+            y.doubles = peaks.map { signal[it] }.toDoubleArray()
+        }
+        scatter {
+            name = "Troughs"
+            mode = ScatterMode.markers
+            x.doubles = troughs.map {  signal.labels[it] }.toDoubleArray()
+            y.doubles = troughs.map { signal[it] }.toDoubleArray()
+        }
+    }.makeFile(resourceLocation = ResourceLocation.REMOTE)
+}

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

@@ -8,11 +8,8 @@ package space.kscience.kmath.series
 import space.kscience.kmath.interpolation.SplineInterpolator
 import space.kscience.kmath.interpolation.interpolate
 import space.kscience.kmath.operations.*
-import space.kscience.kmath.operations.Float64BufferOps.Companion.div
-import space.kscience.kmath.operations.Float64BufferOps.Companion.pow
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.asBuffer
-import kotlin.math.sign
+import space.kscience.kmath.structures.last
 
 /**
  * Empirical mode decomposition of a signal represented as a [Series].
@@ -29,13 +26,13 @@ import kotlin.math.sign
  * @param nModes how many modes should be extracted at most. The algorithm may return fewer modes if it was not
  * possible to extract more modes from the signal.
  */
-public class EmpiricalModeDecomposition<BA, L: Number> (
-    private val seriesAlgebra: SeriesAlgebra<Double, *, BA, L>,
+public class EmpiricalModeDecomposition<T: Comparable<T>, A: Field<T>, BA, L: T> (
+    private val seriesAlgebra: SeriesAlgebra<T, A, BA, L>,
     private val sConditionThreshold: Int = 15,
     private val maxSiftIterations: Int = 20,
-    private val siftingDelta: Double = 1e-2,
+    private val siftingDelta: T,
     private val nModes: Int = 6
-) where BA: BufferAlgebra<Double, *>,  BA: RingOps<Buffer<Double>> {
+) where BA: BufferAlgebra<T, A>,  BA: FieldOps<Buffer<T>> {
 
     /**
      * Take a signal, construct an upper and a lower envelopes, find the mean value of two,
@@ -45,37 +42,38 @@ public class EmpiricalModeDecomposition<BA, L: Number> (
      * @return mean [Series] or `null`. `null` is returned in case
      * the signal does not have enough extrema to construct envelopes.
      */
-    private fun findMean(signal: Series<Double>): Series<Double>? = (seriesAlgebra) {
-        val interpolator = SplineInterpolator(Float64Field)
-        fun generateEnvelope(extrema: List<Int>, paddedExtremeValues: DoubleArray): Series<Double> {
+    private fun findMean(signal: Series<T>): Series<T>? = (seriesAlgebra) {
+        val interpolator = SplineInterpolator(elementAlgebra)
+        val makeBuffer = elementAlgebra.bufferFactory
+        fun generateEnvelope(extrema: List<Int>, paddedExtremeValues: Buffer<T>): Series<T> {
             val envelopeFunction = interpolator.interpolate(
-                Buffer(extrema.size) { signal.labels[extrema[it]].toDouble() },
-                paddedExtremeValues.asBuffer()
+                makeBuffer(extrema.size) { signal.labels[extrema[it]] },
+                paddedExtremeValues
             )
             return signal.mapWithLabel { _, label ->
                 // For some reason PolynomialInterpolator is exclusive and the right boundary
                 // TODO Notify interpolator authors
-                envelopeFunction(label.toDouble()) ?: paddedExtremeValues.last()
+                envelopeFunction(label) ?: paddedExtremeValues.last()
                 // need to make the interpolator yield values outside boundaries?
             }
         }
         // Extrema padding (experimental) TODO padding needs a dedicated function
         val maxima = listOf(0) + signal.peaks() + (signal.size - 1)
-        val maxValues = DoubleArray(maxima.size) { signal[maxima[it]] }
+        val maxValues = makeBuffer(maxima.size) { signal[maxima[it]] }
         if (maxValues[0] < maxValues[1]) {
             maxValues[0] = maxValues[1]
         }
-        if (maxValues.last() < maxValues[maxValues.lastIndex - 1]) {
-            maxValues[maxValues.lastIndex] = maxValues[maxValues.lastIndex - 1]
+        if (maxValues.last() < maxValues[maxValues.size - 2]) {
+            maxValues[maxValues.size - 1] = maxValues[maxValues.size - 2]
         }
 
         val minima = listOf(0) + signal.troughs() + (signal.size - 1)
-        val minValues = DoubleArray(minima.size) { signal[minima[it]] }
+        val minValues = makeBuffer(minima.size) { signal[minima[it]] }
         if (minValues[0] > minValues[1]) {
             minValues[0] = minValues[1]
         }
-        if (minValues.last() > minValues[minValues.lastIndex - 1]) {
-            minValues[minValues.lastIndex] = minValues[minValues.lastIndex - 1]
+        if (minValues.last() > minValues[minValues.size - 2]) {
+            minValues[minValues.size - 1] = minValues[minValues.size - 2]
         }
          return if (maxima.size < 3 || minima.size < 3) null else { // maybe make an early return?
             val upperEnvelope = generateEnvelope(maxima, maxValues)
@@ -92,13 +90,13 @@ public class EmpiricalModeDecomposition<BA, L: Number> (
      * @return [SiftingResult.NotEnoughExtrema] is returned if the signal has too few extrema to extract a mode.
      * Success of an appropriate type (See [SiftingResult.Success] class) is returned otherwise.
      */
-    private fun sift(signal: Series<Double>): SiftingResult = siftInner(signal, 1, 0)
+    private fun sift(signal: Series<T>): SiftingResult = siftInner(signal, 1, 0)
 
     /**
      * Compute a single iteration of the sifting process.
      */
     private tailrec fun siftInner(
-        prevMode: Series<Double>,
+        prevMode: Series<T>,
         iterationNumber: Int,
         sNumber: Int
     ): SiftingResult = (seriesAlgebra) {
@@ -106,11 +104,12 @@ public class EmpiricalModeDecomposition<BA, L: Number> (
             return if (iterationNumber == 1) SiftingResult.NotEnoughExtrema
             else SiftingResult.SignalFlattened(prevMode)
         val mode = prevMode.zip(mean) { p, m -> p - m }
-        val newSNumber = if (mode.sCondition()) sNumber + 1 else sNumber
+        val newSNumber = if (sCondition(mode)) sNumber + 1 else sNumber
         return when {
             iterationNumber >= maxSiftIterations -> SiftingResult.MaxIterationsReached(mode)
             sNumber >= sConditionThreshold -> SiftingResult.SNumberReached(mode)
-            relativeDifference(mode, prevMode) < siftingDelta * mode.size -> SiftingResult.DeltaReached(mode)
+            relativeDifference(mode, prevMode) < (elementAlgebra) { siftingDelta * mode.size } ->
+                SiftingResult.DeltaReached(mode)
             else -> siftInner(mode, iterationNumber + 1, newSNumber)
         }
     }
@@ -123,8 +122,8 @@ public class EmpiricalModeDecomposition<BA, L: Number> (
      * Modes returned in a list which contains as many modes as it was possible
      * to extract before triggering one of the termination conditions.
      */
-     public fun decompose(signal: Series<Double>): EMDecompositionResult = (seriesAlgebra) {
-        val modes = mutableListOf<Series<Double>>()
+     public fun decompose(signal: Series<T>): EMDecompositionResult<T> = (seriesAlgebra) {
+        val modes = mutableListOf<Series<T>>()
         var residual = signal
         repeat(nModes) {
             val nextMode = when(val r = sift(residual)) {
@@ -132,14 +131,15 @@ public class EmpiricalModeDecomposition<BA, L: Number> (
                     return EMDecompositionResult(
                         if (it == 0) EMDTerminationReason.SIGNAL_TOO_FLAT
                         else EMDTerminationReason.ALL_POSSIBLE_MODES_EXTRACTED,
-                        modes
+                        modes,
+                        residual
                     )
-                is SiftingResult.Success -> r.result
+                is SiftingResult.Success<*> -> r.result
             }
-            modes.add(nextMode)
+            modes.add(nextMode as Series<T>) // TODO remove unchecked cast
             residual = residual.zip(nextMode) { l, r -> l - r }
         }
-        return EMDecompositionResult(EMDTerminationReason.MAX_MODES_REACHED, modes)
+        return EMDecompositionResult(EMDTerminationReason.MAX_MODES_REACHED, modes, residual)
     }
 }
 
@@ -157,13 +157,13 @@ public class EmpiricalModeDecomposition<BA, L: Number> (
  * @param nModes how many modes should be extracted at most. The algorithm may return fewer modes if it was not
  * possible to extract more modes from the signal.
  */
-public fun <L: Number, BA> SeriesAlgebra<Double, *, BA, L>.empiricalModeDecomposition(
+public fun <T: Comparable<T>, L: T, A: Field<T>, BA> SeriesAlgebra<T, A, BA, L>.empiricalModeDecomposition(
     sConditionThreshold: Int = 15,
     maxSiftIterations: Int = 20,
-    siftingDelta: Double = 1e-2,
+    siftingDelta: T,
     nModes: Int = 3
-): EmpiricalModeDecomposition<BA, L>
-where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>> = EmpiricalModeDecomposition(
+): EmpiricalModeDecomposition<T, A, BA, L>
+where BA: BufferAlgebra<T, A>, BA: FieldOps<Buffer<T>> = EmpiricalModeDecomposition(
     seriesAlgebra = this,
     sConditionThreshold = sConditionThreshold,
     maxSiftIterations = maxSiftIterations,
@@ -174,12 +174,15 @@ where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>> = EmpiricalModeD
 /**
  * Brute force count all zeros in the series.
  */
-private fun Series<Double>.countZeros(): Int {
-    require(size >= 2) { "Expected series with at least 2 elements, but got $size elements" }
-    data class SignCounter(val prevSign: Double, val zeroCount: Int)
+internal fun <T: Comparable<T>, A: Ring<T>, BA> SeriesAlgebra<T, A, BA, *>.countZeros(
+    signal: Series<T>
+): Int where BA: BufferAlgebra<T, A>, BA: FieldOps<Buffer<T>> {
+    require(signal.size >= 2) { "Expected series with at least 2 elements, but got ${signal.size} elements" }
+    data class SignCounter(val prevSign: Int, val zeroCount: Int)
+    fun strictSign(arg: T): Int = if (arg > elementAlgebra.zero) 1 else -1
 
-    return fold(SignCounter(sign(get(0)), 0)) { acc: SignCounter, it: Double ->
-        val currentSign = sign(it)
+    return signal.fold(SignCounter(strictSign(signal[0]), 0)) { acc, it ->
+        val currentSign = strictSign(it)
         if (acc.prevSign != currentSign) SignCounter(currentSign, acc.zeroCount + 1)
         else SignCounter(currentSign, acc.zeroCount)
     }.zeroCount
@@ -188,18 +191,19 @@ private fun Series<Double>.countZeros(): Int {
 /**
  * Compute relative difference of two series.
  */
-private fun <BA> SeriesAlgebra<Double, *, BA, *>.relativeDifference(
-    current: Series<Double>,
-    previous: Series<Double>
-):Double where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>> =
-    (current - previous).pow(2)
-        .div(previous pow 2)
-        .fold(0.0) { acc, d -> acc + d } // TODO replace with Series<>.sum() method when it's implemented
+private fun <T, A: Ring<T>, BA> SeriesAlgebra<T, A, BA, *>.relativeDifference(
+    current: Series<T>,
+    previous: Series<T>
+): T where BA: BufferAlgebra<T, A>, BA: FieldOps<Buffer<T>> = (bufferAlgebra) {
+    ((current - previous) * (current - previous))
+        .div(previous * previous)
+        .fold(elementAlgebra.zero) { acc, it -> acc + it}
+}
 
 /**
  * Brute force count all extrema of a series.
  */
-private fun Series<Double>.countExtrema(): Int {
+internal fun <T: Comparable<T>> Series<T>.countExtrema(): Int {
     require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" }
     return peaks().size + troughs().size
 }
@@ -208,7 +212,10 @@ private fun Series<Double>.countExtrema(): Int {
  * Check whether the numbers of zeroes and extrema of a series differ by no more than 1.
  * This is a necessary condition of an empirical mode.
  */
-private fun Series<Double>.sCondition(): Boolean = (countExtrema() - countZeros()) in -1..1
+private fun <T: Comparable<T>, A: Ring<T>, BA> SeriesAlgebra<T, A, BA, *>.sCondition(
+    signal: Series<T>
+): Boolean where BA: BufferAlgebra<T, A>, BA: FieldOps<Buffer<T>> =
+    (signal.countExtrema() - countZeros(signal)) in -1..1
 
 internal sealed interface SiftingResult {
 
@@ -216,33 +223,33 @@ internal sealed interface SiftingResult {
      * Represents a condition when a mode has been successfully
      * extracted in a sifting process.
      */
-    open class Success(val result: Series<Double>): SiftingResult
+    open class Success<T>(val result: Series<T>): SiftingResult
 
     /**
      * Returned when no termination condition was reached and the proto-mode
      * has become too flat (with not enough extrema to build envelopes)
      * after several sifting iterations.
      */
-    class SignalFlattened(result: Series<Double>) : Success(result)
+    class SignalFlattened<T>(result: Series<T>) : Success<T>(result)
 
     /**
      * Returned when sifting process has been terminated due to the
      * S-number condition being reached.
      */
-    class SNumberReached(result: Series<Double>) : Success(result)
+    class SNumberReached<T>(result: Series<T>) : Success<T>(result)
 
     /**
      * Returned when sifting process has been terminated due to the
      * delta condition (Cauchy criterion) being reached.
      */
-    class DeltaReached(result: Series<Double>) : Success(result)
+    class DeltaReached<T>(result: Series<T>) : Success<T>(result)
 
     /**
      * Returned when sifting process has been terminated after
      * executing the maximum number of iterations (specified when creating an instance
      * of [EmpiricalModeDecomposition]).
      */
-    class MaxIterationsReached(result: Series<Double>): Success(result)
+    class MaxIterationsReached<T>(result: Series<T>): Success<T>(result)
 
     /**
      * Returned when the submitted signal has not enough extrema to build envelopes,
@@ -274,7 +281,8 @@ public enum class EMDTerminationReason {
     ALL_POSSIBLE_MODES_EXTRACTED
 }
 
-public data class EMDecompositionResult(
+public data class EMDecompositionResult<T>(
     val terminatedBecause: EMDTerminationReason,
-    val modes: List<Series<Double>>
+    val modes: List<Series<T>>,
+    val residual: Series<T>
 )

kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestEmd.kt

@@ -0,0 +1,66 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.series
+
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.asBuffer
+import kotlin.math.sin
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertTrue
+import kotlin.random.Random
+
+
+class TestEmd {
+    companion object{
+        val testAlgebra = (Double.algebra.bufferAlgebra) { SeriesAlgebra(this) { it.toDouble() } }
+    }
+
+    @Test
+    fun testBasic() = (testAlgebra) {
+        val signal = DoubleArray(800) {
+            sin(it.toDouble() / 10.0) + 3.5 * sin(it.toDouble() / 60.0)
+        }.asBuffer().moveTo(0)
+        val emd = empiricalModeDecomposition(
+            sConditionThreshold = 1,
+            maxSiftIterations = 15,
+            siftingDelta = 1e-2,
+            nModes = 4
+        ).decompose(signal)
+
+        assertEquals(emd.modes.size, 3)
+        emd.modes.forEach { imf ->
+            assertTrue(imf.peaks().size - imf.troughs().size in -1..1)
+        }
+    }
+
+    @Test
+    fun testNoiseFiltering() = (testAlgebra) {
+        val signal = DoubleArray(800) {
+            sin(it.toDouble() / 30.0) + 2.0 * (Random.nextDouble() - 0.5)
+        }.asBuffer().moveTo(0)
+        val emd = empiricalModeDecomposition(
+            sConditionThreshold = 10,
+            maxSiftIterations = 15,
+            siftingDelta = 1e-2,
+            nModes = 10
+        ).decompose(signal)
+        // Check whether the signal with the expected frequency is present
+        assertEquals(emd.modes.count { it.countExtrema() in 7..9 }, 1)
+    }
+
+    @Test
+    fun testZeros() = (testAlgebra) {
+        val nSamples = 200
+        // sin(10*x) where x in [0, 1)
+        val signal = DoubleArray(nSamples) {
+            sin(it * 10.0 / (nSamples - 1))
+        }.asBuffer().moveTo(0)
+        assertEquals(countZeros(signal), 4)
+    }
+}

kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestPeakFinding.kt

@@ -0,0 +1,39 @@
+/*
+ * Copyright 2018-2024 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
+ */
+
+package space.kscience.kmath.series
+
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.asBuffer
+import kotlin.math.sin
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class TestPeakFinding {
+    companion object {
+        val testAlgebra = (Double.algebra.bufferAlgebra) { SeriesAlgebra(this) { it.toDouble() } }
+    }
+
+    @Test
+    fun testPeakFinding() = (testAlgebra) {
+        val nSamples = 600
+        val signal = DoubleArray(nSamples) {
+            val x = it * 12.0 / (nSamples - 1)
+            (3.5 * sin(x)).coerceIn(-3.0 .. 3.0)
+        }.asBuffer().moveTo(0)
+
+        val peaksAvg = signal.peaks(PlateauEdgePolicy.AVERAGE)
+        val troughsAvg = signal.troughs(PlateauEdgePolicy.AVERAGE)
+        assertEquals(peaksAvg.size, 2)
+        assertEquals(troughsAvg.size, 2)
+
+        val peaksBoth = signal.peaks(PlateauEdgePolicy.KEEP_ALL_EDGES)
+        val troughsBoth = signal.peaks(PlateauEdgePolicy.KEEP_ALL_EDGES)
+        assertEquals(peaksBoth.size, 4)
+        assertEquals(troughsBoth.size, 4)
+    }
+}