From 8a07140f7c91e3ff0cf18cdc6ef95b9ac02e59da Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Wed, 19 May 2021 22:35:06 +0300 Subject: [PATCH] Minor piecewise rework --- .../kscience/kmath/functions/interpolate.kt | 2 + .../kscience/kmath/functions/Piecewise.kt | 34 +++++++++---- .../kscience/kmath/functions/Polynomial.kt | 48 +++++++++++++++---- .../kmath/interpolation/LinearInterpolator.kt | 3 +- .../kmath/interpolation/SplineInterpolator.kt | 3 +- .../kmath/functions/PolynomialTest.kt | 17 +++++++ .../kmath/histogram/TreeHistogramSpace.kt | 11 ++--- .../kmath/histogram/UnivariateHistogram.kt | 5 +- 8 files changed, 91 insertions(+), 32 deletions(-) create mode 100644 kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/PolynomialTest.kt diff --git a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt b/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt index 69ba80fb6..b05c003a6 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt @@ -10,12 +10,14 @@ import space.kscience.kmath.interpolation.interpolatePolynomials import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.structures.DoubleBuffer import space.kscience.plotly.Plotly +import space.kscience.plotly.UnstablePlotlyAPI import space.kscience.plotly.makeFile import space.kscience.plotly.models.functionXY import space.kscience.plotly.scatter import kotlin.math.PI import kotlin.math.sin +@OptIn(UnstablePlotlyAPI::class) fun main() { val data = (0..10).map { val x = it.toDouble() / 5 * PI diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt index 6bb2a1656..dd9c571b0 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt @@ -22,8 +22,20 @@ public fun interface Piecewise { /** * Represents piecewise-defined function where all the sub-functions are polynomials. + * @param pieces An ordered list of range-polynomial pairs. The list does not in general guarantee that there are no "holes" in it. */ -public fun interface PiecewisePolynomial : Piecewise> +public class PiecewisePolynomial>( + public val pieces: List, Polynomial>>, +) : Piecewise> { + + public override fun findPiece(arg: T): Polynomial? { + return if (arg < pieces.first().first.start || arg >= pieces.last().first.endInclusive) + null + else { + pieces.firstOrNull { arg in it.first }?.second + } + } +} /** * A [Piecewise] builder where all the pieces are ordered by the [Comparable] type instances. @@ -31,7 +43,7 @@ public fun interface PiecewisePolynomial : Piecewise> * @param T the comparable piece key type. * @param delimiter the initial piecewise separator */ -public class OrderedPiecewisePolynomial>(delimiter: T) : PiecewisePolynomial { +public class PiecewiseBuilder>(delimiter: T) { private val delimiters: MutableList = arrayListOf(delimiter) private val pieces: MutableList> = arrayListOf() @@ -59,17 +71,19 @@ public class OrderedPiecewisePolynomial>(delimiter: T) : Piece pieces.add(0, piece) } - public 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") - } + public fun build(): PiecewisePolynomial { + return PiecewisePolynomial(delimiters.zipWithNext { l, r -> l..r }.zip(pieces)) } } +/** + * A builder for [PiecewisePolynomial] + */ +public fun > PiecewisePolynomial( + startingPoint: T, + builder: PiecewiseBuilder.() -> Unit, +): PiecewisePolynomial = PiecewiseBuilder(startingPoint).apply(builder).build() + /** * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise * definition. diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt index 427a17385..4976befcb 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt @@ -5,10 +5,8 @@ package space.kscience.kmath.functions -import space.kscience.kmath.operations.Group -import space.kscience.kmath.operations.Ring -import space.kscience.kmath.operations.ScaleOperations -import space.kscience.kmath.operations.invoke +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.operations.* import kotlin.contracts.InvocationKind import kotlin.contracts.contract import kotlin.math.max @@ -19,7 +17,7 @@ import kotlin.math.pow * * @param coefficients constant is the leftmost coefficient. */ -public class Polynomial(public val coefficients: List){ +public class Polynomial(public val coefficients: List) { override fun toString(): String = "Polynomial$coefficients" } @@ -32,7 +30,9 @@ public fun Polynomial(vararg coefficients: T): Polynomial = Polynomial(co /** * Evaluates the value of the given double polynomial for given double argument. */ -public fun Polynomial.value(): Double = coefficients.reduceIndexed { index, acc, d -> acc + d.pow(index) } +public fun Polynomial.value(arg: Double): Double = coefficients.reduceIndexed { index, acc, c -> + acc + c * arg.pow(index) +} /** * Evaluates the value of the given polynomial for given argument. @@ -50,9 +50,38 @@ public fun > Polynomial.value(ring: C, arg: T): T = ring { /** * Represent the polynomial as a regular context-less function. */ -public fun > Polynomial.asFunction(ring: C): (T) -> T = { value(ring, it) } +public fun > Polynomial.asFunction(ring: C): (T) -> T = { value(ring, it) } -//public fun +/** + * Create a polynomial witch represents differentiated version of this polynomial + */ +@UnstableKMathAPI +public fun Polynomial.differentiate( + algebra: A, +): Polynomial where A : Ring, A : NumericAlgebra = algebra { + Polynomial(coefficients.drop(1).mapIndexed { index, t -> number(index) * t }) +} + +/** + * Create a polynomial witch represents indefinite integral version of this polynomial + */ +@UnstableKMathAPI +public fun Polynomial.integrate( + algebra: A, +): Polynomial where A : Field, A : NumericAlgebra = algebra { + Polynomial(coefficients.mapIndexed { index, t -> t / number(index) }) +} + +/** + * Compute a definite integral of a given polynomial in a [range] + */ +@UnstableKMathAPI +public fun , A> Polynomial.integrate( + algebra: A, + range: ClosedRange, +): T where A : Field, A : NumericAlgebra = algebra { + value(algebra, range.endInclusive) - value(algebra, range.start) +} /** * Space of polynomials. @@ -87,6 +116,9 @@ public class PolynomialSpace( * Evaluates the polynomial for the given value [arg]. */ public operator fun Polynomial.invoke(arg: T): T = value(ring, arg) + + public fun Polynomial.asFunction(): (T) -> T = asFunction(ring) + } public inline fun C.polynomial(block: PolynomialSpace.() -> R): R where C : Ring, C : ScaleOperations { diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt index 3fbf6157e..24c049647 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt @@ -6,7 +6,6 @@ package space.kscience.kmath.interpolation import space.kscience.kmath.data.XYColumnarData -import space.kscience.kmath.functions.OrderedPiecewisePolynomial import space.kscience.kmath.functions.PiecewisePolynomial import space.kscience.kmath.functions.Polynomial import space.kscience.kmath.misc.UnstableKMathAPI @@ -28,7 +27,7 @@ public class LinearInterpolator>(public override val algebra: require(points.size > 0) { "Point array should not be empty" } insureSorted(points) - OrderedPiecewisePolynomial(points.x[0]).apply { + PiecewisePolynomial(points.x[0]) { for (i in 0 until points.size - 1) { val slope = (points.y[i + 1] - points.y[i]) / (points.x[i + 1] - points.x[i]) val const = points.y[i] - slope * points.x[i] diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt index c3d6ffae7..bf291c315 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt @@ -6,7 +6,6 @@ package space.kscience.kmath.interpolation import space.kscience.kmath.data.XYColumnarData -import space.kscience.kmath.functions.OrderedPiecewisePolynomial import space.kscience.kmath.functions.PiecewisePolynomial import space.kscience.kmath.functions.Polynomial import space.kscience.kmath.misc.UnstableKMathAPI @@ -50,7 +49,7 @@ public class SplineInterpolator>( // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants) - OrderedPiecewisePolynomial(points.x[points.size - 1]).apply { + PiecewisePolynomial(points.x[points.size - 1]) { var cOld = zero for (j in n - 1 downTo 0) { diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/PolynomialTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/PolynomialTest.kt new file mode 100644 index 000000000..05c16d17e --- /dev/null +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/PolynomialTest.kt @@ -0,0 +1,17 @@ +/* + * Copyright 2018-2021 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.functions + +import kotlin.test.Test +import kotlin.test.assertEquals + +class PolynomialTest { + @Test + fun testIntegration() { + val polynomial = Polynomial(1.0, -2.0, 1.0) + assertEquals(0.0, polynomial.value(1.0), 0.001) + } +} \ No newline at end of file diff --git a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramSpace.kt b/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramSpace.kt index 5a217f6c2..6ae8b5ee3 100644 --- a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramSpace.kt +++ b/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramSpace.kt @@ -28,7 +28,6 @@ private fun > TreeMap.getBin(val @UnstableKMathAPI public class TreeHistogram( - override val context: TreeHistogramSpace, private val binMap: TreeMap, ) : UnivariateHistogram { override fun get(value: Double): UnivariateBin? = binMap.getBin(value) @@ -79,15 +78,15 @@ public class TreeHistogramSpace( val count = binCounter.counter.value resBins[key] = UnivariateBin(binCounter.domain, count, sqrt(count)) } - return TreeHistogram(this, resBins) + return TreeHistogram(resBins) } override fun add( a: UnivariateHistogram, b: UnivariateHistogram, ): UnivariateHistogram { - require(a.context == this) { "Histogram $a does not belong to this context" } - require(b.context == this) { "Histogram $b does not belong to this context" } +// require(a.context == this) { "Histogram $a does not belong to this context" } +// require(b.context == this) { "Histogram $b does not belong to this context" } val bins = TreeMap().apply { (a.bins.map { it.domain } union b.bins.map { it.domain }).forEach { def -> put(def.center, @@ -100,7 +99,7 @@ public class TreeHistogramSpace( ) } } - return TreeHistogram(this, bins) + return TreeHistogram(bins) } override fun scale(a: UnivariateHistogram, value: Double): UnivariateHistogram { @@ -116,7 +115,7 @@ public class TreeHistogramSpace( } } - return TreeHistogram(this, bins) + return TreeHistogram(bins) } override fun UnivariateHistogram.unaryMinus(): UnivariateHistogram = this * (-1) diff --git a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt b/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt index 645116ade..70125e22e 100644 --- a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt +++ b/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt @@ -7,8 +7,6 @@ package space.kscience.kmath.histogram import space.kscience.kmath.domains.UnivariateDomain import space.kscience.kmath.misc.UnstableKMathAPI -import space.kscience.kmath.operations.Group -import space.kscience.kmath.operations.GroupElement import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.asSequence @@ -35,8 +33,7 @@ public class UnivariateBin( } @OptIn(UnstableKMathAPI::class) -public interface UnivariateHistogram : Histogram, - GroupElement> { +public interface UnivariateHistogram : Histogram{ public operator fun get(value: Double): UnivariateBin? public override operator fun get(point: Buffer): UnivariateBin? = get(point[0])