kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

Feature: Polynomials and rational functions #469

Merged
lounres merged 132 commits from feature/polynomials into dev 2022-07-28 18:04:06 +03:00
Conversation 19 Commits 132 Files Changed 51 +580
8 changed files with 580 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit 5928adfe45 - Show all commits

132
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt Normal file
View File

@ -0,0 +1,132 @@
/*
* 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 space.kscience.kmath.misc.PerformancePitfall
import space.kscience.kmath.operations.Ring
/**
* Represents piecewise-defined function.
*
* @param T the piece key type.
* @param R the sub-function type.
*/
public fun interface Piecewise<in T, out R> {
/**
* Returns the appropriate sub-function for given piece key.
*/
public fun findPiece(arg: T): R?
}
/**
* Represents piecewise-defined function where all the sub-functions are polynomials.
*
* @property pieces An ordered list of range-polynomial pairs. The list does not in general guarantee that there are no
* "holes" in it.
*/
public interface PiecewisePolynomial<T : Comparable<T>> : Piecewise<T, ListPolynomial<T>> {
public val pieces: Collection<Pair<ClosedRange<T>, ListPolynomial<T>>>
override fun findPiece(arg: T): ListPolynomial<T>?
}
/**
* A generic piecewise without constraints on how pieces are placed
*/
@PerformancePitfall("findPiece method of resulting piecewise is slow")
public fun <T : Comparable<T>> PiecewisePolynomial(
pieces: Collection<Pair<ClosedRange<T>, ListPolynomial<T>>>,
): PiecewisePolynomial<T> = object : PiecewisePolynomial<T> {
override val pieces: Collection<Pair<ClosedRange<T>, ListPolynomial<T>>> = pieces
override fun findPiece(arg: T): ListPolynomial<T>? = pieces.firstOrNull { arg in it.first }?.second
}
/**
* An optimized piecewise that uses not separate pieces, but a range separated by delimiters.
* The pieces search is logarithmic.
*/
private class OrderedPiecewisePolynomial<T : Comparable<T>>(
override val pieces: List<Pair<ClosedRange<T>, ListPolynomial<T>>>,
) : PiecewisePolynomial<T> {
override fun findPiece(arg: T): ListPolynomial<T>? {
val index = pieces.binarySearch { (range, _) ->
when {
arg >= range.endInclusive -> -1
arg < range.start -> +1
else -> 0
}
}
return if (index < 0) null else pieces[index].second
}
}
/**
* A [Piecewise] builder where all the pieces are ordered by the [Comparable] type instances.
*
* @param T the comparable piece key type.
* @param delimiter the initial piecewise separator
*/
public class PiecewiseBuilder<T : Comparable<T>>(delimiter: T) {
private val delimiters: MutableList<T> = arrayListOf(delimiter)
private val pieces: MutableList<ListPolynomial<T>> = arrayListOf()
/**
* Dynamically adds a piece to the right side (beyond maximum argument value of previous piece)
*
* @param right new rightmost position. If is less than current rightmost position, an error is thrown.
* @param piece the sub-function.
*/
public fun putRight(right: T, piece: ListPolynomial<T>) {
require(right > delimiters.last()) { "New delimiter should be to the right of old one" }
delimiters += right
pieces += piece
}
/**
* Dynamically adds a piece to the left side (beyond maximum argument value of previous piece)
*
* @param left the new leftmost position. If is less than current rightmost position, an error is thrown.
* @param piece the sub-function.
*/
public fun putLeft(left: T, piece: ListPolynomial<T>) {
require(left < delimiters.first()) { "New delimiter should be to the left of old one" }
delimiters.add(0, left)
pieces.add(0, piece)
}
public fun build(): PiecewisePolynomial<T> = OrderedPiecewisePolynomial(delimiters.zipWithNext { l, r ->
l..r
}.zip(pieces))
}
/**
* A builder for [PiecewisePolynomial]
*/
public fun <T : Comparable<T>> PiecewisePolynomial(
startingPoint: T,
builder: PiecewiseBuilder<T>.() -> Unit,
): PiecewisePolynomial<T> = PiecewiseBuilder(startingPoint).apply(builder).build()
/**
* Return a value of polynomial function with given [ring] a given [arg] or null if argument is outside piecewise
* definition.
*/
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.substitute(ring: C, arg: T): T? =
findPiece(arg)?.substitute(ring, arg)
/**
* Convert this polynomial to a function returning nullable value (null if argument is outside piecewise range).
*/
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { substitute(ring, it) }
/**
* Convert this polynomial to a function using [defaultValue] for arguments outside the piecewise range.
*/
public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C, defaultValue: T): (T) -> T =
{ substitute(ring, it) ?: defaultValue }

