From 65a8d8f581f8b9ec6755530f8ff4a00d0124b9df Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 16 Apr 2021 15:50:33 +0300 Subject: [PATCH] Field element integration --- CHANGELOG.md | 1 + .../kmath/benchmarks/BigIntBenchmark.kt | 14 +-- .../kscience/kmath/functions/integrate.kt | 16 +++ .../kmath/functions/matrixIntegration.kt | 22 ++++ .../kmath/commons/integration/CMIntegrator.kt | 2 +- .../integration/GaussRuleIntegrator.kt | 2 +- kmath-core/api/kmath-core.api | 2 + .../space/kscience/kmath/structures/Buffer.kt | 9 +- kmath-functions/build.gradle.kts | 27 +++-- .../kmath/integration/GaussIntegrator.kt | 105 ++++++++++++------ .../integration/GaussIntegratorRuleFactory.kt | 29 ++++- .../kscience/kmath/integration/Integrand.kt | 4 +- .../kmath/integration/GaussIntegralTest.kt | 4 +- 13 files changed, 172 insertions(+), 65 deletions(-) create mode 100644 examples/src/main/kotlin/space/kscience/kmath/functions/integrate.kt create mode 100644 examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt diff --git a/CHANGELOG.md b/CHANGELOG.md index f7acd08cd..c3bd2641a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ - bindSymbolOrNull - Blocking chains and Statistics - Multiplatform integration +- Integration for any Field element ### Changed - Exponential operations merged with hyperbolic functions diff --git a/examples/src/benchmarks/kotlin/space/kscience/kmath/benchmarks/BigIntBenchmark.kt b/examples/src/benchmarks/kotlin/space/kscience/kmath/benchmarks/BigIntBenchmark.kt index fec80d31a..672efd5c2 100644 --- a/examples/src/benchmarks/kotlin/space/kscience/kmath/benchmarks/BigIntBenchmark.kt +++ b/examples/src/benchmarks/kotlin/space/kscience/kmath/benchmarks/BigIntBenchmark.kt @@ -20,22 +20,22 @@ internal class BigIntBenchmark { val jvmNumber = JBigIntegerField.number(Int.MAX_VALUE) @Benchmark - fun kmAdd(blackhole: Blackhole) = BigIntField{ + fun kmAdd(blackhole: Blackhole) = BigIntField { blackhole.consume(kmNumber + kmNumber + kmNumber) } @Benchmark - fun jvmAdd(blackhole: Blackhole) = JBigIntegerField{ - blackhole.consume(jvmNumber + jvmNumber+ jvmNumber) + fun jvmAdd(blackhole: Blackhole) = JBigIntegerField { + blackhole.consume(jvmNumber + jvmNumber + jvmNumber) } @Benchmark - fun kmMultiply(blackhole: Blackhole) = BigIntField{ - blackhole.consume(kmNumber*kmNumber*kmNumber) + fun kmMultiply(blackhole: Blackhole) = BigIntField { + blackhole.consume(kmNumber * kmNumber * kmNumber) } @Benchmark - fun jvmMultiply(blackhole: Blackhole) = JBigIntegerField{ - blackhole.consume(jvmNumber*jvmNumber*jvmNumber) + fun jvmMultiply(blackhole: Blackhole) = JBigIntegerField { + blackhole.consume(jvmNumber * jvmNumber * jvmNumber) } } \ No newline at end of file diff --git a/examples/src/main/kotlin/space/kscience/kmath/functions/integrate.kt b/examples/src/main/kotlin/space/kscience/kmath/functions/integrate.kt new file mode 100644 index 000000000..761d006d3 --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/functions/integrate.kt @@ -0,0 +1,16 @@ +package space.kscience.kmath.functions + +import space.kscience.kmath.integration.GaussIntegrator +import space.kscience.kmath.integration.value +import kotlin.math.pow + +fun main() { + //Define a function + val function: UnivariateFunction = { x -> 3 * x.pow(2) + 2 * x + 1 } + + //get the result of the integration + val result = GaussIntegrator.legendre(0.0..10.0, function = function) + + //the value is nullable because in some cases the integration could not succeed + println(result.value) +} \ No newline at end of file diff --git a/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt b/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt new file mode 100644 index 000000000..5e92ce22a --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt @@ -0,0 +1,22 @@ +package space.kscience.kmath.functions + +import space.kscience.kmath.integration.GaussIntegrator +import space.kscience.kmath.integration.UnivariateIntegrand +import space.kscience.kmath.integration.value +import space.kscience.kmath.nd.StructureND +import space.kscience.kmath.nd.nd +import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.operations.invoke + +fun main(): Unit = DoubleField { + nd(2, 2) { + //Define a function in a nd space + val function: UnivariateFunction> = { x -> 3 * x.pow(2) + 2 * x + 1 } + + //get the result of the integration + val result: UnivariateIntegrand> = GaussIntegrator.legendre(this, 0.0..10.0, function = function) + + //the value is nullable because in some cases the integration could not succeed + println(result.value) + } +} \ No newline at end of file diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt index 15a734531..535c6b39e 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt @@ -36,7 +36,7 @@ public class CMIntegrator( IntegrandValue(res) + IntegrandAbsoluteAccuracy(integrator.absoluteAccuracy) + IntegrandRelativeAccuracy(integrator.relativeAccuracy) + - IntegrandCalls(integrator.evaluations + integrand.calls) + IntegrandCallsPerformed(integrator.evaluations + integrand.calls) } diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt index 67d032e8b..071bac315 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt @@ -22,7 +22,7 @@ public class GaussRuleIntegrator( val integrator: GaussIntegrator = getIntegrator(range) //TODO check performance val res: Double = integrator.integrate(integrand.function) - return integrand + IntegrandValue(res) + IntegrandCalls(integrand.calls + numpoints) + return integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + numpoints) } private fun getIntegrator(range: ClosedRange): GaussIntegrator { diff --git a/kmath-core/api/kmath-core.api b/kmath-core/api/kmath-core.api index c62abbf12..6b300123c 100644 --- a/kmath-core/api/kmath-core.api +++ b/kmath-core/api/kmath-core.api @@ -1791,6 +1791,7 @@ public final class space/kscience/kmath/structures/IntBufferKt { } public final class space/kscience/kmath/structures/ListBuffer : space/kscience/kmath/structures/Buffer { + public fun (ILkotlin/jvm/functions/Function1;)V public fun (Ljava/util/List;)V public fun get (I)Ljava/lang/Object; public final fun getList ()Ljava/util/List; @@ -1865,6 +1866,7 @@ public final class space/kscience/kmath/structures/MutableBuffer$Companion { public final class space/kscience/kmath/structures/MutableListBuffer : space/kscience/kmath/structures/MutableBuffer { public static final synthetic fun box-impl (Ljava/util/List;)Lspace/kscience/kmath/structures/MutableListBuffer; + public static fun constructor-impl (ILkotlin/jvm/functions/Function1;)Ljava/util/List; public static fun constructor-impl (Ljava/util/List;)Ljava/util/List; public fun copy ()Lspace/kscience/kmath/structures/MutableBuffer; public static fun copy-impl (Ljava/util/List;)Lspace/kscience/kmath/structures/MutableBuffer; diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt index 551602115..d187beab1 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt @@ -194,6 +194,9 @@ public interface MutableBuffer : Buffer { * @property list The underlying list. */ public class ListBuffer(public val list: List) : Buffer { + + public constructor(size: Int, initializer: (Int) -> T) : this(List(size, initializer)) + override val size: Int get() = list.size override operator fun get(index: Int): T = list[index] @@ -213,8 +216,10 @@ public fun List.asBuffer(): ListBuffer = ListBuffer(this) */ @JvmInline public value class MutableListBuffer(public val list: MutableList) : MutableBuffer { - override val size: Int - get() = list.size + + public constructor(size: Int, initializer: (Int) -> T) : this(MutableList(size, initializer)) + + override val size: Int get() = list.size override operator fun get(index: Int): T = list[index] diff --git a/kmath-functions/build.gradle.kts b/kmath-functions/build.gradle.kts index 795137845..ca678bc0e 100644 --- a/kmath-functions/build.gradle.kts +++ b/kmath-functions/build.gradle.kts @@ -15,16 +15,23 @@ kotlin.sourceSets.commonMain { } readme { - description = "Functions and interpolation" - maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE + description = "Functions, integration and interpolation" + maturity = ru.mipt.npm.gradle.Maturity.EXPERIMENTAL propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md")) - feature("piecewise", "src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt", "Piecewise functions.") - feature("polynomials", "src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt", "Polynomial functions.") - feature("linear interpolation", - "src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt", - "Linear XY interpolator.") - feature("spline interpolation", - "src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt", - "Cubic spline XY interpolator.") + feature("piecewise", "src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt") { + "Piecewise functions." + } + feature("polynomials", "src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt") { + "Polynomial functions." + } + feature("linear interpolation", "src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt") { + "Linear XY interpolator." + } + feature("spline interpolation", "src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt") { + "Cubic spline XY interpolator." + } + feature("integration") { + "Univariate and multivariate quadratures" + } } \ No newline at end of file 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 9c94fa7bd..c4b9c572f 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 @@ -4,27 +4,31 @@ */ package space.kscience.kmath.integration +import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.operations.DoubleField -import space.kscience.kmath.operations.Ring -import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.indices +import space.kscience.kmath.operations.Field +import space.kscience.kmath.structures.* + /** * A simple one-pass integrator based on Gauss rule */ -public class GaussIntegrator> internal constructor( - public val algebra: Ring, - private val points: Buffer, - private val weights: Buffer, +public class GaussIntegrator( + public val algebra: Field, + public val bufferFactory: BufferFactory, ) : UnivariateIntegrator { - init { - require(points.size == weights.size) { "Inconsistent points and weights sizes" } - require(points.indices.all { i -> i == 0 || points[i] > points[i - 1] }) { "Integration nodes must be sorted" } + private fun buildRule(integrand: UnivariateIntegrand): Pair, Buffer> { + val factory = integrand.getFeature>() + ?: GenericGaussLegendreRuleFactory(algebra, bufferFactory) + val numPoints = integrand.getFeature()?.maxCalls ?: 100 + val range = integrand.getFeature>()?.range ?: 0.0..1.0 + return factory.build(numPoints, range) } override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand = with(algebra) { val f = integrand.function + val (points, weights) = buildRule(integrand) var res = zero var c = zero for (i in points.indices) { @@ -35,40 +39,73 @@ public class GaussIntegrator> internal constructor( c = t - res - y res = t } - return integrand + IntegrandValue(res) + IntegrandCalls(integrand.calls + points.size) + return integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + points.size) } public companion object { /** - * Integrate given [function] in a [range] with Gauss-Legendre quadrature with [numPoints] points. + * Integrate [T]-valued univariate function using provided set of [IntegrandFeature] + * Following features are evaluated: + * * [GaussIntegratorRuleFactory] - A factory for computing the Gauss integration rule. By default uses [GenericGaussLegendreRuleFactory] + * * [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 100 points. + */ + public fun integrate( + algebra: Field, + bufferFactory: BufferFactory = ::ListBuffer, + vararg features: IntegrandFeature, + function: (T) -> T, + ): UnivariateIntegrand = + GaussIntegrator(algebra, bufferFactory).integrate(UnivariateIntegrand(function, *features)) + + /** + * Integrate in real numbers */ public fun integrate( + vararg features: IntegrandFeature, + function: (Double) -> Double, + ): UnivariateIntegrand = integrate(DoubleField, ::DoubleBuffer, features = features, function) + + /** + * Integrate given [function] in a [range] with Gauss-Legendre quadrature with [numPoints] points. + * The [range] is automatically transformed into [T] using [Field.number] + */ + @UnstableKMathAPI + public fun legendre( + algebra: Field, range: ClosedRange, numPoints: Int = 100, - ruleFactory: GaussIntegratorRuleFactory = GaussLegendreDoubleRuleFactory, - features: List = emptyList(), - function: (Double) -> Double, - ): UnivariateIntegrand { - val (points, weights) = ruleFactory.build(numPoints, range) - return GaussIntegrator(DoubleField, points, weights).integrate( - UnivariateIntegrand(function, IntegrationRange(range)) + bufferFactory: BufferFactory = ::ListBuffer, + vararg features: IntegrandFeature, + function: (T) -> T, + ): UnivariateIntegrand = GaussIntegrator(algebra, bufferFactory).integrate( + UnivariateIntegrand( + function, + IntegrationRange(range), + DoubleGaussLegendreRuleFactory, + IntegrandMaxCalls(numPoints), + *features ) - } + ) -// public fun integrate( -// borders: List, -// numPoints: Int = 10, -// ruleFactory: GaussIntegratorRuleFactory = GaussLegendreDoubleRuleFactory, -// features: List = emptyList(), -// function: (Double) -> Double, -// ): UnivariateIntegrand { -// require(borders.indices.all { i -> i == 0 || borders[i] > borders[i - 1] }){"Borders are not sorted"} -// -// val (points, weights) = ruleFactory.build(numPoints, range) -// return GaussIntegrator(DoubleField, points, weights).integrate( -// UnivariateIntegrand(function, IntegrationRange(range)) -// ) -// } + /** + * Integrate given [function] in a [range] with Gauss-Legendre quadrature with [numPoints] points. + */ + @UnstableKMathAPI + public fun legendre( + range: ClosedRange, + numPoints: Int = 100, + vararg features: IntegrandFeature, + function: (Double) -> Double, + ): UnivariateIntegrand = GaussIntegrator(DoubleField, ::DoubleBuffer).integrate( + UnivariateIntegrand( + function, + IntegrationRange(range), + DoubleGaussLegendreRuleFactory, + IntegrandMaxCalls(numPoints), + *features + ) + ) } } \ No newline at end of file diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegratorRuleFactory.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegratorRuleFactory.kt index 235325dc4..8e961cc62 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegratorRuleFactory.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegratorRuleFactory.kt @@ -12,7 +12,7 @@ import kotlin.jvm.Synchronized import kotlin.math.ulp import kotlin.native.concurrent.ThreadLocal -public interface GaussIntegratorRuleFactory { +public interface GaussIntegratorRuleFactory : IntegrandFeature { public val algebra: Field public val bufferFactory: BufferFactory @@ -20,7 +20,7 @@ public interface GaussIntegratorRuleFactory { public companion object { public fun double(numPoints: Int, range: ClosedRange): Pair, Buffer> = - GaussLegendreDoubleRuleFactory.build(numPoints, range) + DoubleGaussLegendreRuleFactory.build(numPoints, range) } } @@ -28,16 +28,16 @@ public interface GaussIntegratorRuleFactory { * Create an integration rule by scaling existing normalized rule * */ -public fun > GaussIntegratorRuleFactory.build( +public fun GaussIntegratorRuleFactory.build( numPoints: Int, - range: ClosedRange, + range: ClosedRange, ): Pair, Buffer> { val normalized = build(numPoints) with(algebra) { val length = range.endInclusive - range.start val points = normalized.first.map(bufferFactory) { - range.start + length / 2 + length * it / 2 + number(range.start + length / 2) + number(length / 2) * it } val weights = normalized.second.map(bufferFactory) { @@ -56,7 +56,7 @@ public fun > GaussIntegratorRuleFactory.build( * */ @ThreadLocal -public object GaussLegendreDoubleRuleFactory : GaussIntegratorRuleFactory { +public object DoubleGaussLegendreRuleFactory : GaussIntegratorRuleFactory { override val algebra: Field get() = DoubleField @@ -173,3 +173,20 @@ public object GaussLegendreDoubleRuleFactory : GaussIntegratorRuleFactory, Buffer> = getOrBuildRule(numPoints) } + +/** + * A generic Gauss-Legendre rule factory that wraps [DoubleGaussLegendreRuleFactory] in a generic way. + */ +public class GenericGaussLegendreRuleFactory( + override val algebra: Field, + override val bufferFactory: BufferFactory, +) : GaussIntegratorRuleFactory { + + override fun build(numPoints: Int): Pair, Buffer> { + val (doublePoints, doubleWeights) = DoubleGaussLegendreRuleFactory.build(numPoints) + + val points = doublePoints.map(bufferFactory) { algebra.number(it) } + val weights = doubleWeights.map(bufferFactory) { algebra.number(it) } + return points to weights + } +} diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrand.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrand.kt index ef8efb1db..1ff8e422e 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrand.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrand.kt @@ -21,8 +21,8 @@ public class IntegrandRelativeAccuracy(public val accuracy: Double) : IntegrandF public class IntegrandAbsoluteAccuracy(public val accuracy: Double) : IntegrandFeature -public class IntegrandCalls(public val calls: Int) : IntegrandFeature +public class IntegrandCallsPerformed(public val calls: Int) : IntegrandFeature -public val Integrand.calls: Int get() = getFeature()?.calls ?: 0 +public val Integrand.calls: Int get() = getFeature()?.calls ?: 0 public class IntegrandMaxCalls(public val maxCalls: Int) : IntegrandFeature diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/GaussIntegralTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/GaussIntegralTest.kt index 00105c8fb..247318367 100644 --- a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/GaussIntegralTest.kt +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/GaussIntegralTest.kt @@ -13,7 +13,7 @@ import kotlin.test.assertEquals class GaussIntegralTest { @Test fun gaussSin() { - val res = GaussIntegrator.integrate(0.0..2 * PI) { x -> + val res = GaussIntegrator.legendre(0.0..2 * PI) { x -> sin(x) } assertEquals(0.0, res.value!!, 1e-4) @@ -21,7 +21,7 @@ class GaussIntegralTest { @Test fun gaussUniform() { - val res = GaussIntegrator.integrate(0.0..100.0,300) { x -> + val res = GaussIntegrator.legendre(0.0..100.0,300) { x -> if(x in 30.0..50.0){ 1.0 } else {