This commit is contained in:
Igor Dunaev 2024-05-04 20:32:23 +03:00
parent 47a9bf0e9a
commit 2a2c5e8765

View File

@ -8,8 +8,6 @@ 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
@ -29,13 +27,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<BA, L: Number> (
private val seriesAlgebra: SeriesAlgebra<Double, *, BA, L>,
public class EmpiricalModeDecomposition<A: Field<Double>, BA, L: Number> (
private val seriesAlgebra: SeriesAlgebra<Double, A, BA, L>,
private val sConditionThreshold: Int = 15,
private val maxSiftIterations: Int = 20,
private val siftingDelta: Double = 1e-2,
private val nModes: Int = 6
) where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>> {
) where BA: BufferAlgebra<Double, A>, BA: FieldOps<Buffer<Double>> {
/**
* Take a signal, construct an upper and a lower envelopes, find the mean value of two,
@ -46,22 +44,22 @@ public class EmpiricalModeDecomposition<BA, L: Number> (
* the signal does not have enough extrema to construct envelopes.
*/
private fun findMean(signal: Series<Double>): Series<Double>? = (seriesAlgebra) {
val interpolator = SplineInterpolator(Float64Field)
fun generateEnvelope(extrema: List<Int>, paddedExtremeValues: DoubleArray): Series<Double> {
val interpolator = SplineInterpolator(elementAlgebra)
fun generateEnvelope(extrema: List<Int>, paddedExtremeValues: Array<Double>): Series<Double> {
val envelopeFunction = interpolator.interpolate(
Buffer(extrema.size) { signal.labels[extrema[it]].toDouble() },
Buffer(extrema.size) { signal.labels[extrema[it]] as Double },
paddedExtremeValues.asBuffer()
)
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 as Double) ?: 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 = Array(maxima.size) { signal[maxima[it]] }
if (maxValues[0] < maxValues[1]) {
maxValues[0] = maxValues[1]
}
@ -70,7 +68,7 @@ public class EmpiricalModeDecomposition<BA, L: Number> (
}
val minima = listOf(0) + signal.troughs() + (signal.size - 1)
val minValues = DoubleArray(minima.size) { signal[minima[it]] }
val minValues = Array(minima.size) { signal[minima[it]] }
if (minValues[0] > minValues[1]) {
minValues[0] = minValues[1]
}
@ -157,13 +155,13 @@ public class EmpiricalModeDecomposition<BA, L: Number> (
* @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 <L: Number, BA> SeriesAlgebra<Double, *, BA, L>.empiricalModeDecomposition(
public fun <L: Number, A: Field<Double>, BA> SeriesAlgebra<Double, A, BA, L>.empiricalModeDecomposition(
sConditionThreshold: Int = 15,
maxSiftIterations: Int = 20,
siftingDelta: Double = 1e-2,
nModes: Int = 3
): EmpiricalModeDecomposition<BA, L>
where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>> = EmpiricalModeDecomposition(
): EmpiricalModeDecomposition<A, BA, L>
where BA: BufferAlgebra<Double, A>, BA: FieldOps<Buffer<Double>> = EmpiricalModeDecomposition(
seriesAlgebra = this,
sConditionThreshold = sConditionThreshold,
maxSiftIterations = maxSiftIterations,
@ -188,13 +186,14 @@ private fun Series<Double>.countZeros(): Int {
/**
* Compute relative difference of two series.
*/
private fun <BA> SeriesAlgebra<Double, *, BA, *>.relativeDifference(
current: Series<Double>,
previous: Series<Double>
):Double where BA: BufferAlgebra<Double, *>, BA: RingOps<Buffer<Double>> =
(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 <T, A: Ring<T>, BA> SeriesAlgebra<T, A, BA, *>.relativeDifference(
current: Series<T>,
previous: Series<T>
):T where BA: BufferAlgebra<T, A>, BA: FieldOps<Buffer<T>> = (bufferAlgebra) {
((current - previous) * (current - previous))
.div(previous * previous)
.fold(elementAlgebra.zero) { acc, it -> acc + it}
}
/**
* Brute force count all extrema of a series.