From d56b4148bef72f6cffe39e7426e47c5bdbed6240 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sun, 12 Jan 2020 20:51:16 +0300 Subject: [PATCH] Generalized linear interpolation --- .../scientifik.kmath.functions/Piecewise.kt | 55 +++++++++++++++++++ .../scientifik.kmath.functions/Polynomial.kt | 31 ++--------- .../kmath/interpolation/LinearInterpolator.kt | 16 ++++-- .../interpolation/LinearInterpolatorTest.kt | 4 +- 4 files changed, 71 insertions(+), 35 deletions(-) create mode 100644 kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt diff --git a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt new file mode 100644 index 000000000..c05660d52 --- /dev/null +++ b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt @@ -0,0 +1,55 @@ +package scientifik.kmath.functions + +import scientifik.kmath.operations.Ring + +interface Piecewise { + fun findPiece(arg: T): R? +} + +interface PiecewisePolynomial : Piecewise> + +/** + * Ordered list of pieces in piecewise function + */ +class OrderedPiecewisePolynomial>(left: T) : PiecewisePolynomial { + + private val delimiters: ArrayList = arrayListOf(left) + private val pieces: ArrayList> = ArrayList() + + /** + * Dynamically add a piece to the "right" side (beyond maximum argument value of previous piece) + * @param right new rightmost position. If is less then current rightmost position, a error is thrown. + */ + fun putRight(right: T, piece: Polynomial) { + require(right > delimiters.last()) { "New delimiter should be to the right of old one" } + delimiters.add(right) + pieces.add(piece) + } + + fun putLeft(left: T, piece: Polynomial) { + require(left < delimiters.first()) { "New delimiter should be to the left of old one" } + delimiters.add(0, left) + pieces.add(0, piece) + } + + override fun findPiece(arg: T): Polynomial? { + if (arg < delimiters.first() || arg >= delimiters.last()) { + return null + } else { + for (index in 1 until delimiters.size) { + if (arg < delimiters[index]) { + return pieces[index - 1] + } + } + error("Piece not found") + } + } +} + +/** + * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition. + */ +fun , C : Ring> PiecewisePolynomial.value(ring: C, arg: T): T? = + findPiece(arg)?.value(ring, arg) + +fun , C : Ring> PiecewisePolynomial.asFunction(ring: C): (T) -> T? = { value(ring, it) } \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt index a05462863..bd5062661 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt @@ -1,9 +1,7 @@ package scientifik.kmath.functions -import scientifik.kmath.operations.RealField import scientifik.kmath.operations.Ring import scientifik.kmath.operations.Space -import kotlin.jvm.JvmName import kotlin.math.max import kotlin.math.pow @@ -20,11 +18,11 @@ fun Polynomial.value() = fun > Polynomial.value(ring: C, arg: T): T = ring.run { - if( coefficients.isEmpty()) return@run zero + if (coefficients.isEmpty()) return@run zero var res = coefficients.first() var powerArg = arg - for( index in 1 until coefficients.size){ - res += coefficients[index]*powerArg + for (index in 1 until coefficients.size) { + res += coefficients[index] * powerArg //recalculating power on each step to avoid power costs on long polynomials powerArg *= arg } @@ -43,9 +41,6 @@ fun > Polynomial.asMathFunction(): MathFunction> Polynomial.asFunction(ring: C): (T) -> T = { value(ring, it) } -@JvmName("asRealUFunction") -fun Polynomial.asFunction(): (Double) -> Double = asFunction(RealField) - /** * An algebra for polynomials */ @@ -73,22 +68,4 @@ class PolynomialSpace>(val ring: C) : Space> fun , R> C.polynomial(block: PolynomialSpace.() -> R): R { return PolynomialSpace(this).run(block) -} - -class PiecewisePolynomial> internal constructor( - val lowerBoundary: T, - val pieces: List>> -) - -private fun > PiecewisePolynomial.findPiece(arg: T): Polynomial? { - if (arg < lowerBoundary || arg > pieces.last().first) return null - return pieces.first { arg < it.first }.second -} - -/** - * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition. - */ -fun , C : Ring> PiecewisePolynomial.value(ring: C, arg: T): T? = - findPiece(arg)?.value(ring, arg) - -fun , C : Ring> PiecewisePolynomial.asFunction(ring: C): (T) -> T? = { value(ring, it) } \ No newline at end of file +} \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt index 2fc8aaccc..720f1080b 100644 --- a/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt @@ -1,5 +1,6 @@ package scientifik.kmath.interpolation +import scientifik.kmath.functions.OrderedPiecewisePolynomial import scientifik.kmath.functions.PiecewisePolynomial import scientifik.kmath.functions.Polynomial import scientifik.kmath.operations.Field @@ -10,15 +11,18 @@ import scientifik.kmath.operations.Field class LinearInterpolator>(override val algebra: Field) : PolynomialInterpolator { override fun interpolatePolynomials(points: Collection>): PiecewisePolynomial = algebra.run { + require(points.isNotEmpty()) { "Point array should not be empty" } + //sorting points val sorted = points.sortedBy { it.first } - val pairs: List>> = (0 until points.size - 1).map { i -> - val slope = (sorted[i + 1].second - sorted[i].second) / (sorted[i + 1].first - sorted[i].first) - val const = sorted[i].second - slope * sorted[i].first - sorted[i + 1].first to Polynomial(const, slope) + return@run OrderedPiecewisePolynomial(points.first().first).apply { + for (i in 0 until points.size - 1) { + val slope = (sorted[i + 1].second - sorted[i].second) / (sorted[i + 1].first - sorted[i].first) + val const = sorted[i].second - slope * sorted[i].first + val polynomial = Polynomial(const, slope) + putRight(sorted[i + 1].first, polynomial) + } } - - return PiecewisePolynomial(sorted.first().first, pairs) } } \ No newline at end of file diff --git a/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt b/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt index 71303107e..6d35d19f1 100644 --- a/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt +++ b/kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt @@ -18,8 +18,8 @@ class LinearInterpolatorTest { val polynomial = LinearInterpolator(RealField).interpolatePolynomials(data) val function = polynomial.asFunction(RealField) -// assertEquals(null, function(-1.0)) -// assertEquals(0.5, function(0.5)) + assertEquals(null, function(-1.0)) + assertEquals(0.5, function(0.5)) assertEquals(2.0, function(1.5)) assertEquals(3.0, function(2.0)) }