112
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt Normal file
View File

@ -0,0 +1,112 @@
/*
* 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.
*/
/*
* 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.integration
import space.kscience.kmath.functions.PiecewisePolynomial
import space.kscience.kmath.functions.integrate
import space.kscience.kmath.functions.antiderivative
import space.kscience.kmath.interpolation.PolynomialInterpolator
import space.kscience.kmath.interpolation.SplineInterpolator
import space.kscience.kmath.interpolation.interpolatePolynomials
import space.kscience.kmath.misc.PerformancePitfall
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.*
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.DoubleBuffer
import space.kscience.kmath.structures.MutableBufferFactory
/**
* Compute analytical indefinite integral of this [PiecewisePolynomial], keeping all intervals intact
*/
@OptIn(PerformancePitfall::class)
@UnstableKMathAPI
public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(algebra: Field<T>): PiecewisePolynomial<T> =
PiecewisePolynomial(pieces.map { it.first to it.second.antiderivative(algebra) })
/**
* Compute definite integral of given [PiecewisePolynomial] piece by piece in a given [range]
* Requires [UnivariateIntegrationNodes] or [IntegrationRange] and [IntegrandMaxCalls]
*
* TODO use context receiver for algebra
*/
@UnstableKMathAPI
public fun <T : Comparable<T>> PiecewisePolynomial<T>.integrate(
algebra: Field<T>, range: ClosedRange<T>,
): T = algebra.sum(
pieces.map { (region, poly) ->
val intersectedRange = maxOf(range.start, region.start)..minOf(range.endInclusive, region.endInclusive)
//Check if polynomial range is not used
if (intersectedRange.start == intersectedRange.endInclusive) algebra.zero
else poly.integrate(algebra, intersectedRange)
}
)
/**
* A generic spline-interpolation-based analytic integration
* * [IntegrationRange]&mdash;the univariate range of integration. By default, uses `0..1` interval.
* * [IntegrandMaxCalls]&mdash;the maximum number of function calls during integration. For non-iterative rules, always uses
* the maximum number of points. By default, uses 10 points.
*/
@UnstableKMathAPI
public class SplineIntegrator<T : Comparable<T>>(
public val algebra: Field<T>,
public val bufferFactory: MutableBufferFactory<T>,
) : UnivariateIntegrator<T> {
override fun process(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = algebra {
val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
val interpolator: PolynomialInterpolator<T> = SplineInterpolator(algebra, bufferFactory)
val nodes: Buffer<Double> = integrand.getFeature<UnivariateIntegrationNodes>()?.nodes ?: run {
val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
val step = (range.endInclusive - range.start) / (numPoints - 1)
DoubleBuffer(numPoints) { i -> range.start + i * step }
}
val values = nodes.map(bufferFactory) { integrand.function(it) }
val polynomials = interpolator.interpolatePolynomials(
nodes.map(bufferFactory) { number(it) },
values
)
val res = polynomials.integrate(algebra, number(range.start)..number(range.endInclusive))
integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + nodes.size)
}
}
/**
* A simplified double-based spline-interpolation-based analytic integration
* * [IntegrationRange]&mdash;the univariate range of integration. By default, uses `0.0..1.0` interval.
* * [IntegrandMaxCalls]&mdash;the maximum number of function calls during integration. For non-iterative rules, always
* uses the maximum number of points. By default, uses 10 points.
*/
@UnstableKMathAPI
public object DoubleSplineIntegrator : UnivariateIntegrator<Double> {
override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
val interpolator: PolynomialInterpolator<Double> = SplineInterpolator(DoubleField, ::DoubleBuffer)
val nodes: Buffer<Double> = integrand.getFeature<UnivariateIntegrationNodes>()?.nodes ?: run {
val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
val step = (range.endInclusive - range.start) / (numPoints - 1)
DoubleBuffer(numPoints) { i -> range.start + i * step }
}
val values = nodes.map { integrand.function(it) }
val polynomials = interpolator.interpolatePolynomials(nodes, values)
val res = polynomials.integrate(DoubleField, range)
return integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + nodes.size)
}
}
@Suppress("unused")
@UnstableKMathAPI
public inline val DoubleField.splineIntegrator: UnivariateIntegrator<Double>
get() = DoubleSplineIntegrator

92
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt Normal file
View File

