API for integration
This commit is contained in:
parent
b2df860381
commit
874875c13a
@ -0,0 +1,77 @@
|
|||||||
|
package space.kscience.kmath.commons.integration
|
||||||
|
|
||||||
|
import org.apache.commons.math3.analysis.integration.IterativeLegendreGaussIntegrator
|
||||||
|
import org.apache.commons.math3.analysis.integration.SimpsonIntegrator
|
||||||
|
import space.kscience.kmath.integration.*
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Integration wrapper for Common-maths UnivariateIntegrator
|
||||||
|
*/
|
||||||
|
public class CMIntegrator(
|
||||||
|
private val defaultMaxCalls: Int = 200,
|
||||||
|
public val integratorBuilder: (Integrand) -> org.apache.commons.math3.analysis.integration.UnivariateIntegrator,
|
||||||
|
) : UnivariateIntegrator<Double> {
|
||||||
|
|
||||||
|
public class TargetRelativeAccuracy(public val value: Double) : IntegrandFeature
|
||||||
|
public class TargetAbsoluteAccuracy(public val value: Double) : IntegrandFeature
|
||||||
|
|
||||||
|
public class MinIterations(public val value: Int) : IntegrandFeature
|
||||||
|
public class MaxIterations(public val value: Int) : IntegrandFeature
|
||||||
|
|
||||||
|
override fun evaluate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
|
||||||
|
val integrator = integratorBuilder(integrand)
|
||||||
|
val maxCalls = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: defaultMaxCalls
|
||||||
|
val remainingCalls = maxCalls - integrand.calls
|
||||||
|
val range = integrand.getFeature<IntegrationRange<Double>>()?.range
|
||||||
|
?: error("Integration range is not provided")
|
||||||
|
val res = integrator.integrate(remainingCalls, integrand.function, range.start, range.endInclusive)
|
||||||
|
|
||||||
|
return integrand +
|
||||||
|
IntegrandValue(res) +
|
||||||
|
IntegrandAbsoluteAccuracy(integrator.absoluteAccuracy) +
|
||||||
|
IntegrandRelativeAccuracy(integrator.relativeAccuracy) +
|
||||||
|
IntegrandCalls(integrator.evaluations + integrand.calls)
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
public companion object {
|
||||||
|
/**
|
||||||
|
* Create a Simpson integrator based on [SimpsonIntegrator]
|
||||||
|
*/
|
||||||
|
public fun simpson(defaultMaxCalls: Int = 200): CMIntegrator = CMIntegrator(defaultMaxCalls) { integrand ->
|
||||||
|
val absoluteAccuracy = integrand.getFeature<TargetAbsoluteAccuracy>()?.value
|
||||||
|
?: SimpsonIntegrator.DEFAULT_ABSOLUTE_ACCURACY
|
||||||
|
val relativeAccuracy = integrand.getFeature<TargetRelativeAccuracy>()?.value
|
||||||
|
?: SimpsonIntegrator.DEFAULT_ABSOLUTE_ACCURACY
|
||||||
|
val minIterations = integrand.getFeature<MinIterations>()?.value
|
||||||
|
?: SimpsonIntegrator.DEFAULT_MIN_ITERATIONS_COUNT
|
||||||
|
val maxIterations = integrand.getFeature<MaxIterations>()?.value
|
||||||
|
?: SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT
|
||||||
|
|
||||||
|
SimpsonIntegrator(relativeAccuracy, absoluteAccuracy, minIterations, maxIterations)
|
||||||
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Create a Gauss-Legandre integrator based on [IterativeLegendreGaussIntegrator]
|
||||||
|
*/
|
||||||
|
public fun legandre(numPoints: Int, defaultMaxCalls: Int = numPoints * 5): CMIntegrator =
|
||||||
|
CMIntegrator(defaultMaxCalls) { integrand ->
|
||||||
|
val absoluteAccuracy = integrand.getFeature<TargetAbsoluteAccuracy>()?.value
|
||||||
|
?: IterativeLegendreGaussIntegrator.DEFAULT_ABSOLUTE_ACCURACY
|
||||||
|
val relativeAccuracy = integrand.getFeature<TargetRelativeAccuracy>()?.value
|
||||||
|
?: IterativeLegendreGaussIntegrator.DEFAULT_ABSOLUTE_ACCURACY
|
||||||
|
val minIterations = integrand.getFeature<MinIterations>()?.value
|
||||||
|
?: IterativeLegendreGaussIntegrator.DEFAULT_MIN_ITERATIONS_COUNT
|
||||||
|
val maxIterations = integrand.getFeature<MaxIterations>()?.value
|
||||||
|
?: IterativeLegendreGaussIntegrator.DEFAULT_MAX_ITERATIONS_COUNT
|
||||||
|
|
||||||
|
IterativeLegendreGaussIntegrator(
|
||||||
|
numPoints,
|
||||||
|
relativeAccuracy,
|
||||||
|
absoluteAccuracy,
|
||||||
|
minIterations,
|
||||||
|
maxIterations
|
||||||
|
)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
@ -0,0 +1,94 @@
|
|||||||
|
/*
|
||||||
|
* Copyright 2015 Alexander Nozik.
|
||||||
|
*
|
||||||
|
* Licensed under the Apache License, Version 2.0 (the "License");
|
||||||
|
* you may not use this file except in compliance with the License.
|
||||||
|
* You may obtain a copy of the License at
|
||||||
|
*
|
||||||
|
* http://www.apache.org/licenses/LICENSE-2.0
|
||||||
|
*
|
||||||
|
* Unless required by applicable law or agreed to in writing, software
|
||||||
|
* distributed under the License is distributed on an "AS IS" BASIS,
|
||||||
|
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||||
|
* See the License for the specific language governing permissions and
|
||||||
|
* limitations under the License.
|
||||||
|
*/
|
||||||
|
package space.kscience.kmath.commons.integration
|
||||||
|
|
||||||
|
import org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator
|
||||||
|
import org.apache.commons.math3.analysis.integration.gauss.GaussIntegratorFactory
|
||||||
|
import space.kscience.kmath.integration.*
|
||||||
|
|
||||||
|
/**
|
||||||
|
* A simple one-pass integrator based on Gauss rule
|
||||||
|
*/
|
||||||
|
public class GaussRuleIntegrator(
|
||||||
|
private val numpoints: Int,
|
||||||
|
private var type: GaussRule = GaussRule.LEGANDRE,
|
||||||
|
) : UnivariateIntegrator<Double> {
|
||||||
|
|
||||||
|
override fun evaluate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
|
||||||
|
val range = integrand.getFeature<IntegrationRange<Double>>()?.range
|
||||||
|
?: error("Integration range is not provided")
|
||||||
|
val integrator: GaussIntegrator = getIntegrator(range)
|
||||||
|
//TODO check performance
|
||||||
|
val res: Double = integrator.integrate(integrand.function)
|
||||||
|
return integrand + IntegrandValue(res) + IntegrandCalls(integrand.calls + numpoints)
|
||||||
|
}
|
||||||
|
|
||||||
|
private fun getIntegrator(range: ClosedRange<Double>): GaussIntegrator {
|
||||||
|
return when (type) {
|
||||||
|
GaussRule.LEGANDRE -> factory.legendre(
|
||||||
|
numpoints,
|
||||||
|
range.start,
|
||||||
|
range.endInclusive
|
||||||
|
)
|
||||||
|
GaussRule.LEGANDREHP -> factory.legendreHighPrecision(
|
||||||
|
numpoints,
|
||||||
|
range.start,
|
||||||
|
range.endInclusive
|
||||||
|
)
|
||||||
|
GaussRule.UNIFORM -> GaussIntegrator(
|
||||||
|
getUniformRule(
|
||||||
|
range.start,
|
||||||
|
range.endInclusive,
|
||||||
|
numpoints
|
||||||
|
)
|
||||||
|
)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
private fun getUniformRule(
|
||||||
|
min: Double,
|
||||||
|
max: Double,
|
||||||
|
numPoints: Int,
|
||||||
|
): org.apache.commons.math3.util.Pair<DoubleArray, DoubleArray> {
|
||||||
|
assert(numPoints > 2)
|
||||||
|
val points = DoubleArray(numPoints)
|
||||||
|
val weights = DoubleArray(numPoints)
|
||||||
|
val step = (max - min) / (numPoints - 1)
|
||||||
|
points[0] = min
|
||||||
|
for (i in 1 until numPoints) {
|
||||||
|
points[i] = points[i - 1] + step
|
||||||
|
weights[i] = step
|
||||||
|
}
|
||||||
|
return org.apache.commons.math3.util.Pair<DoubleArray, DoubleArray>(points, weights)
|
||||||
|
}
|
||||||
|
|
||||||
|
public enum class GaussRule {
|
||||||
|
UNIFORM, LEGANDRE, LEGANDREHP
|
||||||
|
}
|
||||||
|
|
||||||
|
public companion object {
|
||||||
|
private val factory: GaussIntegratorFactory = GaussIntegratorFactory()
|
||||||
|
|
||||||
|
public fun integrate(
|
||||||
|
range: ClosedRange<Double>,
|
||||||
|
numPoints: Int = 100,
|
||||||
|
type: GaussRule = GaussRule.LEGANDRE,
|
||||||
|
function: (Double) -> Double,
|
||||||
|
): Double = GaussRuleIntegrator(numPoints, type).evaluate(
|
||||||
|
UnivariateIntegrand(function, IntegrationRange(range))
|
||||||
|
).value!!
|
||||||
|
}
|
||||||
|
}
|
@ -0,0 +1,30 @@
|
|||||||
|
package space.kscience.kmath.commons.integration
|
||||||
|
|
||||||
|
import org.junit.jupiter.api.Test
|
||||||
|
import space.kscience.kmath.integration.integrate
|
||||||
|
import space.kscience.kmath.misc.UnstableKMathAPI
|
||||||
|
import space.kscience.kmath.operations.RealField.sin
|
||||||
|
import kotlin.math.PI
|
||||||
|
import kotlin.math.abs
|
||||||
|
import kotlin.test.assertTrue
|
||||||
|
|
||||||
|
@UnstableKMathAPI
|
||||||
|
internal class IntegrationTest {
|
||||||
|
private val function: (Double) -> Double = { sin(it) }
|
||||||
|
|
||||||
|
@Test
|
||||||
|
fun simpson() {
|
||||||
|
val res = CMIntegrator.simpson().integrate(0.0..2 * PI, function)
|
||||||
|
assertTrue { abs(res) < 1e-3 }
|
||||||
|
}
|
||||||
|
|
||||||
|
@Test
|
||||||
|
fun customSimpson() {
|
||||||
|
val res = CMIntegrator.simpson().integrate(0.0..PI, function) {
|
||||||
|
add(CMIntegrator.TargetRelativeAccuracy(1e-4))
|
||||||
|
add(CMIntegrator.TargetAbsoluteAccuracy(1e-4))
|
||||||
|
}
|
||||||
|
assertTrue { abs(res - 2) < 1e-3 }
|
||||||
|
assertTrue { abs(res - 2) > 1e-12 }
|
||||||
|
}
|
||||||
|
}
|
@ -5,5 +5,18 @@ import kotlin.reflect.KClass
|
|||||||
public interface IntegrandFeature
|
public interface IntegrandFeature
|
||||||
|
|
||||||
public interface Integrand {
|
public interface Integrand {
|
||||||
public fun <T : IntegrandFeature> getFeature(type: KClass<T>): T? = null
|
public fun <T : IntegrandFeature> getFeature(type: KClass<T>): T?
|
||||||
}
|
}
|
||||||
|
|
||||||
|
public inline fun <reified T : IntegrandFeature> Integrand.getFeature(): T? = getFeature(T::class)
|
||||||
|
|
||||||
|
public class IntegrandValue<T : Any>(public val value: T) : IntegrandFeature
|
||||||
|
|
||||||
|
public class IntegrandRelativeAccuracy(public val accuracy: Double) : IntegrandFeature
|
||||||
|
|
||||||
|
public class IntegrandAbsoluteAccuracy(public val accuracy: Double) : IntegrandFeature
|
||||||
|
|
||||||
|
public class IntegrandCalls(public val calls: Int) : IntegrandFeature
|
||||||
|
public val Integrand.calls: Int get() = getFeature<IntegrandCalls>()?.calls ?: 0
|
||||||
|
|
||||||
|
public class IntegrandMaxCalls(public val maxCalls: Int) : IntegrandFeature
|
||||||
|
@ -0,0 +1,11 @@
|
|||||||
|
package space.kscience.kmath.integration
|
||||||
|
|
||||||
|
/**
|
||||||
|
* A general interface for all integrators
|
||||||
|
*/
|
||||||
|
public interface Integrator<I: Integrand> {
|
||||||
|
/**
|
||||||
|
* Run one integration pass and return a new [Integrand] with a new set of features
|
||||||
|
*/
|
||||||
|
public fun evaluate(integrand: I): I
|
||||||
|
}
|
@ -0,0 +1,27 @@
|
|||||||
|
package space.kscience.kmath.integration
|
||||||
|
|
||||||
|
import space.kscience.kmath.linear.Point
|
||||||
|
import kotlin.reflect.KClass
|
||||||
|
|
||||||
|
public class MultivariateIntegrand<T : Any> internal constructor(
|
||||||
|
private val features: Map<KClass<*>, IntegrandFeature>,
|
||||||
|
public val function: (Point<T>) -> T,
|
||||||
|
) : Integrand {
|
||||||
|
|
||||||
|
@Suppress("UNCHECKED_CAST")
|
||||||
|
override fun <T : IntegrandFeature> getFeature(type: KClass<T>): T? = features[type] as? T
|
||||||
|
|
||||||
|
public operator fun <F : IntegrandFeature> plus(pair: Pair<KClass<out F>, F>): MultivariateIntegrand<T> =
|
||||||
|
MultivariateIntegrand(features + pair, function)
|
||||||
|
|
||||||
|
public operator fun <F : IntegrandFeature> plus(feature: F): MultivariateIntegrand<T> =
|
||||||
|
plus(feature::class to feature)
|
||||||
|
}
|
||||||
|
|
||||||
|
@Suppress("FunctionName")
|
||||||
|
public fun <T : Any> MultivariateIntegrand(
|
||||||
|
vararg features: IntegrandFeature,
|
||||||
|
function: (Point<T>) -> T,
|
||||||
|
): MultivariateIntegrand<T> = MultivariateIntegrand(features.associateBy { it::class }, function)
|
||||||
|
|
||||||
|
public val <T : Any> MultivariateIntegrand<T>.value: T? get() = getFeature<IntegrandValue<T>>()?.value
|
@ -0,0 +1,64 @@
|
|||||||
|
package space.kscience.kmath.integration
|
||||||
|
|
||||||
|
import space.kscience.kmath.misc.UnstableKMathAPI
|
||||||
|
import kotlin.reflect.KClass
|
||||||
|
|
||||||
|
public class UnivariateIntegrand<T : Any> internal constructor(
|
||||||
|
private val features: Map<KClass<*>, IntegrandFeature>,
|
||||||
|
public val function: (T) -> T,
|
||||||
|
) : Integrand {
|
||||||
|
|
||||||
|
@Suppress("UNCHECKED_CAST")
|
||||||
|
override fun <T : IntegrandFeature> getFeature(type: KClass<T>): T? = features[type] as? T
|
||||||
|
|
||||||
|
public operator fun <F : IntegrandFeature> plus(pair: Pair<KClass<out F>, F>): UnivariateIntegrand<T> =
|
||||||
|
UnivariateIntegrand(features + pair, function)
|
||||||
|
|
||||||
|
public operator fun <F : IntegrandFeature> plus(feature: F): UnivariateIntegrand<T> =
|
||||||
|
plus(feature::class to feature)
|
||||||
|
}
|
||||||
|
|
||||||
|
@Suppress("FunctionName")
|
||||||
|
public fun <T : Any> UnivariateIntegrand(
|
||||||
|
function: (T) -> T,
|
||||||
|
vararg features: IntegrandFeature,
|
||||||
|
): UnivariateIntegrand<T> = UnivariateIntegrand(features.associateBy { it::class }, function)
|
||||||
|
|
||||||
|
public typealias UnivariateIntegrator<T> = Integrator<UnivariateIntegrand<T>>
|
||||||
|
|
||||||
|
public inline class IntegrationRange<T : Comparable<T>>(public val range: ClosedRange<T>) : IntegrandFeature
|
||||||
|
|
||||||
|
public val <T : Any> UnivariateIntegrand<T>.value: T? get() = getFeature<IntegrandValue<T>>()?.value
|
||||||
|
|
||||||
|
/**
|
||||||
|
* A shortcut method to integrate a [function] in [range] with additional [features].
|
||||||
|
* The [function] is placed in the end position to allow passing a lambda.
|
||||||
|
*/
|
||||||
|
@UnstableKMathAPI
|
||||||
|
public fun UnivariateIntegrator<Double>.integrate(
|
||||||
|
range: ClosedRange<Double>,
|
||||||
|
vararg features: IntegrandFeature,
|
||||||
|
function: (Double) -> Double,
|
||||||
|
): Double = evaluate(
|
||||||
|
UnivariateIntegrand(function, IntegrationRange(range), *features)
|
||||||
|
).value ?: error("Unexpected: no value after integration.")
|
||||||
|
|
||||||
|
/**
|
||||||
|
* A shortcut method to integrate a [function] in [range] with additional [features].
|
||||||
|
* The [function] is placed in the end position to allow passing a lambda.
|
||||||
|
*/
|
||||||
|
@UnstableKMathAPI
|
||||||
|
public fun UnivariateIntegrator<Double>.integrate(
|
||||||
|
range: ClosedRange<Double>,
|
||||||
|
function: (Double) -> Double,
|
||||||
|
featureBuilder: (MutableList<IntegrandFeature>.() -> Unit) = {},
|
||||||
|
): Double {
|
||||||
|
//TODO use dedicated feature builder class instead or add extensions to MutableList<IntegrandFeature>
|
||||||
|
val features = buildList {
|
||||||
|
featureBuilder()
|
||||||
|
add(IntegrationRange(range))
|
||||||
|
}
|
||||||
|
return evaluate(
|
||||||
|
UnivariateIntegrand(function, *features.toTypedArray())
|
||||||
|
).value ?: error("Unexpected: no value after integration.")
|
||||||
|
}
|
Loading…
Reference in New Issue
Block a user