forked from kscience/kmath
Minor piecewise rework
This commit is contained in:
parent
0898285542
commit
8a07140f7c
@ -10,12 +10,14 @@ import space.kscience.kmath.interpolation.interpolatePolynomials
|
|||||||
import space.kscience.kmath.operations.DoubleField
|
import space.kscience.kmath.operations.DoubleField
|
||||||
import space.kscience.kmath.structures.DoubleBuffer
|
import space.kscience.kmath.structures.DoubleBuffer
|
||||||
import space.kscience.plotly.Plotly
|
import space.kscience.plotly.Plotly
|
||||||
|
import space.kscience.plotly.UnstablePlotlyAPI
|
||||||
import space.kscience.plotly.makeFile
|
import space.kscience.plotly.makeFile
|
||||||
import space.kscience.plotly.models.functionXY
|
import space.kscience.plotly.models.functionXY
|
||||||
import space.kscience.plotly.scatter
|
import space.kscience.plotly.scatter
|
||||||
import kotlin.math.PI
|
import kotlin.math.PI
|
||||||
import kotlin.math.sin
|
import kotlin.math.sin
|
||||||
|
|
||||||
|
@OptIn(UnstablePlotlyAPI::class)
|
||||||
fun main() {
|
fun main() {
|
||||||
val data = (0..10).map {
|
val data = (0..10).map {
|
||||||
val x = it.toDouble() / 5 * PI
|
val x = it.toDouble() / 5 * PI
|
||||||
|
@ -22,8 +22,20 @@ public fun interface Piecewise<T, R> {
|
|||||||
|
|
||||||
/**
|
/**
|
||||||
* Represents piecewise-defined function where all the sub-functions are polynomials.
|
* Represents piecewise-defined function where all the sub-functions are polynomials.
|
||||||
|
* @param pieces An ordered list of range-polynomial pairs. The list does not in general guarantee that there are no "holes" in it.
|
||||||
*/
|
*/
|
||||||
public fun interface PiecewisePolynomial<T : Any> : Piecewise<T, Polynomial<T>>
|
public class PiecewisePolynomial<T : Comparable<T>>(
|
||||||
|
public val pieces: List<Pair<ClosedRange<T>, Polynomial<T>>>,
|
||||||
|
) : Piecewise<T, Polynomial<T>> {
|
||||||
|
|
||||||
|
public override fun findPiece(arg: T): Polynomial<T>? {
|
||||||
|
return if (arg < pieces.first().first.start || arg >= pieces.last().first.endInclusive)
|
||||||
|
null
|
||||||
|
else {
|
||||||
|
pieces.firstOrNull { arg in it.first }?.second
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* A [Piecewise] builder 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.
|
||||||
@ -31,7 +43,7 @@ public fun interface PiecewisePolynomial<T : Any> : Piecewise<T, Polynomial<T>>
|
|||||||
* @param T the comparable piece key type.
|
* @param T the comparable piece key type.
|
||||||
* @param delimiter the initial piecewise separator
|
* @param delimiter the initial piecewise separator
|
||||||
*/
|
*/
|
||||||
public class OrderedPiecewisePolynomial<T : Comparable<T>>(delimiter: T) : PiecewisePolynomial<T> {
|
public class PiecewiseBuilder<T : Comparable<T>>(delimiter: T) {
|
||||||
private val delimiters: MutableList<T> = arrayListOf(delimiter)
|
private val delimiters: MutableList<T> = arrayListOf(delimiter)
|
||||||
private val pieces: MutableList<Polynomial<T>> = arrayListOf()
|
private val pieces: MutableList<Polynomial<T>> = arrayListOf()
|
||||||
|
|
||||||
@ -59,17 +71,19 @@ public class OrderedPiecewisePolynomial<T : Comparable<T>>(delimiter: T) : Piece
|
|||||||
pieces.add(0, piece)
|
pieces.add(0, piece)
|
||||||
}
|
}
|
||||||
|
|
||||||
public override fun findPiece(arg: T): Polynomial<T>? {
|
public fun build(): PiecewisePolynomial<T> {
|
||||||
if (arg < delimiters.first() || arg >= delimiters.last())
|
return PiecewisePolynomial(delimiters.zipWithNext { l, r -> l..r }.zip(pieces))
|
||||||
return null
|
|
||||||
else {
|
|
||||||
for (index in 1 until delimiters.size)
|
|
||||||
if (arg < delimiters[index]) return pieces[index - 1]
|
|
||||||
error("Piece not found")
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* A builder for [PiecewisePolynomial]
|
||||||
|
*/
|
||||||
|
public fun <T : Comparable<T>> PiecewisePolynomial(
|
||||||
|
startingPoint: T,
|
||||||
|
builder: PiecewiseBuilder<T>.() -> Unit,
|
||||||
|
): PiecewisePolynomial<T> = PiecewiseBuilder(startingPoint).apply(builder).build()
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise
|
* Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise
|
||||||
* definition.
|
* definition.
|
||||||
|
@ -5,10 +5,8 @@
|
|||||||
|
|
||||||
package space.kscience.kmath.functions
|
package space.kscience.kmath.functions
|
||||||
|
|
||||||
import space.kscience.kmath.operations.Group
|
import space.kscience.kmath.misc.UnstableKMathAPI
|
||||||
import space.kscience.kmath.operations.Ring
|
import space.kscience.kmath.operations.*
|
||||||
import space.kscience.kmath.operations.ScaleOperations
|
|
||||||
import space.kscience.kmath.operations.invoke
|
|
||||||
import kotlin.contracts.InvocationKind
|
import kotlin.contracts.InvocationKind
|
||||||
import kotlin.contracts.contract
|
import kotlin.contracts.contract
|
||||||
import kotlin.math.max
|
import kotlin.math.max
|
||||||
@ -19,7 +17,7 @@ import kotlin.math.pow
|
|||||||
*
|
*
|
||||||
* @param coefficients constant is the leftmost coefficient.
|
* @param coefficients constant is the leftmost coefficient.
|
||||||
*/
|
*/
|
||||||
public class Polynomial<T>(public val coefficients: List<T>){
|
public class Polynomial<T>(public val coefficients: List<T>) {
|
||||||
override fun toString(): String = "Polynomial$coefficients"
|
override fun toString(): String = "Polynomial$coefficients"
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -32,7 +30,9 @@ public fun <T> Polynomial(vararg coefficients: T): Polynomial<T> = Polynomial(co
|
|||||||
/**
|
/**
|
||||||
* Evaluates the value of the given double polynomial for given double argument.
|
* Evaluates the value of the given double polynomial for given double argument.
|
||||||
*/
|
*/
|
||||||
public fun Polynomial<Double>.value(): Double = coefficients.reduceIndexed { index, acc, d -> acc + d.pow(index) }
|
public fun Polynomial<Double>.value(arg: Double): Double = coefficients.reduceIndexed { index, acc, c ->
|
||||||
|
acc + c * arg.pow(index)
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Evaluates the value of the given polynomial for given argument.
|
* Evaluates the value of the given polynomial for given argument.
|
||||||
@ -50,9 +50,38 @@ public fun <T, C : Ring<T>> Polynomial<T>.value(ring: C, arg: T): T = ring {
|
|||||||
/**
|
/**
|
||||||
* Represent the polynomial as a regular context-less function.
|
* Represent the polynomial as a regular context-less function.
|
||||||
*/
|
*/
|
||||||
public fun <T , C : Ring<T>> Polynomial<T>.asFunction(ring: C): (T) -> T = { value(ring, it) }
|
public fun <T, C : Ring<T>> Polynomial<T>.asFunction(ring: C): (T) -> T = { value(ring, it) }
|
||||||
|
|
||||||
//public fun <T: Any>
|
/**
|
||||||
|
* Create a polynomial witch represents differentiated version of this polynomial
|
||||||
|
*/
|
||||||
|
@UnstableKMathAPI
|
||||||
|
public fun <T, A> Polynomial<T>.differentiate(
|
||||||
|
algebra: A,
|
||||||
|
): Polynomial<T> where A : Ring<T>, A : NumericAlgebra<T> = algebra {
|
||||||
|
Polynomial(coefficients.drop(1).mapIndexed { index, t -> number(index) * t })
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Create a polynomial witch represents indefinite integral version of this polynomial
|
||||||
|
*/
|
||||||
|
@UnstableKMathAPI
|
||||||
|
public fun <T, A> Polynomial<T>.integrate(
|
||||||
|
algebra: A,
|
||||||
|
): Polynomial<T> where A : Field<T>, A : NumericAlgebra<T> = algebra {
|
||||||
|
Polynomial(coefficients.mapIndexed { index, t -> t / number(index) })
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Compute a definite integral of a given polynomial in a [range]
|
||||||
|
*/
|
||||||
|
@UnstableKMathAPI
|
||||||
|
public fun <T : Comparable<T>, A> Polynomial<T>.integrate(
|
||||||
|
algebra: A,
|
||||||
|
range: ClosedRange<T>,
|
||||||
|
): T where A : Field<T>, A : NumericAlgebra<T> = algebra {
|
||||||
|
value(algebra, range.endInclusive) - value(algebra, range.start)
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Space of polynomials.
|
* Space of polynomials.
|
||||||
@ -87,6 +116,9 @@ public class PolynomialSpace<T, C>(
|
|||||||
* Evaluates the polynomial for the given value [arg].
|
* Evaluates the polynomial for the given value [arg].
|
||||||
*/
|
*/
|
||||||
public operator fun Polynomial<T>.invoke(arg: T): T = value(ring, arg)
|
public operator fun Polynomial<T>.invoke(arg: T): T = value(ring, arg)
|
||||||
|
|
||||||
|
public fun Polynomial<T>.asFunction(): (T) -> T = asFunction(ring)
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
public inline fun <T, C, R> C.polynomial(block: PolynomialSpace<T, C>.() -> R): R where C : Ring<T>, C : ScaleOperations<T> {
|
public inline fun <T, C, R> C.polynomial(block: PolynomialSpace<T, C>.() -> R): R where C : Ring<T>, C : ScaleOperations<T> {
|
||||||
|
@ -6,7 +6,6 @@
|
|||||||
package space.kscience.kmath.interpolation
|
package space.kscience.kmath.interpolation
|
||||||
|
|
||||||
import space.kscience.kmath.data.XYColumnarData
|
import space.kscience.kmath.data.XYColumnarData
|
||||||
import space.kscience.kmath.functions.OrderedPiecewisePolynomial
|
|
||||||
import space.kscience.kmath.functions.PiecewisePolynomial
|
import space.kscience.kmath.functions.PiecewisePolynomial
|
||||||
import space.kscience.kmath.functions.Polynomial
|
import space.kscience.kmath.functions.Polynomial
|
||||||
import space.kscience.kmath.misc.UnstableKMathAPI
|
import space.kscience.kmath.misc.UnstableKMathAPI
|
||||||
@ -28,7 +27,7 @@ public class LinearInterpolator<T : Comparable<T>>(public override val algebra:
|
|||||||
require(points.size > 0) { "Point array should not be empty" }
|
require(points.size > 0) { "Point array should not be empty" }
|
||||||
insureSorted(points)
|
insureSorted(points)
|
||||||
|
|
||||||
OrderedPiecewisePolynomial(points.x[0]).apply {
|
PiecewisePolynomial(points.x[0]) {
|
||||||
for (i in 0 until points.size - 1) {
|
for (i in 0 until points.size - 1) {
|
||||||
val slope = (points.y[i + 1] - points.y[i]) / (points.x[i + 1] - points.x[i])
|
val slope = (points.y[i + 1] - points.y[i]) / (points.x[i + 1] - points.x[i])
|
||||||
val const = points.y[i] - slope * points.x[i]
|
val const = points.y[i] - slope * points.x[i]
|
||||||
|
@ -6,7 +6,6 @@
|
|||||||
package space.kscience.kmath.interpolation
|
package space.kscience.kmath.interpolation
|
||||||
|
|
||||||
import space.kscience.kmath.data.XYColumnarData
|
import space.kscience.kmath.data.XYColumnarData
|
||||||
import space.kscience.kmath.functions.OrderedPiecewisePolynomial
|
|
||||||
import space.kscience.kmath.functions.PiecewisePolynomial
|
import space.kscience.kmath.functions.PiecewisePolynomial
|
||||||
import space.kscience.kmath.functions.Polynomial
|
import space.kscience.kmath.functions.Polynomial
|
||||||
import space.kscience.kmath.misc.UnstableKMathAPI
|
import space.kscience.kmath.misc.UnstableKMathAPI
|
||||||
@ -50,7 +49,7 @@ public class SplineInterpolator<T : Comparable<T>>(
|
|||||||
|
|
||||||
// cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)
|
// cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)
|
||||||
|
|
||||||
OrderedPiecewisePolynomial(points.x[points.size - 1]).apply {
|
PiecewisePolynomial(points.x[points.size - 1]) {
|
||||||
var cOld = zero
|
var cOld = zero
|
||||||
|
|
||||||
for (j in n - 1 downTo 0) {
|
for (j in n - 1 downTo 0) {
|
||||||
|
@ -0,0 +1,17 @@
|
|||||||
|
/*
|
||||||
|
* 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 kotlin.test.Test
|
||||||
|
import kotlin.test.assertEquals
|
||||||
|
|
||||||
|
class PolynomialTest {
|
||||||
|
@Test
|
||||||
|
fun testIntegration() {
|
||||||
|
val polynomial = Polynomial(1.0, -2.0, 1.0)
|
||||||
|
assertEquals(0.0, polynomial.value(1.0), 0.001)
|
||||||
|
}
|
||||||
|
}
|
@ -28,7 +28,6 @@ private fun <B : ClosedFloatingPointRange<Double>> TreeMap<Double, B>.getBin(val
|
|||||||
|
|
||||||
@UnstableKMathAPI
|
@UnstableKMathAPI
|
||||||
public class TreeHistogram(
|
public class TreeHistogram(
|
||||||
override val context: TreeHistogramSpace,
|
|
||||||
private val binMap: TreeMap<Double, out UnivariateBin>,
|
private val binMap: TreeMap<Double, out UnivariateBin>,
|
||||||
) : UnivariateHistogram {
|
) : UnivariateHistogram {
|
||||||
override fun get(value: Double): UnivariateBin? = binMap.getBin(value)
|
override fun get(value: Double): UnivariateBin? = binMap.getBin(value)
|
||||||
@ -79,15 +78,15 @@ public class TreeHistogramSpace(
|
|||||||
val count = binCounter.counter.value
|
val count = binCounter.counter.value
|
||||||
resBins[key] = UnivariateBin(binCounter.domain, count, sqrt(count))
|
resBins[key] = UnivariateBin(binCounter.domain, count, sqrt(count))
|
||||||
}
|
}
|
||||||
return TreeHistogram(this, resBins)
|
return TreeHistogram(resBins)
|
||||||
}
|
}
|
||||||
|
|
||||||
override fun add(
|
override fun add(
|
||||||
a: UnivariateHistogram,
|
a: UnivariateHistogram,
|
||||||
b: UnivariateHistogram,
|
b: UnivariateHistogram,
|
||||||
): UnivariateHistogram {
|
): UnivariateHistogram {
|
||||||
require(a.context == this) { "Histogram $a does not belong to this context" }
|
// require(a.context == this) { "Histogram $a does not belong to this context" }
|
||||||
require(b.context == this) { "Histogram $b does not belong to this context" }
|
// require(b.context == this) { "Histogram $b does not belong to this context" }
|
||||||
val bins = TreeMap<Double, UnivariateBin>().apply {
|
val bins = TreeMap<Double, UnivariateBin>().apply {
|
||||||
(a.bins.map { it.domain } union b.bins.map { it.domain }).forEach { def ->
|
(a.bins.map { it.domain } union b.bins.map { it.domain }).forEach { def ->
|
||||||
put(def.center,
|
put(def.center,
|
||||||
@ -100,7 +99,7 @@ public class TreeHistogramSpace(
|
|||||||
)
|
)
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
return TreeHistogram(this, bins)
|
return TreeHistogram(bins)
|
||||||
}
|
}
|
||||||
|
|
||||||
override fun scale(a: UnivariateHistogram, value: Double): UnivariateHistogram {
|
override fun scale(a: UnivariateHistogram, value: Double): UnivariateHistogram {
|
||||||
@ -116,7 +115,7 @@ public class TreeHistogramSpace(
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
return TreeHistogram(this, bins)
|
return TreeHistogram(bins)
|
||||||
}
|
}
|
||||||
|
|
||||||
override fun UnivariateHistogram.unaryMinus(): UnivariateHistogram = this * (-1)
|
override fun UnivariateHistogram.unaryMinus(): UnivariateHistogram = this * (-1)
|
||||||
|
@ -7,8 +7,6 @@ package space.kscience.kmath.histogram
|
|||||||
|
|
||||||
import space.kscience.kmath.domains.UnivariateDomain
|
import space.kscience.kmath.domains.UnivariateDomain
|
||||||
import space.kscience.kmath.misc.UnstableKMathAPI
|
import space.kscience.kmath.misc.UnstableKMathAPI
|
||||||
import space.kscience.kmath.operations.Group
|
|
||||||
import space.kscience.kmath.operations.GroupElement
|
|
||||||
import space.kscience.kmath.structures.Buffer
|
import space.kscience.kmath.structures.Buffer
|
||||||
import space.kscience.kmath.structures.asSequence
|
import space.kscience.kmath.structures.asSequence
|
||||||
|
|
||||||
@ -35,8 +33,7 @@ public class UnivariateBin(
|
|||||||
}
|
}
|
||||||
|
|
||||||
@OptIn(UnstableKMathAPI::class)
|
@OptIn(UnstableKMathAPI::class)
|
||||||
public interface UnivariateHistogram : Histogram<Double, UnivariateBin>,
|
public interface UnivariateHistogram : Histogram<Double, UnivariateBin>{
|
||||||
GroupElement<UnivariateHistogram, Group<UnivariateHistogram>> {
|
|
||||||
public operator fun get(value: Double): UnivariateBin?
|
public operator fun get(value: Double): UnivariateBin?
|
||||||
public override operator fun get(point: Buffer<Double>): UnivariateBin? = get(point[0])
|
public override operator fun get(point: Buffer<Double>): UnivariateBin? = get(point[0])
|
||||||
|
|
||||||
|
Loading…
Reference in New Issue
Block a user