@ -0,0 +1,92 @@
/*
* 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.
*/
@file:OptIn(UnstableKMathAPI::class)
package space.kscience.kmath.interpolation
import space.kscience.kmath.data.XYColumnarData
import space.kscience.kmath.functions.PiecewisePolynomial
import space.kscience.kmath.functions.asFunction
import space.kscience.kmath.functions.substitute
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.Ring
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.asBuffer
/**
* And interpolator for data with x column type [X], y column type [Y].
*/
public fun interface Interpolator<T, in X : T, Y : T> {
public fun interpolate(points: XYColumnarData<T, X, Y>): (X) -> Y
}
/**
* And interpolator returning [PiecewisePolynomial] function
*/
public interface PolynomialInterpolator<T : Comparable<T>> : Interpolator<T, T, T> {
public val algebra: Ring<T>
public fun getDefaultValue(): T = error("Out of bounds")
public fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T>
override fun interpolate(points: XYColumnarData<T, T, T>): (T) -> T = { x ->
interpolatePolynomials(points).substitute(algebra, x) ?: getDefaultValue()
}
}
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
x: Buffer<T>,
y: Buffer<T>,
): PiecewisePolynomial<T> {
val pointSet = XYColumnarData.of(x, y)
return interpolatePolynomials(pointSet)
}
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
data: Map<T, T>,
): PiecewisePolynomial<T> {
val pointSet = XYColumnarData.of(data.keys.toList().asBuffer(), data.values.toList().asBuffer())
return interpolatePolynomials(pointSet)
}
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
data: List<Pair<T, T>>,
): PiecewisePolynomial<T> {
val pointSet = XYColumnarData.of(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer())
return interpolatePolynomials(pointSet)
}
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
x: Buffer<T>,
y: Buffer<T>,
): (T) -> T? = interpolatePolynomials(x, y).asFunction(algebra)
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
data: Map<T, T>,
): (T) -> T? = interpolatePolynomials(data).asFunction(algebra)
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
data: List<Pair<T, T>>,
): (T) -> T? = interpolatePolynomials(data).asFunction(algebra)
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
x: Buffer<T>,
y: Buffer<T>,
defaultValue: T,
): (T) -> T = interpolatePolynomials(x, y).asFunction(algebra, defaultValue)
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
data: Map<T, T>,
defaultValue: T,
): (T) -> T = interpolatePolynomials(data).asFunction(algebra, defaultValue)
public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolate(
data: List<Pair<T, T>>,
defaultValue: T,
): (T) -> T = interpolatePolynomials(data).asFunction(algebra, defaultValue)

43
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt Normal file
View File

@ -0,0 +1,43 @@
/*
* 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.interpolation
import space.kscience.kmath.data.XYColumnarData
import space.kscience.kmath.functions.PiecewisePolynomial
import space.kscience.kmath.functions.ListPolynomial
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.Field
import space.kscience.kmath.operations.invoke
@OptIn(UnstableKMathAPI::class)
internal fun <T : Comparable<T>> insureSorted(points: XYColumnarData<*, T, *>) {
for (i in 0 until points.size - 1)
require(points.x[i + 1] > points.x[i]) { "Input data is not sorted at index $i" }
}
/**
* Reference JVM implementation: https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/analysis/interpolation/LinearInterpolator.java
*/
public class LinearInterpolator<T : Comparable<T>>(override val algebra: Field<T>) : PolynomialInterpolator<T> {
@OptIn(UnstableKMathAPI::class)
override fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T> = algebra {
require(points.size > 0) { "Point array should not be empty" }
insureSorted(points)
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]
val polynomial = ListPolynomial(const, slope)
putRight(points.x[i + 1], polynomial)
}
}
}
}
public val <T : Comparable<T>> Field<T>.linearInterpolator: LinearInterpolator<T>
get() = LinearInterpolator(this)

88
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt Normal file
View File

