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

@@ -6,11 +6,12 @@
 package space.kscience.kmath.series
 
 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 = with(Double.seriesAlgebra()) {
+fun main(): Unit = (Double.seriesAlgebra()) {
     val signal = DoubleArray(800) {
         sin(it.toDouble() / 10.0) + 3.5 * sin(it.toDouble() / 60.0)
     }.asBuffer().moveTo(0)

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

@@ -12,7 +12,6 @@ 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.last
-import space.kscience.kmath.structures.slice
 import kotlin.math.sign
 
 /**
@@ -30,22 +29,23 @@ 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> (
+public class EmpiricalModeDecomposition<BA, L: Number> (
     private val seriesAlgebra: SeriesAlgebra<Double, *, BA, L>,
     private val sConditionThreshold: Int = 15,
     private val maxSiftIterations: Int = 20,
     private val siftingDelta: Double = 1e-2,
     private val nModes: Int = 6
-) where BA: BufferAlgebra<Double, *>,  BA: RingOps<Buffer<Double>>,  L : Number {
+) where BA: BufferAlgebra<Double, *>,  BA: RingOps<Buffer<Double>> {
 
     /**
-     * Take a signal, construct an upper and a lower envelope, find the mean value of two,
+     * Take a signal, construct an upper and a lower envelopes, find the mean value of two,
      * represented as a series.
      * @param signal Signal to compute on
      *
-     * @return null is returned if the signal has not enough extrema to construct envelopes.
+     * @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>? = with(seriesAlgebra) {
+    private fun findMean(signal: Series<Double>): Series<Double>? = (seriesAlgebra) {
         val interpolator = SplineInterpolator(Float64Field)
         fun generateEnvelope(extrema: List<Int>): Series<Double> {
             val envelopeFunction = interpolator.interpolate(
@@ -59,12 +59,12 @@ public class EmpiricalModeDecomposition<BA, L> (
                 // need to make the interpolator yield values outside boundaries?
             }
         }
-        val maxima = signal.maxima()
-        val minima = signal.minima()
+        val maxima = signal.paddedMaxima()
+        val minima = signal.paddedMinima()
         return if (maxima.size < 3 || minima.size < 3) null else {
             val upperEnvelope = generateEnvelope(maxima)
             val lowerEnvelope = generateEnvelope(minima)
-            return (upperEnvelope + lowerEnvelope).map { it * 0.5 }
+            return upperEnvelope.zip(lowerEnvelope) { upper, lower -> upper + lower / 2.0 }
         }
     }
 
@@ -76,15 +76,16 @@ public class EmpiricalModeDecomposition<BA, L> (
      * @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 = with(seriesAlgebra) {
+    private fun sift(signal: Series<Double>): SiftingResult = (seriesAlgebra) {
         val mean = findMean(signal) ?: return SiftingResult.NotEnoughExtrema
         val protoMode = signal.zip(mean) { s, m -> s - m }
         val sNumber = if (protoMode.sCondition()) 1 else 0
-        return if (maxSiftIterations == 1) SiftingResult.MaxIterationsReached(protoMode)
-        else if (sNumber >= sConditionThreshold) SiftingResult.SNumberReached(protoMode)
-        else if (relativeDifference(signal, protoMode) < siftingDelta * signal.size) {
-            SiftingResult.DeltaReached(protoMode)
-        } else siftInner(protoMode, 2, sNumber)
+        return when {
+            maxSiftIterations == 1 -> SiftingResult.MaxIterationsReached(protoMode)
+            sNumber >= sConditionThreshold -> SiftingResult.SNumberReached(protoMode)
+            relativeDifference(protoMode, signal) < siftingDelta * signal.size -> SiftingResult.DeltaReached(protoMode)
+            else -> siftInner(protoMode, 2, sNumber)
+        }
     }
 
     /**
@@ -94,14 +95,16 @@ public class EmpiricalModeDecomposition<BA, L> (
         prevMode: Series<Double>,
         iterationNumber: Int,
         sNumber: Int
-    ): SiftingResult = with(seriesAlgebra) {
+    ): SiftingResult = (seriesAlgebra) {
         val mean = findMean(prevMode) ?: return SiftingResult.SignalFlattened(prevMode)
         val mode = prevMode.zip(mean) { p, m -> p - m }
         val newSNumber = if (mode.sCondition()) sNumber + 1 else sNumber
-        return if (iterationNumber >= maxSiftIterations) SiftingResult.MaxIterationsReached(mode)
-        else if (sNumber >= sConditionThreshold) SiftingResult.SNumberReached(mode)
-        else if (relativeDifference(prevMode, mode) < siftingDelta * mode.size) SiftingResult.DeltaReached(mode)
-        else siftInner(mode, iterationNumber + 1, newSNumber)
+        return when {
+            iterationNumber >= maxSiftIterations -> SiftingResult.MaxIterationsReached(mode)
+            sNumber >= sConditionThreshold -> SiftingResult.SNumberReached(mode)
+            relativeDifference(mode, prevMode) < siftingDelta * mode.size -> SiftingResult.DeltaReached(mode)
+            else -> siftInner(mode, iterationNumber + 1, newSNumber)
+        }
     }
 
     /**
@@ -112,20 +115,17 @@ public class EmpiricalModeDecomposition<BA, L> (
      * 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 = with(seriesAlgebra) {
+     public fun decompose(signal: Series<Double>): EMDecompositionResult = (seriesAlgebra) {
         val modes = mutableListOf<Series<Double>>()
-        val mode = when(val r = sift(signal)) {
-            SiftingResult.NotEnoughExtrema ->
-                return EMDecompositionResult(EMDTerminationReason.SIGNAL_TOO_FLAT, modes)
-            is SiftingResult.Success -> r.result
-        }
-        modes.add(mode)
-
-        var residual = signal.zip(mode) { l, r -> l - r }
-        repeat(nModes - 1) {
+        var residual = signal
+        repeat(nModes) {
             val nextMode = when(val r = sift(residual)) {
                 SiftingResult.NotEnoughExtrema ->
-                    return EMDecompositionResult(EMDTerminationReason.ALL_POSSIBLE_MODES_EXTRACTED, modes)
+                    return EMDecompositionResult(
+                        if (it == 0) EMDTerminationReason.SIGNAL_TOO_FLAT
+                        else EMDTerminationReason.ALL_POSSIBLE_MODES_EXTRACTED,
+                        modes
+                    )
                 is SiftingResult.Success -> r.result
             }
             modes.add(nextMode)
@@ -149,13 +149,13 @@ public class EmpiricalModeDecomposition<BA, L> (
  * @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, BA> SeriesAlgebra<Double, *, BA, L>.empiricalModeDecomposition(
+public fun <L: Number, BA> SeriesAlgebra<Double, *, BA, L>.empiricalModeDecomposition(
     sConditionThreshold: Int = 15,
     maxSiftIterations: Int = 20,
     siftingDelta: Double = 1e-2,
     nModes: Int = 3
 ): EmpiricalModeDecomposition<BA, L>
-where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>>, L: Number = EmpiricalModeDecomposition(
+where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>> = EmpiricalModeDecomposition(
     seriesAlgebra = this,
     sConditionThreshold = sConditionThreshold,
     maxSiftIterations = maxSiftIterations,
@@ -168,22 +168,25 @@ where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>>, L: Number = Emp
  */
 private fun Series<Double>.countZeros(): Int {
     require(size >= 2) { "Expected series with at least 2 elements, but got $size elements" }
-    return fold(Pair(get(0), 0)) { acc: Pair<Double, Int>, it: Double ->
-        if (sign(acc.first) != sign(it)) Pair(it, acc.second + 1) else Pair(it, acc.second)
-    }.second
+    data class SignCounter(val prevSign: Double, val zeroCount: Int)
+
+    return fold(SignCounter(sign(get(0)), 0)) { acc: SignCounter, it: Double ->
+        val currentSign = sign(it)
+        if (acc.prevSign != currentSign) SignCounter(currentSign, acc.zeroCount + 1)
+        else SignCounter(currentSign, acc.zeroCount)
+    }.zeroCount
 }
 
 /**
  * Compute relative difference of two series.
  */
-private fun <L, BA> SeriesAlgebra<Double, *, BA, L>.relativeDifference(
+private fun <BA> SeriesAlgebra<Double, *, BA, *>.relativeDifference(
     current: Series<Double>,
     previous: Series<Double>
-):Double where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>>, L: Number {
-    return (current - previous).pow(2)
+):Double where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>> =
+    (current - previous).pow(2)
         .div(previous pow 2)
-        .fold(0.0) { acc, d -> acc + d } // to avoid unnecessary boxing, but i may be wrong
-}
+        .fold(0.0) { acc, d -> acc + d } // TODO replace with Series<>.sum() method when it's implemented
 
 private fun <T: Comparable<T>> isExtreme(prev: T, elem: T, next: T): Boolean =
     (elem > prev && elem > next) || (elem < prev && elem < next)
@@ -193,43 +196,40 @@ private fun <T: Comparable<T>> isExtreme(prev: T, elem: T, next: T): Boolean =
  */
 private fun Series<Double>.countExtrema(): Int {
     require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" }
-    return slice(1..< size - 1).asIterable().foldIndexed(0) { index, acc, it ->
-        if (isExtreme(get(index), it, get(index + 2))) acc + 1 else acc
-    }
+    return (1 .. size - 2).count { isExtreme(this[it - 1], this[it], this[it + 1]) }
 }
 
 
+/**
+ * Retrieve indices of knot points for spline interpolation matching the predicate.
+ * The first and the last points of a series are always included.
+ */
+private fun <T: Comparable<T>> Series<T>.knotPoints(predicate: (T, T, T) -> Boolean): List<Int> {
+    require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" }
+    val points = mutableListOf(0)
+    for (index in 1 .. size - 2) {
+        val left = this[index - 1]
+        val middle = this[index]
+        val right = this[index + 1]
+        if (predicate(left, middle, right)) points.add(index)
+    }
+    points.add(size - 1)
+    return points
+}
+
 /**
  * Retrieve indices of knot points used to construct an upper envelope,
  * namely maxima together with the first last point in a series.
  */
-private fun<T> Series<T>.maxima(): List<Int> where T: Comparable<T> {
-    require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" }
-    val maxima = mutableListOf(0)
-    slice(1..< size - 1).asIterable().forEachIndexed { index, it ->
-        if (it > get(index) && it > get(index + 2)) { // weird offset, is there a way to do it better?
-            maxima.add(index + 1)
-        }
-    }
-    maxima.add(size - 1)
-    return maxima
-}
+private fun <T: Comparable<T>> Series<T>.paddedMaxima(): List<Int> =
+    knotPoints { left, middle, right -> (middle > left && middle > right) }
 
 /**
  * Retrieve indices of knot points used to construct a lower envelope,
  * namely minima together with the first last point in a series.
  */
-private fun<T> Series<T>.minima(): List<Int> where T: Comparable<T> {
-    require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" }
-    val minima = mutableListOf(0)
-    slice(1..< size - 1).asIterable().forEachIndexed { index, it ->
-        if (it < get(index) && it < get(index + 2)) { // weird offset, is there a way to do it better?
-            minima.add(index + 1)
-        }
-    }
-    minima.add(size - 1)
-    return minima
-}
+private fun <T: Comparable<T>> Series<T>.paddedMinima(): List<Int> =
+    knotPoints { left, middle, right -> (middle < left && middle < right) }
 
 /**
  * Check whether the numbers of zeroes and extrema of a series differ by no more than 1.