From 4f2ebb17ff285891c072f9f1597d22975b685ade Mon Sep 17 00:00:00 2001 From: Igor Dunaev Date: Thu, 25 Apr 2024 20:55:54 +0300 Subject: [PATCH] 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.