@ -0,0 +1,88 @@
/*
* 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.
*/
/*
* 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.interpolation
import space.kscience.kmath.data.XYColumnarData
import space.kscience.kmath.functions.PiecewisePolynomial
import space.kscience.kmath.functions.ListPolynomial
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.DoubleField
import space.kscience.kmath.operations.Field
import space.kscience.kmath.operations.invoke
import space.kscience.kmath.structures.DoubleBuffer
import space.kscience.kmath.structures.MutableBufferFactory
/**
* Generic spline interpolator. Not recommended for performance critical places, use platform-specific and type
* specific ones.
*
* Based on
* https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/SplineInterpolator.java
*/
public class SplineInterpolator<T : Comparable<T>>(
override val algebra: Field<T>,
public val bufferFactory: MutableBufferFactory<T>,
) : PolynomialInterpolator<T> {
//TODO possibly optimize zeroed buffers
@OptIn(UnstableKMathAPI::class)
override fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T> = algebra {
require(points.size >= 3) { "Can't use spline interpolator with less than 3 points" }
insureSorted(points)
// Number of intervals. The number of data points is n + 1.
val n = points.size - 1
// Differences between knot points
val h = bufferFactory(n) { i -> points.x[i + 1] - points.x[i] }
val mu = bufferFactory(n) { zero }
val z = bufferFactory(n + 1) { zero }
for (i in 1 until n) {
val g = 2.0 * (points.x[i + 1] - points.x[i - 1]) - h[i - 1] * mu[i - 1]
mu[i] = h[i] / g
z[i] =
((points.y[i + 1] * h[i - 1] - points.y[i] * (points.x[i + 1] - points.x[i - 1]) + points.y[i - 1] * h[i]) * 3.0 /
(h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g
}
// cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)
PiecewisePolynomial(points.x[points.size - 1]) {
var cOld = zero
for (j in n - 1 downTo 0) {
val c = z[j] - mu[j] * cOld
val a = points.y[j]
val b = (points.y[j + 1] - points.y[j]) / h[j] - h[j] * (cOld + 2.0 * c) / 3.0
val d = (cOld - c) / (3.0 * h[j])
val x0 = points.x[j]
val x02 = x0 * x0
val x03 = x02 * x0
//Shift coefficients to represent absolute polynomial instead of one with an offset
val polynomial = ListPolynomial(
a - b * x0 + c * x02 - d * x03,
b - 2 * c * x0 + 3 * d * x02,
c - 3 * d * x0,
d
)
cOld = c
putLeft(x0, polynomial)
}
}
}
}
public fun <T : Comparable<T>> Field<T>.splineInterpolator(
bufferFactory: MutableBufferFactory<T>,
): SplineInterpolator<T> = SplineInterpolator(this, bufferFactory)
public val DoubleField.splineInterpolator: SplineInterpolator<Double>
get() = SplineInterpolator(this, ::DoubleBuffer)

48
kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SplineIntegralTest.kt Normal file
View File

@ -0,0 +1,48 @@
/*
* 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.integration
import space.kscience.kmath.functions.ListPolynomial
import space.kscience.kmath.functions.integrate
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.DoubleField
import kotlin.math.PI
import kotlin.math.sin
import kotlin.test.Test
import kotlin.test.assertEquals
@OptIn(UnstableKMathAPI::class)
class SplineIntegralTest {
@Test
fun integratePolynomial(){
val polynomial = ListPolynomial(1.0, 2.0, 3.0)
val integral = polynomial.integrate(DoubleField,1.0..2.0)
assertEquals(11.0, integral, 0.001)
}
@Test
fun gaussSin() {
val res = DoubleField.splineIntegrator.integrate(0.0..2 * PI, IntegrandMaxCalls(5)) { x ->
sin(x)
}
assertEquals(0.0, res.value, 1e-2)
}
@Test
fun gaussUniform() {
val res = DoubleField.splineIntegrator.integrate(35.0..100.0, IntegrandMaxCalls(20)) { x ->
if(x in 30.0..50.0){
1.0
} else {
0.0
}
}
assertEquals(15.0, res.value, 0.5)
}
}

29
kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/LinearInterpolatorTest.kt Normal file
View File

@ -0,0 +1,29 @@
/*
* 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.interpolation
import space.kscience.kmath.operations.DoubleField
import kotlin.test.Test
import kotlin.test.assertEquals
internal class LinearInterpolatorTest {
@Test
fun testInterpolation() {
val data = listOf(
0.0 to 0.0,
1.0 to 1.0,
2.0 to 3.0,
3.0 to 4.0
)
//val polynomial: PiecewisePolynomial<Double> = DoubleField.linearInterpolator.interpolatePolynomials(data)
val function = DoubleField.linearInterpolator.interpolate(data)
assertEquals(null, function(-1.0))
assertEquals(0.5, function(0.5))
assertEquals(2.0, function(1.5))
assertEquals(3.0, function(2.0))
}
}

36
kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/SplineInterpolatorTest.kt Normal file
View File

@ -0,0 +1,36 @@
/*
* 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.
*/
/*
* 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.interpolation
import space.kscience.kmath.operations.DoubleField
import kotlin.math.PI
import kotlin.math.sin
import kotlin.test.Test
import kotlin.test.assertEquals
internal class SplineInterpolatorTest {
@Test
fun testInterpolation() {
val data = (0..10).map {
val x = it.toDouble() / 5 * PI
x to sin(x)
}
//val polynomial: PiecewisePolynomial<Double> = DoubleField.splineInterpolator.interpolatePolynomials(data)
val function = DoubleField.splineInterpolator.interpolate(data, Double.NaN)
assertEquals(Double.NaN, function(-1.0))
assertEquals(sin(0.5), function(0.5), 0.1)
assertEquals(sin(1.5), function(1.5), 0.1)
assertEquals(sin(2.0), function(2.0), 0.1)
}
}