[WIP] Interpolation fix

examples/build.gradle.kts

@@ -43,7 +43,7 @@ dependencies {
 
     implementation("org.slf4j:slf4j-simple:1.7.30")
     // plotting
-    implementation("space.kscience:plotlykt-server:0.4.0-dev-2")
+    implementation("space.kscience:plotlykt-server:0.4.0")
 }
 
 kotlin.sourceSets.all {

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

@@ -0,0 +1,56 @@
+/*
+ * 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.structures.DoubleBuffer
+import space.kscience.plotly.Plotly
+import space.kscience.plotly.makeFile
+import space.kscience.plotly.models.functionXY
+import space.kscience.plotly.scatter
+import kotlin.math.PI
+import kotlin.math.sin
+
+fun main() {
+    val data = (0..10).map {
+        val x = it.toDouble() / 5 * PI
+        x to sin(x)
+    }
+
+    val polynomial: PiecewisePolynomial<Double> = SplineInterpolator(
+        DoubleField, ::DoubleBuffer
+    ).interpolatePolynomials(data)
+
+    val function = polynomial.asFunction(DoubleField, 0.0)
+
+    val cmInterpolate = org.apache.commons.math3.analysis.interpolation.SplineInterpolator().interpolate(
+        data.map { it.first }.toDoubleArray(),
+        data.map { it.second }.toDoubleArray()
+    )
+
+    println(function(2.0))
+    println(cmInterpolate.value(2.0))
+
+//
+//    Plotly.plot {
+//        scatter {
+//            name = "interpolated"
+//            x.numbers = data.map { it.first }
+//            y.numbers = x.doubles.map { function(it) }
+//        }
+//        scatter {
+//            name = "original"
+//            functionXY(0.0..(2 * PI), 0.1) { sin(it) }
+//        }
+//        scatter {
+//            name = "cm"
+//            x.numbers = data.map { it.first }
+//            y.numbers = x.doubles.map { cmInterpolate.value(it) }
+//        }
+//    }.makeFile()
+}

kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt

@@ -8,28 +8,18 @@
 package space.kscience.kmath.ejml
 
 import org.ejml.data.*
-import org.ejml.dense.row.CommonOps_DDRM
-import org.ejml.dense.row.CommonOps_FDRM
-import org.ejml.dense.row.factory.DecompositionFactory_DDRM
-import org.ejml.dense.row.factory.DecompositionFactory_FDRM
-import org.ejml.sparse.FillReducing
-import org.ejml.sparse.csc.CommonOps_DSCC
-import org.ejml.sparse.csc.CommonOps_FSCC
-import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC
-import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC
-import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC
-import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC
 import space.kscience.kmath.linear.*
+import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.*
+import space.kscience.kmath.misc.*
+import kotlin.reflect.*
+import org.ejml.dense.row.*
+import org.ejml.dense.row.factory.*
+import org.ejml.sparse.*
+import org.ejml.sparse.csc.*
+import org.ejml.sparse.csc.factory.*
+import space.kscience.kmath.nd.*
 import space.kscience.kmath.linear.Matrix
-import space.kscience.kmath.misc.UnstableKMathAPI
-import space.kscience.kmath.nd.StructureFeature
-import space.kscience.kmath.operations.DoubleField
-import space.kscience.kmath.operations.FloatField
-import space.kscience.kmath.operations.invoke
-import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.FloatBuffer
-import kotlin.reflect.KClass
-import kotlin.reflect.cast
 
 /**
  * [EjmlVector] specialization for [Double].

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

@@ -26,12 +26,12 @@ public fun interface Piecewise<T, R> {
 public fun interface PiecewisePolynomial<T : Any> : Piecewise<T, Polynomial<T>>
 
 /**
- * Basic [Piecewise] implementation where all the pieces are ordered by the [Comparable] type instances.
+ * A [Piecewise]  builder where all the pieces are ordered by the [Comparable] type instances.
  *
  * @param T the comparable piece key type.
+ * @param delimiter the initial piecewise separator
  */
-public class OrderedPiecewisePolynomial<T : Comparable<T>>(delimiter: T) :
+public class OrderedPiecewisePolynomial<T : Comparable<T>>(delimiter: T) : PiecewisePolynomial<T> {
-    PiecewisePolynomial<T> {
     private val delimiters: MutableList<T> = arrayListOf(delimiter)
     private val pieces: MutableList<Polynomial<T>> = arrayListOf()
 
@@ -64,9 +64,7 @@ public class OrderedPiecewisePolynomial<T : Comparable<T>>(delimiter: T) :
             return null
         else {
             for (index in 1 until delimiters.size)
-                if (arg < delimiters[index])
+                if (arg < delimiters[index]) return pieces[index - 1]
-                    return pieces[index - 1]
-
             error("Piece not found")
         }
     }

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

@@ -19,7 +19,9 @@ import kotlin.math.pow
  *
  * @param coefficients constant is the leftmost coefficient.
  */
-public class Polynomial<T : Any>(public val coefficients: List<T>)
+public class Polynomial<T : Any>(public val coefficients: List<T>){
+    override fun toString(): String = "Polynomial$coefficients"
+}
 
 /**
  * Returns a [Polynomial] instance with given [coefficients].
@@ -34,19 +36,15 @@ public fun Polynomial<Double>.value(): Double = coefficients.reduceIndexed { ind
 
 /**
  * Evaluates the value of the given polynomial for given argument.
+ * https://en.wikipedia.org/wiki/Horner%27s_method
  */
 public fun <T : Any, C : Ring<T>> Polynomial<T>.value(ring: C, arg: T): T = ring {
     if (coefficients.isEmpty()) return@ring zero
-    var res = coefficients.first()
+    var result: T = coefficients.last()
-    var powerArg = arg
+    for (j in coefficients.size - 2 downTo 0) {
-
+        result = (arg * result) + coefficients[j]
-    for (index in 1 until coefficients.size) {
-        res += coefficients[index] * powerArg
-        // recalculating power on each step to avoid power costs on long polynomials
-        powerArg *= arg
     }
-
+    return result
-    res
 }
 
 /**

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt

@@ -34,19 +34,16 @@ public class SplineInterpolator<T : Comparable<T>>(
         // Number of intervals.  The number of data points is n + 1.
         val n = points.size - 1
         // Differences between knot points
-        val h = bufferFactory(points.size) { i -> points.x[i + 1] - points.x[i] }
+        val h = bufferFactory(n) { i -> points.x[i + 1] - points.x[i] }
-        val mu = bufferFactory(points.size - 1) { zero }
+        val mu = bufferFactory(n) { zero }
-        val z = bufferFactory(points.size) { zero }
+        val z = bufferFactory(n + 1) { zero }
 
         for (i in 1 until n) {
             val g = 2.0 * (points.x[i + 1] - points.x[i - 1]) - h[i - 1] * mu[i - 1]
             mu[i] = h[i] / g
-
+            z[i] =
-            z[i] = (3.0 * (points.y[i + 1] * h[i - 1]
+                ((points.y[i + 1] * h[i - 1] - points.y[i] * (points.x[i + 1] - points.x[i - 1]) + points.y[i - 1] * h[i]) * 3.0 /
-                    - points.x[i] * (points.x[i + 1] - points.x[i - 1])
+                        (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g
-                    + points.y[i - 1] * h[i]) / (h[i - 1] * h[i])
-                    - h[i - 1] * z[i - 1]
-                    ) / g
         }
 
         // cubic spline coefficients --  b is linear, c quadratic, d is cubic (original y's are constants)

kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/SplineInterpolatorTest.kt

@@ -0,0 +1,35 @@
+/*
+ * 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.interpolation
+
+import space.kscience.kmath.functions.PiecewisePolynomial
+import space.kscience.kmath.functions.asFunction
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.structures.DoubleBuffer
+import kotlin.math.PI
+import kotlin.math.sin
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+internal class SplineInterpolatorTest {
+    @Test
+    fun testInterpolation() {
+        val data = (0..10).map {
+            val x = it.toDouble() / 5 * PI
+            x to sin(x)
+        }
+
+        val polynomial: PiecewisePolynomial<Double> = SplineInterpolator(
+            DoubleField, ::DoubleBuffer
+        ).interpolatePolynomials(data)
+
+        val function = polynomial.asFunction(DoubleField)
+        assertEquals(null, function(-1.0))
+        assertEquals(sin(0.5), function(0.5)!!, 0.1)
+        assertEquals(sin(1.5), function(1.5)!!, 0.1)
+        assertEquals(sin(2.0), function(2.0)!!, 0.1)
+    }
+}