Simpson integrator

examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt

@@ -35,7 +35,6 @@ fun main() {
         data.map { it.second }.toDoubleArray()
     )
 
-
     Plotly.plot {
         scatter {
             name = "interpolated"

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMGaussRuleIntegrator.kt

@@ -11,7 +11,7 @@ import space.kscience.kmath.integration.*
 /**
  * A simple one-pass integrator based on Gauss rule
  */
-public class GaussRuleIntegrator(
+public class CMGaussRuleIntegrator(
     private val numpoints: Int,
     private var type: GaussRule = GaussRule.LEGANDRE,
 ) : UnivariateIntegrator<Double> {
@@ -76,7 +76,7 @@ public class GaussRuleIntegrator(
             numPoints: Int = 100,
             type: GaussRule = GaussRule.LEGANDRE,
             function: (Double) -> Double,
-        ): Double = GaussRuleIntegrator(numPoints, type).integrate(
+        ): Double = CMGaussRuleIntegrator(numPoints, type).integrate(
             UnivariateIntegrand(function, IntegrationRange(range))
         ).valueOrNull!!
     }

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

@@ -5,13 +5,37 @@
 
 package space.kscience.kmath.integration
 
-import space.kscience.kmath.operations.Ring
-
+import space.kscience.kmath.operations.Field
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.BufferFactory
 
+/**
+ * [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: Ring<T>,
+    public val algebra: Field<T>,
+    public val bufferFactory: BufferFactory<T>,
 ) : UnivariateIntegrator<T> {
-    override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> {
-        TODO("Not yet implemented")
+    override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = with(algebra) {
+        val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
+        require(numPoints >= 4)
+        val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
+        val h: Double = (range.endInclusive - range.start) / (numPoints - 1)
+        val points: Buffer<T> = bufferFactory(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])
+        var res = zero
+        res += simpson(1) / 1.5 //border points with 1.5 factor
+        for (i in 2 until (points.size - 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)
     }
 }
+
+public inline val <reified T : Any> Field<T>.simpsonIntegrator: SimpsonIntegrator<T>
+    get() = SimpsonIntegrator(this, Buffer.Companion::auto)

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

@@ -0,0 +1,12 @@
+/*
+ * 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
+
+public class SplineIntegrator<T>: UnivariateIntegrator<T> {
+    override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> {
+        TODO("Not yet implemented")
+    }
+}

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

@@ -9,7 +9,7 @@ import space.kscience.kmath.misc.UnstableKMathAPI
 import kotlin.jvm.JvmInline
 import kotlin.reflect.KClass
 
-public class UnivariateIntegrand<T : Any> internal constructor(
+public class UnivariateIntegrand<T> internal constructor(
     private val features: Map<KClass<*>, IntegrandFeature>,
     public val function: (Double) -> T,
 ) : Integrand {

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

@@ -19,7 +19,7 @@ class GaussIntegralTest {
         val res = DoubleField.gaussIntegrator.integrate(0.0..2 * PI) { x ->
             sin(x)
         }
-        assertEquals(0.0, res.valueOrNull!!, 1e-2)
+        assertEquals(0.0, res.value, 1e-2)
     }
 
     @Test
@@ -31,7 +31,7 @@ class GaussIntegralTest {
                 0.0
             }
         }
-        assertEquals(15.0, res.valueOrNull!!, 0.5)
+        assertEquals(15.0, res.value, 0.5)
     }
 
 

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

@@ -0,0 +1,38 @@
+/*
+ * 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.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 SimpsonIntegralTest {
+    @Test
+    fun gaussSin() {
+        val res = DoubleField.simpsonIntegrator.integrate(0.0..2 * PI, IntegrandMaxCalls(5)) { x ->
+            sin(x)
+        }
+        assertEquals(0.0, res.value, 1e-2)
+    }
+
+    @Test
+    fun gaussUniform() {
+        val res = DoubleField.simpsonIntegrator.integrate(35.0..100.0) { x ->
+            if(x in 30.0..50.0){
+                1.0
+            } else {
+                0.0
+            }
+        }
+        assertEquals(15.0, res.value, 0.5)
+    }
+
+
+}