diff --git a/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt b/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt new file mode 100644 index 000000000..11529b904 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt @@ -0,0 +1,40 @@ +/* + * 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.structures.* +import space.kscience.plotly.* +import space.kscience.plotly.models.Scatter +import kotlin.math.sin + +fun main(): Unit = with(Double.seriesAlgebra()) { + 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, + nModes = 4 + ).decompose(signal) + println("EMD: ${emd.modes.size} modes extracted, terminated because ${emd.terminatedBecause}") + + fun Plot.series(name: String, buffer: Buffer, block: Scatter.() -> Unit = {}) { + this.scatter { + this.name = name + this.x.numbers = buffer.offsetIndices + 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) +} \ No newline at end of file diff --git a/kmath-stat/build.gradle.kts b/kmath-stat/build.gradle.kts index 8116a8925..50fcbac00 100644 --- a/kmath-stat/build.gradle.kts +++ b/kmath-stat/build.gradle.kts @@ -14,6 +14,7 @@ kotlin.sourceSets { dependencies { api(projects.kmathCoroutines) //implementation(spclibs.atomicfu) + api(project(":kmath-functions")) } } diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt new file mode 100644 index 000000000..d7f53179a --- /dev/null +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt @@ -0,0 +1,307 @@ +/* + * 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.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 + +/** + * 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 ( + private val seriesAlgebra: SeriesAlgebra, + private val sConditionThreshold: Int = 15, + private val maxSiftIterations: Int = 20, + private val siftingDelta: Double = 1e-2, + private val nModes: Int = 6 +) where BA: BufferAlgebra, BA: RingOps>, L : Number { + + /** + * Take a signal, construct an upper and a lower envelope, 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. + */ + private fun findMean(signal: Series): Series? = with(seriesAlgebra) { + val interpolator = SplineInterpolator(Float64Field) + fun generateEnvelope(extrema: List): Series { + val envelopeFunction = interpolator.interpolate( + Buffer(extrema.size) { signal.labels[extrema[it]].toDouble() }, + Buffer(extrema.size) { signal[extrema[it]] } + ) + return signal.mapWithLabel { _, label -> + // For some reason PolynomialInterpolator is exclusive and the right boundary + // TODO Notify interpolator authors + envelopeFunction(label.toDouble()) ?: signal.last() + // need to make the interpolator yield values outside boundaries? + } + } + val maxima = signal.maxima() + val minima = signal.minima() + 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 } + } + } + + /** + * 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): SiftingResult = with(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) + } + + /** + * Compute a single iteration of the sifting process. + */ + private fun siftInner( + prevMode: Series, + iterationNumber: Int, + sNumber: Int + ): SiftingResult = with(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) + } + + /** + * 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): EMDecompositionResult = with(seriesAlgebra) { + val modes = mutableListOf>() + 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) { + val nextMode = when(val r = sift(residual)) { + SiftingResult.NotEnoughExtrema -> + return EMDecompositionResult(EMDTerminationReason.ALL_POSSIBLE_MODES_EXTRACTED, modes) + is SiftingResult.Success -> r.result + } + modes.add(nextMode) + residual = residual.zip(nextMode) { l, r -> l - r } + } + return EMDecompositionResult(EMDTerminationReason.MAX_MODES_REACHED, modes) + } +} + +/** + * 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 SeriesAlgebra.empiricalModeDecomposition( + sConditionThreshold: Int = 15, + maxSiftIterations: Int = 20, + siftingDelta: Double = 1e-2, + nModes: Int = 3 +): EmpiricalModeDecomposition +where BA: BufferAlgebra, BA: RingOps>, L: Number = EmpiricalModeDecomposition( + seriesAlgebra = this, + sConditionThreshold = sConditionThreshold, + maxSiftIterations = maxSiftIterations, + siftingDelta = siftingDelta, + nModes = nModes +) + +/** + * Brute force count all zeros in the series. + */ +private fun Series.countZeros(): Int { + require(size >= 2) { "Expected series with at least 2 elements, but got $size elements" } + return fold(Pair(get(0), 0)) { acc: Pair, it: Double -> + if (sign(acc.first) != sign(it)) Pair(it, acc.second + 1) else Pair(it, acc.second) + }.second +} + +/** + * Compute relative difference of two series. + */ +private fun SeriesAlgebra.relativeDifference( + current: Series, + previous: Series +):Double where BA: BufferAlgebra, BA: RingOps>, L: Number { + return (current - previous).pow(2) + .div(previous pow 2) + .fold(0.0) { acc, d -> acc + d } // to avoid unnecessary boxing, but i may be wrong +} + +private fun > isExtreme(prev: T, elem: T, next: T): Boolean = + (elem > prev && elem > next) || (elem < prev && elem < next) + +/** + * Brute force count all extrema of a series. + */ +private fun Series.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 + } +} + + +/** + * Retrieve indices of knot points used to construct an upper envelope, + * namely maxima together with the first last point in a series. + */ +private fun Series.maxima(): List where T: Comparable { + 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 +} + +/** + * Retrieve indices of knot points used to construct a lower envelope, + * namely minima together with the first last point in a series. + */ +private fun Series.minima(): List where T: Comparable { + 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 +} + +/** + * 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.sCondition(): Boolean = (countExtrema() - countZeros()) in -1..1 + +internal sealed interface SiftingResult { + + /** + * Represents a condition when a mode has been successfully + * extracted in a sifting process. + */ + open class Success(val result: Series): 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) : Success(result) + + /** + * Returned when sifting process has been terminated due to the + * S-number condition being reached. + */ + class SNumberReached(result: Series) : Success(result) + + /** + * Returned when sifting process has been terminated due to the + * delta condition (Cauchy criterion) being reached. + */ + class DeltaReached(result: Series) : Success(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): Success(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( + val terminatedBecause: EMDTerminationReason, + val modes: List> +) \ No newline at end of file