From a24c8dcbce3256bb0f707fcf373cbaccbec1760a Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sat, 22 May 2021 20:08:49 +0300 Subject: [PATCH] Simpson and spline integration --- .../kmath/functions/interpolateSquare.kt | 45 +++++++++ .../kmath/operations/algebraExtensions.kt | 2 + .../kmath/structures/NumberNDFieldTest.kt | 4 +- .../kscience/kmath/functions/Polynomial.kt | 15 ++- .../kmath/integration/GaussIntegrator.kt | 8 +- .../kmath/integration/SimpsonIntegrator.kt | 45 ++++++--- .../kmath/integration/SplineIntegrator.kt | 98 ++++++++++++++++++- .../kmath/integration/UnivariateIntegrand.kt | 15 +++ .../kmath/interpolation/Interpolator.kt | 6 ++ .../kmath/integration/SimpsonIntegralTest.kt | 6 +- .../kmath/integration/SplineIntegralTest.kt | 48 +++++++++ settings.gradle.kts | 8 +- 12 files changed, 260 insertions(+), 40 deletions(-) create mode 100644 examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt create mode 100644 kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SplineIntegralTest.kt diff --git a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt b/examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt new file mode 100644 index 000000000..33973c880 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt @@ -0,0 +1,45 @@ +/* + * 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.interpolation.SplineInterpolator +import space.kscience.kmath.interpolation.interpolatePolynomials +import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.real.step +import space.kscience.kmath.structures.map +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 + +@OptIn(UnstablePlotlyAPI::class) +fun main() { + val function: UnivariateFunction = { x -> + if (x in 30.0..50.0) { + 1.0 + } else { + 0.0 + } + } + val xs = 0.0..100.0 step 0.5 + val ys = xs.map(function) + + val polynomial: PiecewisePolynomial = SplineInterpolator.double.interpolatePolynomials(xs, ys) + + val polyFunction = polynomial.asFunction(DoubleField, 0.0) + + Plotly.plot { + scatter { + name = "interpolated" + functionXY(25.0..55.0, 0.1) { polyFunction(it) } + } + scatter { + name = "original" + functionXY(25.0..55.0, 0.1) { function(it) } + } + }.makeFile() +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/algebraExtensions.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/algebraExtensions.kt index 7967aeadb..d52be943a 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/algebraExtensions.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/algebraExtensions.kt @@ -16,6 +16,8 @@ public fun Ring.sum(data: Iterable): T = data.fold(zero) { left, right add(left, right) } +//TODO replace by sumOf with multi-receivers + /** * Returns the sum of all elements in the sequence in this [Ring]. * diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt index 57393fadd..fb51553f7 100644 --- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt @@ -74,9 +74,9 @@ class NumberNDFieldTest { val division = array1.combine(array2, Double::div) } - object L2Norm : Norm, Double> { + object L2Norm : Norm, Double> { @OptIn(PerformancePitfall::class) - override fun norm(arg: StructureND): Double = + override fun norm(arg: StructureND): Double = kotlin.math.sqrt(arg.elements().sumOf { it.second.toDouble() }) } 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 4976befcb..ba77d7b25 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 @@ -69,18 +69,23 @@ public fun Polynomial.differentiate( public fun Polynomial.integrate( algebra: A, ): Polynomial where A : Field, A : NumericAlgebra = algebra { - Polynomial(coefficients.mapIndexed { index, t -> t / number(index) }) + val integratedCoefficients = buildList(coefficients.size + 1) { + add(zero) + coefficients.forEachIndexed{ index, t -> add(t / (number(index) + one)) } + } + Polynomial(integratedCoefficients) } /** * Compute a definite integral of a given polynomial in a [range] */ @UnstableKMathAPI -public fun , A> Polynomial.integrate( - algebra: A, +public fun > Polynomial.integrate( + algebra: Field, range: ClosedRange, -): T where A : Field, A : NumericAlgebra = algebra { - value(algebra, range.endInclusive) - value(algebra, range.start) +): T = algebra { + val integral = integrate(algebra) + integral.value(algebra, range.endInclusive) - integral.value(algebra, range.start) } /** diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt index 26757ae68..283f97557 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt @@ -10,13 +10,7 @@ import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.asBuffer import space.kscience.kmath.structures.indices -/** - * Set of univariate integration ranges. First components correspond to ranges themselves, second components to number of - * integration nodes per range - */ -public class UnivariateIntegrandRanges(public val ranges: List, Int>>) : IntegrandFeature { - public constructor(vararg pairs: Pair, Int>) : this(pairs.toList()) -} + /** * A simple one-pass integrator based on Gauss rule diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SimpsonIntegrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SimpsonIntegrator.kt index 7eca90f88..56647bace 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SimpsonIntegrator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SimpsonIntegrator.kt @@ -6,36 +6,53 @@ package space.kscience.kmath.integration import space.kscience.kmath.operations.Field -import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.BufferFactory +import space.kscience.kmath.operations.invoke +import space.kscience.kmath.operations.sum /** + * Use double pass Simpson rule integration with a fixed number of points. + * Requires [UnivariateIntegrandRanges] or [IntegrationRange] and [IntegrandMaxCalls] * [IntegrationRange] - the univariate range of integration. By default uses 0..1 interval. * [IntegrandMaxCalls] - the maximum number of function calls during integration. For non-iterative rules, always uses the maximum number of points. By default uses 10 points. */ public class SimpsonIntegrator( public val algebra: Field, - public val bufferFactory: BufferFactory, ) : UnivariateIntegrator { - override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand = with(algebra) { - val numPoints = integrand.getFeature()?.maxCalls ?: 100 - require(numPoints >= 4) - val range = integrand.getFeature()?.range ?: 0.0..1.0 + + private fun integrateRange( + integrand: UnivariateIntegrand, range: ClosedRange, numPoints: Int, + ): T = algebra { val h: Double = (range.endInclusive - range.start) / (numPoints - 1) - val points: Buffer = bufferFactory(numPoints) { i -> + val values: List = List(numPoints) { i -> integrand.function(range.start + i * h) }// equally distributed point - fun simpson(index: Int) = h / 3 * (points[index - 1] + 4 * points[index] + points[index + 1]) + //TODO don't use list, reassign values instead + fun simpson(index: Int) = h / 3 * (values[index - 1] + 4 * values[index] + values[index + 1]) + var res = zero res += simpson(1) / 1.5 //border points with 1.5 factor - for (i in 2 until (points.size - 2)) { + for (i in 2 until (values.size - 2)) { + //each half-interval is computed twice, therefore /2 res += simpson(i) / 2 } - res += simpson(points.size - 2) / 1.5 //border points with 1.5 factor - return integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + points.size) + res += simpson(values.size - 2) / 1.5 //border points with 1.5 factor + return res + } + + override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand { + val ranges = integrand.getFeature() + return if (ranges != null) { + val res = algebra.sum(ranges.ranges.map { integrateRange(integrand, it.first, it.second) }) + integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second }) + } else { + val numPoints = integrand.getFeature()?.maxCalls ?: 100 + require(numPoints >= 4) { "Simpson integrator requires at least 4 nodes" } + val range = integrand.getFeature()?.range ?: 0.0..1.0 + val res = integrateRange(integrand, range, numPoints) + integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + numPoints) + } } } -public inline val Field.simpsonIntegrator: SimpsonIntegrator - get() = SimpsonIntegrator(this, Buffer.Companion::auto) \ No newline at end of file +public val Field.simpsonIntegrator: SimpsonIntegrator get() = SimpsonIntegrator(this) \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt index da5730497..23d7bdd8d 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt @@ -5,8 +5,98 @@ package space.kscience.kmath.integration -public class SplineIntegrator: UnivariateIntegrator { - override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand { - TODO("Not yet implemented") +import space.kscience.kmath.functions.PiecewisePolynomial +import space.kscience.kmath.functions.integrate +import space.kscience.kmath.interpolation.PolynomialInterpolator +import space.kscience.kmath.interpolation.SplineInterpolator +import space.kscience.kmath.interpolation.interpolatePolynomials +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.operations.sum +import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.DoubleBuffer +import space.kscience.kmath.structures.MutableBufferFactory +import space.kscience.kmath.structures.map + +/** + * Compute analytical indefinite integral of this [PiecewisePolynomial], keeping all intervals intact + */ +@UnstableKMathAPI +public fun > PiecewisePolynomial.integrate(algebra: Field): PiecewisePolynomial = + PiecewisePolynomial(pieces.map { it.first to it.second.integrate(algebra) }) + +/** + * Compute definite integral of given [PiecewisePolynomial] piece by piece in a given [range] + * Requires [UnivariateIntegrationNodes] or [IntegrationRange] and [IntegrandMaxCalls] + */ +@UnstableKMathAPI +public fun > PiecewisePolynomial.integrate( + algebra: Field, range: ClosedRange, +): 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) } -} \ No newline at end of file +) + +/** + * A generic spline-interpolation-based analytic integration + * [IntegrationRange] - the univariate range of integration. By default uses 0..1 interval. + * [IntegrandMaxCalls] - 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>( + public val algebra: Field, + public val bufferFactory: MutableBufferFactory, +) : UnivariateIntegrator { + override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand = algebra { + val range = integrand.getFeature()?.range ?: 0.0..1.0 + + val interpolator: PolynomialInterpolator = SplineInterpolator(algebra, bufferFactory) + val nodes: Buffer = integrand.getFeature()?.nodes ?: run { + val numPoints = integrand.getFeature()?.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] - the univariate range of integration. By default uses 0..1 interval. + * [IntegrandMaxCalls] - 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 { + override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand { + val range = integrand.getFeature()?.range ?: 0.0..1.0 + + val interpolator: PolynomialInterpolator = SplineInterpolator(DoubleField, ::DoubleBuffer) + val nodes: Buffer = integrand.getFeature()?.nodes ?: run { + val numPoints = integrand.getFeature()?.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) + } +} + +@UnstableKMathAPI +public inline val DoubleField.splineIntegrator: UnivariateIntegrator + get() = DoubleSplineIntegrator \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt index d3d0220f6..8c515065f 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt @@ -6,6 +6,8 @@ package space.kscience.kmath.integration import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.DoubleBuffer import kotlin.jvm.JvmInline import kotlin.reflect.KClass @@ -35,6 +37,19 @@ public typealias UnivariateIntegrator = Integrator> @JvmInline public value class IntegrationRange(public val range: ClosedRange) : IntegrandFeature +/** + * Set of univariate integration ranges. First components correspond to ranges themselves, second components to number of + * integration nodes per range + */ +public class UnivariateIntegrandRanges(public val ranges: List, Int>>) : IntegrandFeature { + public constructor(vararg pairs: Pair, Int>) : this(pairs.toList()) +} + +public class UnivariateIntegrationNodes(public val nodes: Buffer) : IntegrandFeature { + public constructor(vararg nodes: Double) : this(DoubleBuffer(nodes)) +} + + /** * Value of the integrand if it is present or null */ diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt index 08090a522..c9ec0d527 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt @@ -15,10 +15,16 @@ 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 { public fun interpolate(points: XYColumnarData): (X) -> Y } +/** + * And interpolator returning [PiecewisePolynomial] function + */ public interface PolynomialInterpolator> : Interpolator { public val algebra: Ring diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SimpsonIntegralTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SimpsonIntegralTest.kt index b0c0cd098..9f2d71554 100644 --- a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SimpsonIntegralTest.kt +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SimpsonIntegralTest.kt @@ -24,8 +24,8 @@ class SimpsonIntegralTest { @Test fun gaussUniform() { - val res = DoubleField.simpsonIntegrator.integrate(35.0..100.0) { x -> - if(x in 30.0..50.0){ + val res = DoubleField.simpsonIntegrator.integrate(35.0..100.0, IntegrandMaxCalls(20)) { x -> + if (x in 30.0..50.0) { 1.0 } else { 0.0 @@ -33,6 +33,4 @@ class SimpsonIntegralTest { } assertEquals(15.0, res.value, 0.5) } - - } \ No newline at end of file diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SplineIntegralTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SplineIntegralTest.kt new file mode 100644 index 000000000..afeba0be4 --- /dev/null +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SplineIntegralTest.kt @@ -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.Polynomial +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 = Polynomial(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) + } + + +} \ No newline at end of file diff --git a/settings.gradle.kts b/settings.gradle.kts index cc9eb4860..065dd3ac4 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -1,11 +1,11 @@ pluginManagement { repositories { + maven("https://repo.kotlin.link") mavenCentral() gradlePluginPortal() - maven("https://repo.kotlin.link") } - val toolsVersion = "0.9.8" + val toolsVersion = "0.9.9" val kotlinVersion = "1.5.0" plugins { @@ -15,8 +15,8 @@ pluginManagement { kotlin("multiplatform") version kotlinVersion kotlin("jvm") version kotlinVersion kotlin("plugin.allopen") version kotlinVersion - id("org.jetbrains.kotlinx.benchmark") version "0.3.0" - kotlin("jupyter.api") version "0.9.1-61" + id("org.jetbrains.kotlinx.benchmark") version "0.3.1" + kotlin("jupyter.api") version "0.10.0-25" } }