From e0125f0b26295ae14b83a97a20b3b49078af029b Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Wed, 24 Apr 2024 16:56:52 +0300 Subject: [PATCH 01/10] mapWithLabel fix --- .../kotlin/space/kscience/kmath/series/SeriesAlgebra.kt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt index e47e4f69f..92be97572 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt @@ -169,7 +169,7 @@ public open class SeriesAlgebra, out BA : BufferAlgebra public inline fun Buffer.mapWithLabel(crossinline transform: A.(arg: T, label: L) -> T): Series { 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) } -- 2.34.1 From 8fa42450b242a3e6a5a9da8c452057ece1d03cbb Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Wed, 24 Apr 2024 16:58:35 +0300 Subject: [PATCH 02/10] Empirical Mode Decomposition for Double series --- .../space/kscience/kmath/series/emdDemo.kt | 40 +++ kmath-stat/build.gradle.kts | 1 + .../series/EmpiricalModeDecomposition.kt | 307 ++++++++++++++++++ 3 files changed, 348 insertions(+) create mode 100644 examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt create mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt 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 -- 2.34.1 From 4f2ebb17ff285891c072f9f1597d22975b685ade Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Thu, 25 Apr 2024 20:55:54 +0300 Subject: [PATCH 03/10] Submit corrections --- .../space/kscience/kmath/series/emdDemo.kt | 3 +- .../series/EmpiricalModeDecomposition.kt | 132 +++++++++--------- 2 files changed, 68 insertions(+), 67 deletions(-) diff --git a/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt b/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt index 11529b904..c69b163cb 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt +++ b/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) 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 index d7f53179a..179f4ae43 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt +++ b/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 ( +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 { +) where BA: BufferAlgebra, BA: RingOps> { /** - * 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): Series? = with(seriesAlgebra) { + private fun findMean(signal: Series): Series? = (seriesAlgebra) { val interpolator = SplineInterpolator(Float64Field) fun generateEnvelope(extrema: List): Series { val envelopeFunction = interpolator.interpolate( @@ -59,12 +59,12 @@ public class EmpiricalModeDecomposition ( // 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 ( * @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) { + private fun sift(signal: Series): 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 ( prevMode: Series, 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 ( * 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) { + public fun decompose(signal: Series): EMDecompositionResult = (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) { + 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 ( * @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( +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( +where BA: BufferAlgebra, BA: RingOps> = EmpiricalModeDecomposition( seriesAlgebra = this, sConditionThreshold = sConditionThreshold, maxSiftIterations = maxSiftIterations, @@ -168,22 +168,25 @@ where BA: BufferAlgebra, BA: RingOps>, L: Number = Emp */ 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 + 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 SeriesAlgebra.relativeDifference( +private fun SeriesAlgebra.relativeDifference( current: Series, previous: Series -):Double where BA: BufferAlgebra, BA: RingOps>, L: Number { - return (current - previous).pow(2) +):Double where BA: BufferAlgebra, BA: RingOps> = + (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 > isExtreme(prev: T, elem: T, next: T): Boolean = (elem > prev && elem > next) || (elem < prev && elem < next) @@ -193,43 +196,40 @@ private fun > isExtreme(prev: T, elem: T, next: T): Boolean = */ 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 - } + 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 > Series.knotPoints(predicate: (T, T, T) -> Boolean): List { + 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 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 -} +private fun > Series.paddedMaxima(): List = + 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 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 -} +private fun > Series.paddedMinima(): List = + 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. -- 2.34.1 From 1c5113da29e9ab35fa4d63b936cd8a9e5eef582d Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Thu, 25 Apr 2024 22:50:07 +0300 Subject: [PATCH 04/10] Allow tail call optimization --- .../space/kscience/kmath/series/EmpiricalModeDecomposition.kt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 index 179f4ae43..d1f4f43d2 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt @@ -91,7 +91,7 @@ public class EmpiricalModeDecomposition ( /** * Compute a single iteration of the sifting process. */ - private fun siftInner( + private tailrec fun siftInner( prevMode: Series, iterationNumber: Int, sNumber: Int -- 2.34.1 From 40ffa8d6fc5380a4f2911178f48ab07c4084997f Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Sun, 28 Apr 2024 18:36:54 +0300 Subject: [PATCH 05/10] Remove unnecessary code from outer sift() function --- .../kmath/series/EmpiricalModeDecomposition.kt | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) 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 index d1f4f43d2..467ea3201 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt @@ -76,17 +76,7 @@ public class EmpiricalModeDecomposition ( * @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 = (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 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) - } - } + private fun sift(signal: Series): SiftingResult = siftInner(signal, 1, 0) /** * Compute a single iteration of the sifting process. @@ -96,7 +86,9 @@ public class EmpiricalModeDecomposition ( iterationNumber: Int, sNumber: Int ): SiftingResult = (seriesAlgebra) { - val mean = findMean(prevMode) ?: return SiftingResult.SignalFlattened(prevMode) + 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 (mode.sCondition()) sNumber + 1 else sNumber return when { -- 2.34.1 From 597c27e615de773594623bef6a283fe0420a28f4 Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Fri, 3 May 2024 17:04:08 +0300 Subject: [PATCH 06/10] Improved peak finding algorithm --- .../kscience/kmath/series/peakFinding.kt | 89 +++++++++++++++++++ 1 file changed, 89 insertions(+) create mode 100644 kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/peakFinding.kt diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/peakFinding.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/peakFinding.kt new file mode 100644 index 000000000..cd340b044 --- /dev/null +++ b/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 > Series.peaks( + plateauEdgePolicy: PlateauEdgePolicy = PlateauEdgePolicy.AVERAGE +): List = findPeaks(plateauEdgePolicy, { other -> this > other }, { other -> this >= other }) + +public fun > Series.troughs( + plateauEdgePolicy: PlateauEdgePolicy = PlateauEdgePolicy.AVERAGE +): List = findPeaks(plateauEdgePolicy, { other -> this < other }, { other -> this <= other }) + + +private fun > Series.findPeaks( + plateauPolicy: PlateauEdgePolicy = PlateauEdgePolicy.AVERAGE, + cmpStrong: T.(T) -> Boolean, + cmpWeak: T.(T) -> Boolean +): List { + 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() + 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 > Series.peaksWithPlateau(cmpStrong: T.(T) -> Boolean): List { + val peaks = mutableListOf() + 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 +} \ No newline at end of file -- 2.34.1 From b8bf56bc42264d2ef99b95fcf5d67edc8636ca93 Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Fri, 3 May 2024 22:07:21 +0300 Subject: [PATCH 07/10] Extrema padding (experimental) --- .../series/EmpiricalModeDecomposition.kt | 69 ++++++++----------- 1 file changed, 27 insertions(+), 42 deletions(-) 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 index 467ea3201..301247698 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt @@ -11,7 +11,7 @@ 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.asBuffer import kotlin.math.sign /** @@ -47,24 +47,40 @@ public class EmpiricalModeDecomposition ( */ private fun findMean(signal: Series): Series? = (seriesAlgebra) { val interpolator = SplineInterpolator(Float64Field) - fun generateEnvelope(extrema: List): Series { + fun generateEnvelope(extrema: List, paddedExtremeValues: DoubleArray): Series { val envelopeFunction = interpolator.interpolate( Buffer(extrema.size) { signal.labels[extrema[it]].toDouble() }, - Buffer(extrema.size) { signal[extrema[it]] } + paddedExtremeValues.asBuffer() ) return signal.mapWithLabel { _, label -> // For some reason PolynomialInterpolator is exclusive and the right boundary // TODO Notify interpolator authors - envelopeFunction(label.toDouble()) ?: signal.last() + envelopeFunction(label.toDouble()) ?: paddedExtremeValues.last() // need to make the interpolator yield values outside boundaries? } } - 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.zip(lowerEnvelope) { upper, lower -> upper + lower / 2.0 } + // 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]] } + if (maxValues[0] < maxValues[1]) { + maxValues[0] = maxValues[1] + } + if (maxValues.last() < maxValues[maxValues.lastIndex - 1]) { + maxValues[maxValues.lastIndex] = maxValues[maxValues.lastIndex - 1] + } + + val minima = listOf(0) + signal.troughs() + (signal.size - 1) + val minValues = DoubleArray(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] + } + 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 } } } @@ -186,43 +202,12 @@ private fun > isExtreme(prev: T, elem: T, next: T): Boolean = /** * Brute force count all extrema of a series. */ +@Deprecated("Does not match the algorithm currently in use.") private fun Series.countExtrema(): Int { require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" } 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 > Series.knotPoints(predicate: (T, T, T) -> Boolean): List { - 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 > Series.paddedMaxima(): List = - 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 > Series.paddedMinima(): List = - 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. * This is a necessary condition of an empirical mode. -- 2.34.1 From 47a9bf0e9a3bb64dfea92f16f82485b9f0468dd6 Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Fri, 3 May 2024 22:37:59 +0300 Subject: [PATCH 08/10] Suboptimal implementation for countExtrema --- .../kscience/kmath/series/EmpiricalModeDecomposition.kt | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) 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 index 301247698..ea8a0aa62 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt @@ -196,16 +196,12 @@ private fun SeriesAlgebra.relativeDifference( .div(previous pow 2) .fold(0.0) { acc, d -> acc + d } // TODO replace with Series<>.sum() method when it's implemented -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. */ -@Deprecated("Does not match the algorithm currently in use.") private fun Series.countExtrema(): Int { require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" } - return (1 .. size - 2).count { isExtreme(this[it - 1], this[it], this[it + 1]) } + return peaks().size + troughs().size } /** -- 2.34.1 From 2e084edc9b4fb7e22d06c6050bd74538dcc2c5b6 Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Mon, 20 May 2024 17:02:45 +0300 Subject: [PATCH 09/10] Implement generic algorithm --- .../space/kscience/kmath/series/emdDemo.kt | 10 +- .../series/EmpiricalModeDecomposition.kt | 118 ++++++++++-------- 2 files changed, 70 insertions(+), 58 deletions(-) diff --git a/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt b/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt index c69b163cb..a46f3768a 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/series/emdDemo.kt +++ b/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, block: Scatter.() -> Unit = {}) { this.scatter { this.name = name - this.x.numbers = buffer.offsetIndices + this.x.numbers = buffer.labels this.y.doubles = buffer.toDoubleArray() block() } 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 index ea8a0aa62..700b15757 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt +++ b/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 ( - private val seriesAlgebra: SeriesAlgebra, +public class EmpiricalModeDecomposition, A: Field, BA, L: T> ( + private val seriesAlgebra: SeriesAlgebra, 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, BA: RingOps> { +) where BA: BufferAlgebra, BA: FieldOps> { /** * Take a signal, construct an upper and a lower envelopes, find the mean value of two, @@ -45,37 +42,38 @@ public class EmpiricalModeDecomposition ( * @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): Series? = (seriesAlgebra) { - val interpolator = SplineInterpolator(Float64Field) - fun generateEnvelope(extrema: List, paddedExtremeValues: DoubleArray): Series { + private fun findMean(signal: Series): Series? = (seriesAlgebra) { + val interpolator = SplineInterpolator(elementAlgebra) + val makeBuffer = elementAlgebra.bufferFactory + fun generateEnvelope(extrema: List, paddedExtremeValues: Buffer): Series { 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 ( * @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 = siftInner(signal, 1, 0) + private fun sift(signal: Series): SiftingResult = siftInner(signal, 1, 0) /** * Compute a single iteration of the sifting process. */ private tailrec fun siftInner( - prevMode: Series, + prevMode: Series, iterationNumber: Int, sNumber: Int ): SiftingResult = (seriesAlgebra) { @@ -106,11 +104,12 @@ public class EmpiricalModeDecomposition ( 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 ( * 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 = (seriesAlgebra) { - val modes = mutableListOf>() + public fun decompose(signal: Series): EMDecompositionResult = (seriesAlgebra) { + val modes = mutableListOf>() var residual = signal repeat(nModes) { val nextMode = when(val r = sift(residual)) { @@ -132,14 +131,15 @@ public class EmpiricalModeDecomposition ( 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) // 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 ( * @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( +public fun , L: T, A: Field, BA> SeriesAlgebra.empiricalModeDecomposition( sConditionThreshold: Int = 15, maxSiftIterations: Int = 20, - siftingDelta: Double = 1e-2, + siftingDelta: T, nModes: Int = 3 -): EmpiricalModeDecomposition -where BA: BufferAlgebra, BA: RingOps> = EmpiricalModeDecomposition( +): EmpiricalModeDecomposition +where BA: BufferAlgebra, BA: FieldOps> = EmpiricalModeDecomposition( seriesAlgebra = this, sConditionThreshold = sConditionThreshold, maxSiftIterations = maxSiftIterations, @@ -174,12 +174,15 @@ where BA: BufferAlgebra, BA: RingOps> = EmpiricalModeD /** * 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" } - data class SignCounter(val prevSign: Double, val zeroCount: Int) +private fun , A: Ring, BA> SeriesAlgebra.countZeros( + signal: Series +): Int where BA: BufferAlgebra, BA: FieldOps> { + 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.countZeros(): Int { /** * Compute relative difference of two series. */ -private fun SeriesAlgebra.relativeDifference( - current: Series, - previous: Series -):Double where BA: BufferAlgebra, BA: RingOps> = - (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 , BA> SeriesAlgebra.relativeDifference( + current: Series, + previous: Series +): T where BA: BufferAlgebra, BA: FieldOps> = (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.countExtrema(): Int { +private fun > Series.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.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.sCondition(): Boolean = (countExtrema() - countZeros()) in -1..1 +private fun , A: Ring, BA> SeriesAlgebra.sCondition( + signal: Series +): Boolean where BA: BufferAlgebra, BA: FieldOps> = + (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): SiftingResult + 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) + 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) + 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) + 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) + class MaxIterationsReached(result: Series): Success(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( val terminatedBecause: EMDTerminationReason, - val modes: List> + val modes: List>, + val residual: Series ) \ No newline at end of file -- 2.34.1 From 6f98c3fbe9d0a50f195507d409a0c2509a3a3ac8 Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Mon, 20 May 2024 17:38:54 +0300 Subject: [PATCH 10/10] Tests and peak finding examples --- .../kscience/kmath/series/peakFinding.kt | 85 +++++++++++++++++++ .../series/EmpiricalModeDecomposition.kt | 4 +- .../space/kscience/kmath/series/TestEmd.kt | 66 ++++++++++++++ .../kscience/kmath/series/TestPeakFinding.kt | 39 +++++++++ 4 files changed, 192 insertions(+), 2 deletions(-) create mode 100644 examples/src/main/kotlin/space/kscience/kmath/series/peakFinding.kt create mode 100644 kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestEmd.kt create mode 100644 kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestPeakFinding.kt diff --git a/examples/src/main/kotlin/space/kscience/kmath/series/peakFinding.kt b/examples/src/main/kotlin/space/kscience/kmath/series/peakFinding.kt new file mode 100644 index 000000000..ba5e334eb --- /dev/null +++ b/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, 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, 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) +} \ No newline at end of file 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 index 700b15757..67a45c19c 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/EmpiricalModeDecomposition.kt @@ -174,7 +174,7 @@ where BA: BufferAlgebra, BA: FieldOps> = EmpiricalModeDecomposit /** * Brute force count all zeros in the series. */ -private fun , A: Ring, BA> SeriesAlgebra.countZeros( +internal fun , A: Ring, BA> SeriesAlgebra.countZeros( signal: Series ): Int where BA: BufferAlgebra, BA: FieldOps> { require(signal.size >= 2) { "Expected series with at least 2 elements, but got ${signal.size} elements" } @@ -203,7 +203,7 @@ private fun , BA> SeriesAlgebra.relativeDifference( /** * Brute force count all extrema of a series. */ -private fun > Series.countExtrema(): Int { +internal fun > Series.countExtrema(): Int { require(size >= 3) { "Expected series with at least 3 elements, but got $size elements" } return peaks().size + troughs().size } diff --git a/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestEmd.kt b/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestEmd.kt new file mode 100644 index 000000000..7db7e5364 --- /dev/null +++ b/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) + } +} \ No newline at end of file diff --git a/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestPeakFinding.kt b/kmath-stat/src/commonTest/kotlin/space/kscience/kmath/series/TestPeakFinding.kt new file mode 100644 index 000000000..5edbd6f58 --- /dev/null +++ b/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) + } +} \ No newline at end of file -- 2.34.1