Part 3

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

@@ -0,0 +1,45 @@
+/*
+ * 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.structures.*
+import space.kscience.kmath.operations.invoke
+import space.kscience.plotly.*
+import space.kscience.plotly.models.Scatter
+import kotlin.math.sin
+
+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}")
+
+    fun Plot.series(name: String, buffer: Buffer<Double>, block: Scatter.() -> Unit = {}) {
+        this.scatter {
+            this.name = name
+            this.x.numbers = buffer.labels
+            this.y.doubles = buffer.toDoubleArray()
+            block()
+        }
+    }
+
+    Plotly.plot {
+        series("Signal", signal)
+        emd.modes.forEachIndexed { index, it ->
+            series("Mode ${index+1}", it)
+        }
+    }.makeFile(resourceLocation = ResourceLocation.REMOTE)
+}

kmath-stat/build.gradle.kts

@@ -14,6 +14,7 @@ kotlin.sourceSets {
         dependencies {
             api(projects.kmathCoroutines)
             //implementation(spclibs.atomicfu)
+            api(project(":kmath-functions"))
         }
     }
 

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

@@ -0,0 +1,288 @@
+/*
+ * 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.interpolation.SplineInterpolator
+import space.kscience.kmath.interpolation.interpolate
+import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.last
+
+/**
+ * Empirical mode decomposition of a signal represented as a [Series].
+ *
+ * @param seriesAlgebra context to perform operations in.
+ * @param maxSiftIterations number of iterations in the mode extraction (sifting) process after which
+ * the result is returned even if no other conditions are met.
+ * @param sConditionThreshold one of the possible criteria for sifting process termination:
+ * how many times in a row should a proto-mode satisfy the s-condition (the number of zeros differs from
+ * the number of extrema no more than by 1) to be considered an empirical mode.
+ * @param siftingDelta one of the possible criteria for sifting process termination: if relative difference of
+ * two proto-modes obtained in a sequence is less that this number, the last proto-mode is considered
+ * an empirical mode and returned.
+ * @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<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: T,
+    private val nModes: Int = 6
+) 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,
+     * represented as a series.
+     * @param signal Signal to compute on
+     *
+     * @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<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(
+                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) ?: 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 = makeBuffer(maxima.size) { signal[maxima[it]] }
+        if (maxValues[0] < maxValues[1]) {
+            maxValues[0] = maxValues[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 = makeBuffer(minima.size) { signal[minima[it]] }
+        if (minValues[0] > minValues[1]) {
+            minValues[0] = minValues[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)
+            val lowerEnvelope = generateEnvelope(minima, minValues)
+            return (upperEnvelope + lowerEnvelope).map { it * 0.5 }
+        }
+    }
+
+    /**
+     * Extract a single empirical mode from a signal. This process is called sifting, hence the name.
+     * @param signal Signal to extract a mode from. The first mode is extracted from the initial signal,
+     * subsequent modes are extracted from the residuals between the signal and all previous modes.
+     *
+     * @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<T>): SiftingResult = siftInner(signal, 1, 0)
+
+    /**
+     * Compute a single iteration of the sifting process.
+     */
+    private tailrec fun siftInner(
+        prevMode: Series<T>,
+        iterationNumber: Int,
+        sNumber: Int
+    ): SiftingResult = (seriesAlgebra) {
+        val mean = findMean(prevMode) ?:
+            return if (iterationNumber == 1) SiftingResult.NotEnoughExtrema
+            else SiftingResult.SignalFlattened(prevMode)
+        val mode = prevMode.zip(mean) { p, m -> p - m }
+        val newSNumber = if (sCondition(mode)) sNumber + 1 else sNumber
+        return when {
+            iterationNumber >= maxSiftIterations -> SiftingResult.MaxIterationsReached(mode)
+            sNumber >= sConditionThreshold -> SiftingResult.SNumberReached(mode)
+            relativeDifference(mode, prevMode) < (elementAlgebra) { siftingDelta * mode.size } ->
+                SiftingResult.DeltaReached(mode)
+            else -> siftInner(mode, iterationNumber + 1, newSNumber)
+        }
+    }
+
+    /**
+     * Extract [nModes] empirical modes from a signal represented by a time series.
+     * @param signal Signal to decompose.
+     * @return [EMDecompositionResult] with an appropriate reason for algorithm termination
+     * (see [EMDTerminationReason] for possible reasons).
+     * 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<T>): EMDecompositionResult<T> = (seriesAlgebra) {
+        val modes = mutableListOf<Series<T>>()
+        var residual = signal
+        repeat(nModes) {
+            val nextMode = when(val r = sift(residual)) {
+                SiftingResult.NotEnoughExtrema ->
+                    return EMDecompositionResult(
+                        if (it == 0) EMDTerminationReason.SIGNAL_TOO_FLAT
+                        else EMDTerminationReason.ALL_POSSIBLE_MODES_EXTRACTED,
+                        modes,
+                        residual
+                    )
+                is SiftingResult.Success<*> -> r.result
+            }
+            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, residual)
+    }
+}
+
+/**
+ * Shortcut to retrieve a decomposition factory from a [SeriesAlgebra] scope.
+ * @receiver scope to perform operations in.
+ * @param maxSiftIterations number of iterations in the mode extraction (sifting) process after which
+ * the result is returned even if no other conditions are met.
+ * @param sConditionThreshold one of the possible criteria for sifting process termination:
+ * how many times in a row should a proto-mode satisfy the s-condition (the number of zeros differs from
+ * the number of extrema no more than by 1) to be considered an empirical mode.
+ * @param siftingDelta one of the possible criteria for sifting process termination: if relative difference of
+ * two proto-modes obtained in a sequence is less that this number, the last proto-mode is considered
+ * an empirical mode and returned.
+ * @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 <T: Comparable<T>, L: T, A: Field<T>, BA> SeriesAlgebra<T, A, BA, L>.empiricalModeDecomposition(
+    sConditionThreshold: Int = 15,
+    maxSiftIterations: Int = 20,
+    siftingDelta: T,
+    nModes: Int = 3
+): EmpiricalModeDecomposition<T, A, BA, L>
+where BA: BufferAlgebra<T, A>, BA: FieldOps<Buffer<T>> = EmpiricalModeDecomposition(
+    seriesAlgebra = this,
+    sConditionThreshold = sConditionThreshold,
+    maxSiftIterations = maxSiftIterations,
+    siftingDelta = siftingDelta,
+    nModes = nModes
+)
+
+/**
+ * Brute force count all zeros in the series.
+ */
+private 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 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
+}
+
+/**
+ * Compute relative difference of two series.
+ */
+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 <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
+}
+
+/**
+ * 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 <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 {
+
+    /**
+     * Represents a condition when a mode has been successfully
+     * extracted in a sifting process.
+     */
+    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<T>(result: Series<T>) : Success<T>(result)
+
+    /**
+     * Returned when sifting process has been terminated due to the
+     * S-number condition being reached.
+     */
+    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<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<T>(result: Series<T>): Success<T>(result)
+
+    /**
+     * Returned when the submitted signal has not enough extrema to build envelopes,
+     * i.e. when [SignalFlattened] condition has already been reached before the first sifting iteration.
+     */
+    data object NotEnoughExtrema : SiftingResult
+
+}
+
+public enum class EMDTerminationReason {
+    /**
+     * Returned when the signal is too flat, i.e. there are too few extrema
+     * to build envelopes necessary to extract modes.
+     */
+    SIGNAL_TOO_FLAT,
+
+    /**
+     * Returned when there has been extracted as many modes as
+     * specified when creating the instance of [EmpiricalModeDecomposition]
+     */
+    MAX_MODES_REACHED,
+
+    /**
+     * Returned when the algorithm terminated after finding impossible
+     * to extract more modes from the signal and the maximum number
+     * of modes (specified when creating an instance of [EmpiricalModeDecomposition])
+     * has not yet been reached.
+     */
+    ALL_POSSIBLE_MODES_EXTRACTED
+}
+
+public data class EMDecompositionResult<T>(
+    val terminatedBecause: EMDTerminationReason,
+    val modes: List<Series<T>>,
+    val residual: Series<T>
+)

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

