Field element integration

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

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{
+    fun jvmAdd(blackhole: Blackhole) = JBigIntegerField {
-        blackhole.consume(jvmNumber + jvmNumber+ jvmNumber)
+        blackhole.consume(jvmNumber + jvmNumber + jvmNumber)
     }
 
     @Benchmark
-    fun kmMultiply(blackhole: Blackhole) = BigIntField{
+    fun kmMultiply(blackhole: Blackhole) = BigIntField {
-        blackhole.consume(kmNumber*kmNumber*kmNumber)
+        blackhole.consume(kmNumber * kmNumber * kmNumber)
     }
 
     @Benchmark
-    fun jvmMultiply(blackhole: Blackhole) = JBigIntegerField{
+    fun jvmMultiply(blackhole: Blackhole) = JBigIntegerField {
-        blackhole.consume(jvmNumber*jvmNumber*jvmNumber)
+        blackhole.consume(jvmNumber * jvmNumber * jvmNumber)
     }
 }

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<Double> = { 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)
+}

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<StructureND<Double>> = { x -> 3 * x.pow(2) + 2 * x + 1 }
+
+        //get the result of the integration
+        val result: UnivariateIntegrand<StructureND<Double>> = 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)
+    }
+}

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)
     }
 
 

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<Double>): GaussIntegrator {

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 <init> (ILkotlin/jvm/functions/Function1;)V
 	public fun <init> (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;

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt

@@ -194,6 +194,9 @@ public interface MutableBuffer<T> : Buffer<T> {
  * @property list The underlying list.
  */
 public class ListBuffer<T>(public val list: List<T>) : Buffer<T> {
+
+    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 <T> List<T>.asBuffer(): ListBuffer<T> = ListBuffer(this)
  */
 @JvmInline
 public value class MutableListBuffer<T>(public val list: MutableList<T>) : MutableBuffer<T> {
-    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]
 

kmath-functions/build.gradle.kts

@@ -15,16 +15,23 @@ kotlin.sourceSets.commonMain {
 }
 
 readme {
-    description = "Functions and interpolation"
+    description = "Functions, integration and interpolation"
-    maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE
+    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("piecewise", "src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt") {
-    feature("polynomials", "src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt", "Polynomial functions.")
+        "Piecewise functions."
-    feature("linear interpolation",
+    }
-        "src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt",
+    feature("polynomials", "src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt") {
-        "Linear XY interpolator.")
+        "Polynomial functions."
-    feature("spline interpolation",
+    }
-        "src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt",
+    feature("linear interpolation", "src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt") {
-        "Cubic spline XY interpolator.")
+        "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"
+    }
 }

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.operations.Field
-import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.*
-import space.kscience.kmath.structures.indices
+
 
 /**
  * A simple one-pass integrator based on Gauss rule
  */
-public class GaussIntegrator<T : Comparable<T>> internal constructor(
+public class GaussIntegrator<T : Any>(
-    public val algebra: Ring<T>,
+    public val algebra: Field<T>,
-    private val points: Buffer<T>,
+    public val bufferFactory: BufferFactory<T>,
-    private val weights: Buffer<T>,
 ) : UnivariateIntegrator<T> {
 
-    init {
+    private fun buildRule(integrand: UnivariateIntegrand<T>): Pair<Buffer<T>, Buffer<T>> {
-        require(points.size == weights.size) { "Inconsistent points and weights sizes" }
+        val factory = integrand.getFeature<GaussIntegratorRuleFactory<T>>()
-        require(points.indices.all { i -> i == 0 || points[i] > points[i - 1] }) { "Integration nodes must be sorted" }
+            ?: GenericGaussLegendreRuleFactory(algebra, bufferFactory)
+        val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
+        val range = integrand.getFeature<IntegrationRange<Double>>()?.range ?: 0.0..1.0
+        return factory.build(numPoints, range)
     }
 
     override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = 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<T : Comparable<T>> 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 <T : Any> integrate(
+            algebra: Field<T>,
+            bufferFactory: BufferFactory<T> = ::ListBuffer,
+            vararg features: IntegrandFeature,
+            function: (T) -> T,
+        ): UnivariateIntegrand<T> =
+            GaussIntegrator(algebra, bufferFactory).integrate(UnivariateIntegrand(function, *features))
+
+        /**
+         * Integrate in real numbers
          */
         public fun integrate(
+            vararg features: IntegrandFeature,
+            function: (Double) -> Double,
+        ): UnivariateIntegrand<Double> = 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  <T : Any> legendre(
+            algebra: Field<T>,
             range: ClosedRange<Double>,
             numPoints: Int = 100,
-            ruleFactory: GaussIntegratorRuleFactory<Double> = GaussLegendreDoubleRuleFactory,
+            bufferFactory: BufferFactory<T> = ::ListBuffer,
-            features: List<IntegrandFeature> = emptyList(),
+            vararg features: IntegrandFeature,
+            function: (T) -> T,
+        ): UnivariateIntegrand<T> = GaussIntegrator(algebra, bufferFactory).integrate(
+            UnivariateIntegrand(
+                function,
+                IntegrationRange(range),
+                DoubleGaussLegendreRuleFactory,
+                IntegrandMaxCalls(numPoints),
+                *features
+            )
+        )
+
+        /**
+         * Integrate given [function] in a [range] with Gauss-Legendre quadrature with [numPoints] points.
+         */
+        @UnstableKMathAPI
+        public fun legendre(
+            range: ClosedRange<Double>,
+            numPoints: Int = 100,
+            vararg features: IntegrandFeature,
             function: (Double) -> Double,
-        ): UnivariateIntegrand<Double> {
+        ): UnivariateIntegrand<Double> = GaussIntegrator(DoubleField, ::DoubleBuffer).integrate(
-            val (points, weights) = ruleFactory.build(numPoints, range)
+            UnivariateIntegrand(
-            return GaussIntegrator(DoubleField, points, weights).integrate(
+                function,
-                UnivariateIntegrand(function, IntegrationRange(range))
+                IntegrationRange(range),
+                DoubleGaussLegendreRuleFactory,
+                IntegrandMaxCalls(numPoints),
+                *features
+            )
         )
     }
-
-//        public fun integrate(
-//            borders: List<Double>,
-//            numPoints: Int = 10,
-//            ruleFactory: GaussIntegratorRuleFactory<Double> = GaussLegendreDoubleRuleFactory,
-//            features: List<IntegrandFeature> = emptyList(),
-//            function: (Double) -> Double,
-//        ): UnivariateIntegrand<Double> {
-//            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))
-//            )
-//        }
-    }
 }

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<T : Any> {
+public interface GaussIntegratorRuleFactory<T : Any> : IntegrandFeature {
     public val algebra: Field<T>
     public val bufferFactory: BufferFactory<T>
 
@@ -20,7 +20,7 @@ public interface GaussIntegratorRuleFactory<T : Any> {
 
     public companion object {
         public fun double(numPoints: Int, range: ClosedRange<Double>): Pair<Buffer<Double>, Buffer<Double>> =
-            GaussLegendreDoubleRuleFactory.build(numPoints, range)
+            DoubleGaussLegendreRuleFactory.build(numPoints, range)
     }
 }
 
@@ -28,16 +28,16 @@ public interface GaussIntegratorRuleFactory<T : Any> {
  * Create an integration rule by scaling existing normalized rule
  *
  */
-public fun <T : Comparable<T>> GaussIntegratorRuleFactory<T>.build(
+public fun <T : Any> GaussIntegratorRuleFactory<T>.build(
     numPoints: Int,
-    range: ClosedRange<T>,
+    range: ClosedRange<Double>,
 ): Pair<Buffer<T>, Buffer<T>> {
     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 <T : Comparable<T>> GaussIntegratorRuleFactory<T>.build(
  *
  */
 @ThreadLocal
-public object GaussLegendreDoubleRuleFactory : GaussIntegratorRuleFactory<Double> {
+public object DoubleGaussLegendreRuleFactory : GaussIntegratorRuleFactory<Double> {
 
     override val algebra: Field<Double> get() = DoubleField
 
@@ -173,3 +173,20 @@ public object GaussLegendreDoubleRuleFactory : GaussIntegratorRuleFactory<Double
     override fun build(numPoints: Int): Pair<Buffer<Double>, Buffer<Double>> = getOrBuildRule(numPoints)
 }
 
+
+/**
+ * A generic Gauss-Legendre rule factory that wraps [DoubleGaussLegendreRuleFactory] in a generic way.
+ */
+public class GenericGaussLegendreRuleFactory<T : Any>(
+    override val algebra: Field<T>,
+    override val bufferFactory: BufferFactory<T>,
+) : GaussIntegratorRuleFactory<T> {
+
+    override fun build(numPoints: Int): Pair<Buffer<T>, Buffer<T>> {
+        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
+    }
+}

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<IntegrandCalls>()?.calls ?: 0
+public val Integrand.calls: Int get() = getFeature<IntegrandCallsPerformed>()?.calls ?: 0
 
 public class IntegrandMaxCalls(public val maxCalls: Int) : IntegrandFeature

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 {