Simpson and spline integration

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

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/algebraExtensions.kt

@@ -16,6 +16,8 @@ public fun <T> Ring<T>.sum(data: Iterable<T>): 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].
  *

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<StructureND<out Number>, Double> {
+    object L2Norm : Norm<StructureND<Number>, Double> {
         @OptIn(PerformancePitfall::class)
-        override fun norm(arg: StructureND<out Number>): Double =
+        override fun norm(arg: StructureND<Number>): Double =
             kotlin.math.sqrt(arg.elements().sumOf { it.second.toDouble() })
     }
 

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt

@@ -69,18 +69,23 @@ public fun <T, A> Polynomial<T>.differentiate(
 public fun <T, A> Polynomial<T>.integrate(
     algebra: A,
 ): Polynomial<T> where  A : Field<T>, A : NumericAlgebra<T> = algebra {
-    Polynomial(coefficients.mapIndexed { index, t -> t / number(index) })
+    val integratedCoefficients = buildList<T>(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 <T : Comparable<T>, A> Polynomial<T>.integrate(
+public fun <T : Comparable<T>> Polynomial<T>.integrate(
-    algebra: A,
+    algebra: Field<T>,
     range: ClosedRange<T>,
-): T where  A : Field<T>, A : NumericAlgebra<T> = algebra {
+): T = algebra {
-    value(algebra, range.endInclusive) - value(algebra, range.start)
+    val integral = integrate(algebra)
+    integral.value(algebra, range.endInclusive) - integral.value(algebra, range.start)
 }
 
 /**

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<Pair<ClosedRange<Double>, Int>>) : IntegrandFeature {
-    public constructor(vararg pairs: Pair<ClosedRange<Double>, Int>) : this(pairs.toList())
-}
 
 /**
  * A simple one-pass integrator based on Gauss rule

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.operations.invoke
-import space.kscience.kmath.structures.BufferFactory
+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<T : Any>(
     public val algebra: Field<T>,
-    public val bufferFactory: BufferFactory<T>,
 ) : UnivariateIntegrator<T> {
-    override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = with(algebra) {
+
-        val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
+    private fun integrateRange(
-        require(numPoints >= 4)
+        integrand: UnivariateIntegrand<T>, range: ClosedRange<Double>, numPoints: Int,
-        val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
+    ): T = algebra {
         val h: Double = (range.endInclusive - range.start) / (numPoints - 1)
-        val points: Buffer<T> = bufferFactory(numPoints) { i ->
+        val values: List<T> = 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
+        res += simpson(values.size - 2) / 1.5 //border points with 1.5 factor
-        return integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + points.size)
+        return res
+    }
+
+    override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> {
+        val ranges = integrand.getFeature<UnivariateIntegrandRanges>()
+        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<IntegrandMaxCalls>()?.maxCalls ?: 100
+            require(numPoints >= 4) { "Simpson integrator requires at least 4 nodes" }
+            val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
+            val res = integrateRange(integrand, range, numPoints)
+            integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + numPoints)
+        }
     }
 }
 
-public inline val <reified T : Any> Field<T>.simpsonIntegrator: SimpsonIntegrator<T>
+public val <T : Any> Field<T>.simpsonIntegrator: SimpsonIntegrator<T> get() = SimpsonIntegrator(this)
-    get() = SimpsonIntegrator(this, Buffer.Companion::auto)

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt

@@ -5,8 +5,98 @@
 
 package space.kscience.kmath.integration
 
-public class SplineIntegrator<T>: UnivariateIntegrator<T> {
+import space.kscience.kmath.functions.PiecewisePolynomial
-    override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> {
+import space.kscience.kmath.functions.integrate
-        TODO("Not yet implemented")
+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 <T : Comparable<T>> PiecewisePolynomial<T>.integrate(algebra: Field<T>): PiecewisePolynomial<T> =
+    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 <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] - 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<T : Comparable<T>>(
+    public val algebra: Field<T>,
+    public val bufferFactory: MutableBufferFactory<T>,
+) : UnivariateIntegrator<T> {
+    override fun integrate(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] - 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<Double> {
+    override fun integrate(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)
+    }
+}
+
+@UnstableKMathAPI
+public inline val DoubleField.splineIntegrator: UnivariateIntegrator<Double>
+    get() = DoubleSplineIntegrator

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<T> = Integrator<UnivariateIntegrand<T>>
 @JvmInline
 public value class IntegrationRange(public val range: ClosedRange<Double>) : 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<Pair<ClosedRange<Double>, Int>>) : IntegrandFeature {
+    public constructor(vararg pairs: Pair<ClosedRange<Double>, Int>) : this(pairs.toList())
+}
+
+public class UnivariateIntegrationNodes(public val nodes: Buffer<Double>) : IntegrandFeature {
+    public constructor(vararg nodes: Double) : this(DoubleBuffer(nodes))
+}
+
+
 /**
  * Value of the integrand if it is present or null
  */

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<T, 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>
 

kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SimpsonIntegralTest.kt

@@ -24,7 +24,7 @@ class SimpsonIntegralTest {
 
     @Test
     fun gaussUniform() {
-        val res = DoubleField.simpsonIntegrator.integrate(35.0..100.0) { x ->
+        val res = DoubleField.simpsonIntegrator.integrate(35.0..100.0, IntegrandMaxCalls(20)) { x ->
             if (x in 30.0..50.0) {
                 1.0
             } else {
@@ -33,6 +33,4 @@ class SimpsonIntegralTest {
         }
         assertEquals(15.0, res.value, 0.5)
     }
-
-
 }

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

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"
+        id("org.jetbrains.kotlinx.benchmark") version "0.3.1"
-        kotlin("jupyter.api") version "0.9.1-61"
+        kotlin("jupyter.api") version "0.10.0-25"
 
     }
 }