Fix Interpolation

This commit is contained in:
Alexander Nozik 2021-05-17 21:55:32 +03:00
parent 53a6d3543f
commit d45d44f96d
4 changed files with 59 additions and 40 deletions

View File

@ -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()
}

View File

@ -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].

View File

@ -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<T : Comparable<T>>(
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<Double> = SplineInterpolator(DoubleField, ::DoubleBuffer)
}
}

View File

@ -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<Double> = SplineInterpolator(
DoubleField, ::DoubleBuffer
).interpolatePolynomials(data)
val polynomial: PiecewisePolynomial<Double> = 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)
}
}