From d45d44f96d87596594b9f60935b742f2a407b6df Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Mon, 17 May 2021 21:55:32 +0300 Subject: [PATCH] Fix Interpolation --- .../kscience/kmath/functions/interpolate.kt | 35 +++++++++---------- .../space/kscience/kmath/ejml/_generated.kt | 30 ++++++++++------ .../kmath/interpolation/SplineInterpolator.kt | 19 ++++++++-- .../interpolation/SplineInterpolatorTest.kt | 15 ++++---- 4 files changed, 59 insertions(+), 40 deletions(-) diff --git a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt b/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt index 54120ac50..69ba80fb6 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt @@ -33,24 +33,21 @@ fun main() { 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() + 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() } \ No newline at end of file diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt index 6f74ab24f..139c55697 100644 --- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt +++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt @@ -8,18 +8,28 @@ 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]. diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt index fc341e019..c3d6ffae7 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt @@ -10,8 +10,10 @@ import space.kscience.kmath.functions.OrderedPiecewisePolynomial import space.kscience.kmath.functions.PiecewisePolynomial import space.kscience.kmath.functions.Polynomial 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.structures.DoubleBuffer import space.kscience.kmath.structures.MutableBufferFactory /** @@ -56,10 +58,23 @@ public class SplineInterpolator>( val a = points.y[j] val b = (points.y[j + 1] - points.y[j]) / h[j] - h[j] * (cOld + 2.0 * c) / 3.0 val d = (cOld - c) / (3.0 * h[j]) - val polynomial = Polynomial(a, b, c, d) + val x0 = points.x[j] + val x02 = x0 * x0 + val x03 = x02 * x0 + //Shift coefficients to represent absolute polynomial instead of one with an offset + val polynomial = Polynomial( + a - b * x0 + c * x02 - d * x03, + b - 2*c*x0 + 3*d*x02, + c - 3*d*x0, + d + ) cOld = c - putLeft(points.x[j], polynomial) + putLeft(x0, polynomial) } } } + + public companion object { + public val double: SplineInterpolator = SplineInterpolator(DoubleField, ::DoubleBuffer) + } } diff --git a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/SplineInterpolatorTest.kt b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/SplineInterpolatorTest.kt index c3bd8195a..3adaab2d1 100644 --- a/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/SplineInterpolatorTest.kt +++ b/kmath-functions/src/commonTest/kotlin/space/kscience/kmath/interpolation/SplineInterpolatorTest.kt @@ -8,7 +8,6 @@ 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 @@ -22,14 +21,12 @@ internal class SplineInterpolatorTest { x to sin(x) } - val polynomial: PiecewisePolynomial = SplineInterpolator( - DoubleField, ::DoubleBuffer - ).interpolatePolynomials(data) + val polynomial: PiecewisePolynomial = SplineInterpolator.double.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) + val function = polynomial.asFunction(DoubleField, Double.NaN) + assertEquals(Double.NaN, 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) } }