Generalized linear interpolation

kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Piecewise.kt

@@ -0,0 +1,55 @@
+package scientifik.kmath.functions
+
+import scientifik.kmath.operations.Ring
+
+interface Piecewise<T, R> {
+    fun findPiece(arg: T): R?
+}
+
+interface PiecewisePolynomial<T : Any> : Piecewise<T, Polynomial<T>>
+
+/**
+ * Ordered list of pieces in piecewise function
+ */
+class OrderedPiecewisePolynomial<T : Comparable<T>>(left: T) : PiecewisePolynomial<T> {
+
+    private val delimiters: ArrayList<T> = arrayListOf(left)
+    private val pieces: ArrayList<Polynomial<T>> = ArrayList()
+
+    /**
+     * Dynamically add a piece to the "right" side (beyond maximum argument value of previous piece)
+     * @param right new rightmost position. If is less then current rightmost position, a error is thrown.
+     */
+    fun putRight(right: T, piece: Polynomial<T>) {
+        require(right > delimiters.last()) { "New delimiter should be to the right of old one" }
+        delimiters.add(right)
+        pieces.add(piece)
+    }
+
+    fun putLeft(left: T, piece: Polynomial<T>) {
+        require(left < delimiters.first()) { "New delimiter should be to the left of old one" }
+        delimiters.add(0, left)
+        pieces.add(0, piece)
+    }
+
+    override fun findPiece(arg: T): Polynomial<T>? {
+        if (arg < delimiters.first() || arg >= delimiters.last()) {
+            return null
+        } else {
+            for (index in 1 until delimiters.size) {
+                if (arg < delimiters[index]) {
+                    return pieces[index - 1]
+                }
+            }
+            error("Piece not found")
+        }
+    }
+}
+
+/**
+ * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition.
+ */
+fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.value(ring: C, arg: T): T? =
+    findPiece(arg)?.value(ring, arg)
+
+fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { value(ring, it) }

kmath-functions/src/commonMain/kotlin/scientifik.kmath.functions/Polynomial.kt

@@ -1,9 +1,7 @@
 package scientifik.kmath.functions
 
-import scientifik.kmath.operations.RealField
 import scientifik.kmath.operations.Ring
 import scientifik.kmath.operations.Space
-import kotlin.jvm.JvmName
 import kotlin.math.max
 import kotlin.math.pow
 
@@ -20,11 +18,11 @@ fun Polynomial<Double>.value() =
 
 
 fun <T : Any, C : Ring<T>> Polynomial<T>.value(ring: C, arg: T): T = ring.run {
-    if( coefficients.isEmpty()) return@run zero
+    if (coefficients.isEmpty()) return@run zero
     var res = coefficients.first()
     var powerArg = arg
-    for( index in 1 until coefficients.size){
+    for (index in 1 until coefficients.size) {
-        res += coefficients[index]*powerArg
+        res += coefficients[index] * powerArg
         //recalculating power on each step to avoid power costs on long polynomials
         powerArg *= arg
     }
@@ -43,9 +41,6 @@ fun <T : Any, C : Ring<T>> Polynomial<T>.asMathFunction(): MathFunction<T, out C
  */
 fun <T : Any, C : Ring<T>> Polynomial<T>.asFunction(ring: C): (T) -> T = { value(ring, it) }
 
-@JvmName("asRealUFunction")
-fun Polynomial<Double>.asFunction(): (Double) -> Double = asFunction(RealField)
-
 /**
  * An algebra for polynomials
  */
@@ -74,21 +69,3 @@ class PolynomialSpace<T : Any, C : Ring<T>>(val ring: C) : Space<Polynomial<T>>
 fun <T : Any, C : Ring<T>, R> C.polynomial(block: PolynomialSpace<T, C>.() -> R): R {
     return PolynomialSpace(this).run(block)
 }
-
-class PiecewisePolynomial<T : Comparable<T>> internal constructor(
-    val lowerBoundary: T,
-    val pieces: List<Pair<T, Polynomial<T>>>
-)
-
-private fun <T : Comparable<T>> PiecewisePolynomial<T>.findPiece(arg: T): Polynomial<T>? {
-    if (arg < lowerBoundary || arg > pieces.last().first) return null
-    return pieces.first { arg < it.first }.second
-}
-
-/**
- * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition.
- */
-fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.value(ring: C, arg: T): T? =
-    findPiece(arg)?.value(ring, arg)
-
-fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.asFunction(ring: C): (T) -> T? = { value(ring, it) }

kmath-functions/src/commonMain/kotlin/scientifik/kmath/interpolation/LinearInterpolator.kt

@@ -1,5 +1,6 @@
 package scientifik.kmath.interpolation
 
+import scientifik.kmath.functions.OrderedPiecewisePolynomial
 import scientifik.kmath.functions.PiecewisePolynomial
 import scientifik.kmath.functions.Polynomial
 import scientifik.kmath.operations.Field
@@ -10,15 +11,18 @@ import scientifik.kmath.operations.Field
 class LinearInterpolator<T : Comparable<T>>(override val algebra: Field<T>) : PolynomialInterpolator<T> {
 
     override fun interpolatePolynomials(points: Collection<Pair<T, T>>): PiecewisePolynomial<T> = algebra.run {
+        require(points.isNotEmpty()) { "Point array should not be empty" }
+
         //sorting points
         val sorted = points.sortedBy { it.first }
 
-        val pairs: List<Pair<T, Polynomial<T>>> = (0 until points.size - 1).map { i ->
+        return@run OrderedPiecewisePolynomial(points.first().first).apply {
+            for (i in 0 until points.size - 1) {
                 val slope = (sorted[i + 1].second - sorted[i].second) / (sorted[i + 1].first - sorted[i].first)
                 val const = sorted[i].second - slope * sorted[i].first
-            sorted[i + 1].first to Polynomial(const, slope)
+                val polynomial = Polynomial(const, slope)
+                putRight(sorted[i + 1].first, polynomial)
+            }
         }
-
-        return PiecewisePolynomial(sorted.first().first, pairs)
     }
 }

kmath-functions/src/commonTest/kotlin/scientifik/kmath/interpolation/LinearInterpolatorTest.kt

@@ -18,8 +18,8 @@ class LinearInterpolatorTest {
         val polynomial = LinearInterpolator(RealField).interpolatePolynomials(data)
         val function = polynomial.asFunction(RealField)
 
-//        assertEquals(null, function(-1.0))
+        assertEquals(null, function(-1.0))
-//        assertEquals(0.5, function(0.5))
+        assertEquals(0.5, function(0.5))
         assertEquals(2.0, function(1.5))
         assertEquals(3.0, function(2.0))
     }