@@ -169,7 +169,7 @@ public open class SeriesAlgebra<T, out A : Ring<T>, out BA : BufferAlgebra<T, A>
     public inline fun Buffer<T>.mapWithLabel(crossinline transform: A.(arg: T, label: L) -> T): Series<T> {
         val labels = labels
         val buf = elementAlgebra.bufferFactory(size) {
-            elementAlgebra.transform(getByOffset(it), labels[it])
+            elementAlgebra.transform(get(it), labels[it])
         }
         return buf.moveTo(offsetIndices.first)
     }

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

@@ -0,0 +1,89 @@
+/*
+ * 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
+
+
+public enum class PlateauEdgePolicy  {
+    /**
+     * Midpoints of plateau are returned, edges not belonging to a plateau are ignored.
+     *
+     * A midpoint is the index closest to the average of indices of the left and right edges
+     * of the plateau:
+     *
+     * `val midpoint = ((leftEdge + rightEdge) / 2).toInt`
+     */
+     AVERAGE,
+
+    /**
+     * Both left and right edges are returned.
+     */
+    KEEP_ALL_EDGES,
+
+    /**
+     * Only right edges are returned.
+     */
+    KEEP_RIGHT_EDGES,
+
+    /**
+     * Only left edges are returned.
+     */
+    KEEP_LEFT_EDGES,
+
+    /**
+     * Ignore plateau, only peaks (troughs) with values strictly greater (less)
+     * than values of the adjacent points are returned.
+     */
+    IGNORE
+}
+
+
+public fun <T: Comparable<T>> Series<T>.peaks(
+    plateauEdgePolicy: PlateauEdgePolicy = PlateauEdgePolicy.AVERAGE
+): List<Int> = findPeaks(plateauEdgePolicy, { other -> this > other },  { other -> this >= other })
+
+public fun <T: Comparable<T>> Series<T>.troughs(
+    plateauEdgePolicy: PlateauEdgePolicy = PlateauEdgePolicy.AVERAGE
+): List<Int> = findPeaks(plateauEdgePolicy, { other -> this < other }, { other -> this <= other })
+
+
+private fun <T: Comparable<T>> Series<T>.findPeaks(
+    plateauPolicy: PlateauEdgePolicy = PlateauEdgePolicy.AVERAGE,
+    cmpStrong: T.(T) -> Boolean,
+    cmpWeak: T.(T) -> Boolean
+): List<Int> {
+    require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" }
+    if (plateauPolicy == PlateauEdgePolicy.AVERAGE) return peaksWithPlateau(cmpStrong)
+    fun peakCriterion(left: T, middle: T, right: T): Boolean = when(plateauPolicy) {
+        PlateauEdgePolicy.KEEP_LEFT_EDGES -> middle.cmpStrong(left) && middle.cmpWeak(right)
+        PlateauEdgePolicy.KEEP_RIGHT_EDGES -> middle.cmpWeak(left) && middle.cmpStrong(right)
+        PlateauEdgePolicy.KEEP_ALL_EDGES ->
+            (middle.cmpStrong(left) && middle.cmpWeak(right)) || (middle.cmpWeak(left) && middle.cmpStrong(right))
+        else -> middle.cmpStrong(right) && middle.cmpStrong(left)
+    }
+    val indices = mutableListOf<Int>()
+    for (index in 1 .. size - 2) {
+        val left = this[index - 1]
+        val middle = this[index]
+        val right = this[index + 1]
+        if (peakCriterion(left, middle, right)) indices.add(index)
+    }
+    return indices
+}
+
+
+private fun <T: Comparable<T>> Series<T>.peaksWithPlateau(cmpStrong: T.(T) -> Boolean): List<Int> {
+    val peaks = mutableListOf<Int>()
+    tailrec fun peaksPlateauInner(index: Int) {
+        val nextUnequal = (index + 1 ..< size).firstOrNull { this[it] != this[index] } ?: (size - 1)
+        val newIndex = if (this[index].cmpStrong(this[index - 1]) && this[index].cmpStrong(this[nextUnequal])) {
+            peaks.add((index + nextUnequal) / 2)
+            nextUnequal
+        } else index + 1
+        if (newIndex < size - 1) peaksPlateauInner(newIndex)
+    }
+    peaksPlateauInner(1)
+    return peaks
+}