teldufalsari/kmath
2
0
Fork 0
forked from kscience/kmath
Code Pull Requests 1 Activity

Part 3

Browse Source
This commit is contained in:
Igor Dunaev 2024-05-04 23:34:49 +03:00
parent d2c8423b6f
commit 5efd5e8304
2 changed files with 31 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -5,20 +5,24 @@
package space.kscience.kmath.series 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.structures.*
import space.kscience.kmath.operations.invoke import space.kscience.kmath.operations.invoke
import space.kscience.plotly.* import space.kscience.plotly.*
import space.kscience.plotly.models.Scatter import space.kscience.plotly.models.Scatter
import kotlin.math.sin 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) { val signal = DoubleArray(800) {
sin(it.toDouble() / 10.0) + 3.5 * sin(it.toDouble() / 60.0) sin(it.toDouble() / 10.0) + 3.5 * sin(it.toDouble() / 60.0)
}.asBuffer().moveTo(0) }.asBuffer().moveTo(0)
val emd = empiricalModeDecomposition( val emd = empiricalModeDecomposition(
sConditionThreshold = 1, sConditionThreshold = 1,
maxSiftIterations = 15, maxSiftIterations = 15,
siftingDelta = 1e-2,
nModes = 4 nModes = 4
).decompose(signal) ).decompose(signal)
println("EMD: ${emd.modes.size} modes extracted, terminated because ${emd.terminatedBecause}") 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<Double>, block: Scatter.() -> Unit = {}) { fun Plot.series(name: String, buffer: Buffer<Double>, block: Scatter.() -> Unit = {}) {
this.scatter { this.scatter {
this.name = name this.name = name
this.x.numbers = buffer.offsetIndices this.x.numbers = buffer.labels
this.y.doubles = buffer.toDoubleArray() this.y.doubles = buffer.toDoubleArray()
block() block()
} }

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

@ -26,13 +26,13 @@ import space.kscience.kmath.structures.last
* @param nModes how many modes should be extracted at most. The algorithm may return fewer modes if it was not * @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. * possible to extract more modes from the signal.
*/ */
public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> ( public class EmpiricalModeDecomposition<T: Comparable<T>, A: Field<T>, BA, L: T> (
private val seriesAlgebra: SeriesAlgebra<Double, A, BA, L>, private val seriesAlgebra: SeriesAlgebra<T, A, BA, L>,
private val sConditionThreshold: Int = 15, private val sConditionThreshold: Int = 15,
private val maxSiftIterations: Int = 20, private val maxSiftIterations: Int = 20,
private val siftingDelta: Double = 1e-2, private val siftingDelta: T,
private val nModes: Int = 6 private val nModes: Int = 6
) where BA: BufferAlgebra<Double, A>, BA: FieldOps<Buffer<Double>> { ) 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, * Take a signal, construct an upper and a lower envelopes, find the mean value of two,
@ -42,18 +42,18 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
* @return mean [Series] or `null`. `null` is returned in case * @return mean [Series] or `null`. `null` is returned in case
* the signal does not have enough extrema to construct envelopes. * the signal does not have enough extrema to construct envelopes.
*/ */
private fun findMean(signal: Series<Double>): Series<Double>? = (seriesAlgebra) { private fun findMean(signal: Series<T>): Series<T>? = (seriesAlgebra) {
val interpolator = SplineInterpolator(elementAlgebra) val interpolator = SplineInterpolator(elementAlgebra)
val makeBuffer = elementAlgebra.bufferFactory val makeBuffer = elementAlgebra.bufferFactory
fun generateEnvelope(extrema: List<Int>, paddedExtremeValues: Buffer<Double>): Series<Double> { fun generateEnvelope(extrema: List<Int>, paddedExtremeValues: Buffer<T>): Series<T> {
val envelopeFunction = interpolator.interpolate( val envelopeFunction = interpolator.interpolate(
Buffer(extrema.size) { signal.labels[extrema[it]].toDouble() }, makeBuffer(extrema.size) { signal.labels[extrema[it]] },
paddedExtremeValues paddedExtremeValues
) )
return signal.mapWithLabel { _, label -> return signal.mapWithLabel { _, label ->
// For some reason PolynomialInterpolator is exclusive and the right boundary // For some reason PolynomialInterpolator is exclusive and the right boundary
// TODO Notify interpolator authors // TODO Notify interpolator authors
envelopeFunction(label.toDouble()) ?: paddedExtremeValues.last() envelopeFunction(label) ?: paddedExtremeValues.last()
// need to make the interpolator yield values outside boundaries? // need to make the interpolator yield values outside boundaries?
} }
} }
@ -90,13 +90,13 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
* @return [SiftingResult.NotEnoughExtrema] is returned if the signal has too few extrema to extract a mode. * @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. * Success of an appropriate type (See [SiftingResult.Success] class) is returned otherwise.
*/ */
private fun sift(signal: Series<Double>): SiftingResult = siftInner(signal, 1, 0) private fun sift(signal: Series<T>): SiftingResult = siftInner(signal, 1, 0)
/** /**
* Compute a single iteration of the sifting process. * Compute a single iteration of the sifting process.
*/ */
private tailrec fun siftInner( private tailrec fun siftInner(
prevMode: Series<Double>, prevMode: Series<T>,
iterationNumber: Int, iterationNumber: Int,
sNumber: Int sNumber: Int
): SiftingResult = (seriesAlgebra) { ): SiftingResult = (seriesAlgebra) {
@ -108,7 +108,8 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
return when { return when {
iterationNumber >= maxSiftIterations -> SiftingResult.MaxIterationsReached(mode) iterationNumber >= maxSiftIterations -> SiftingResult.MaxIterationsReached(mode)
sNumber >= sConditionThreshold -> SiftingResult.SNumberReached(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) else -> siftInner(mode, iterationNumber + 1, newSNumber)
} }
} }
@ -121,8 +122,8 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
* Modes returned in a list which contains as many modes as it was possible * Modes returned in a list which contains as many modes as it was possible
* to extract before triggering one of the termination conditions. * to extract before triggering one of the termination conditions.
*/ */
public fun decompose(signal: Series<Double>): EMDecompositionResult<Double> = (seriesAlgebra) { public fun decompose(signal: Series<T>): EMDecompositionResult<T> = (seriesAlgebra) {
val modes = mutableListOf<Series<Double>>() val modes = mutableListOf<Series<T>>()
var residual = signal var residual = signal
repeat(nModes) { repeat(nModes) {
val nextMode = when(val r = sift(residual)) { val nextMode = when(val r = sift(residual)) {
@ -130,14 +131,15 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
return EMDecompositionResult( return EMDecompositionResult(
if (it == 0) EMDTerminationReason.SIGNAL_TOO_FLAT if (it == 0) EMDTerminationReason.SIGNAL_TOO_FLAT
else EMDTerminationReason.ALL_POSSIBLE_MODES_EXTRACTED, else EMDTerminationReason.ALL_POSSIBLE_MODES_EXTRACTED,
modes modes,
residual
) )
is SiftingResult.Success<*> -> r.result is SiftingResult.Success<*> -> r.result
} }
modes.add(nextMode as Series<Double>) // TODO remove unchecked cast modes.add(nextMode as Series<T>) // TODO remove unchecked cast
residual = residual.zip(nextMode) { l, r -> l - r } residual = residual.zip(nextMode) { l, r -> l - r }
} }
return EMDecompositionResult(EMDTerminationReason.MAX_MODES_REACHED, modes) return EMDecompositionResult(EMDTerminationReason.MAX_MODES_REACHED, modes, residual)
} }
} }
@ -155,13 +157,13 @@ public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
* @param nModes how many modes should be extracted at most. The algorithm may return fewer modes if it was not * @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. * possible to extract more modes from the signal.
*/ */
public fun <L: Number, A: Field<Double>, BA> SeriesAlgebra<Double, A, BA, L>.empiricalModeDecomposition( public fun <T: Comparable<T>, L: T, A: Field<T>, BA> SeriesAlgebra<T, A, BA, L>.empiricalModeDecomposition(
sConditionThreshold: Int = 15, sConditionThreshold: Int = 15,
maxSiftIterations: Int = 20, maxSiftIterations: Int = 20,
siftingDelta: Double = 1e-2, siftingDelta: T,
nModes: Int = 3 nModes: Int = 3
): EmpiricalModeDecomposition<A, BA, L> ): EmpiricalModeDecomposition<T, A, BA, L>
where BA: BufferAlgebra<Double, A>, BA: FieldOps<Buffer<Double>> = EmpiricalModeDecomposition( where BA: BufferAlgebra<T, A>, BA: FieldOps<Buffer<T>> = EmpiricalModeDecomposition(
seriesAlgebra = this, seriesAlgebra = this,
sConditionThreshold = sConditionThreshold, sConditionThreshold = sConditionThreshold,
maxSiftIterations = maxSiftIterations, maxSiftIterations = maxSiftIterations,
@ -281,5 +283,6 @@ public enum class EMDTerminationReason {
public data class EMDecompositionResult<T>( public data class EMDecompositionResult<T>(
val terminatedBecause: EMDTerminationReason, val terminatedBecause: EMDTerminationReason,
val modes: List<Series<T>> val modes: List<Series<T>>,
val residual: Series<T>
) )