kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

Merge pull request #250 from mipt-npm/feature/integration

Browse Source 
Feature/integration
This commit is contained in:
Alexander Nozik 2021-03-16 20:44:41 +03:00 committed by GitHub
commit 7c0d805018
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 341 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt Normal file
View File

@ -0,0 +1,92 @@
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.*
import space.kscience.kmath.misc.UnstableKMathAPI
/**
* 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 integrate(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
)
}
}
}
@UnstableKMathAPI
public var MutableList<IntegrandFeature>.targetAbsoluteAccuracy: Double?
get() = filterIsInstance<CMIntegrator.TargetAbsoluteAccuracy>().lastOrNull()?.value
set(value){
value?.let { add(CMIntegrator.TargetAbsoluteAccuracy(value))}
}
@UnstableKMathAPI
public var MutableList<IntegrandFeature>.targetRelativeAccuracy: Double?
get() = filterIsInstance<CMIntegrator.TargetRelativeAccuracy>().lastOrNull()?.value
set(value){
value?.let { add(CMIntegrator.TargetRelativeAccuracy(value))}
}

94
kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/GaussRuleIntegrator.kt Normal file
View File

@ -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 integrate(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).integrate(
UnivariateIntegrand(function, IntegrationRange(range))
).value!!
}
}

2
kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimizationProblem.kt
View File

@ -19,7 +19,7 @@ import kotlin.reflect.KClass
public operator fun PointValuePair.component1(): DoubleArray = point
public operator fun PointValuePair.component2(): Double = value
public class CMOptimizationProblem(override val symbols: List<Symbol>, ) :
public class CMOptimizationProblem(override val symbols: List<Symbol>) :
OptimizationProblem<Double>, SymbolIndexer, OptimizationFeature {
private val optimizationData: HashMap<KClass<out OptimizationData>, OptimizationData> = HashMap()
private var optimizatorBuilder: (() -> MultivariateOptimizer)? = null

30
kmath-commons/src/test/kotlin/space/kscience/kmath/commons/integration/IntegrationTest.kt Normal file
View File

@ -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) {
targetRelativeAccuracy = 1e-4
targetAbsoluteAccuracy = 1e-4
}
assertTrue { abs(res - 2) < 1e-3 }
assertTrue { abs(res - 2) > 1e-12 }
}
}

22
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrand.kt Normal file
View File

@ -0,0 +1,22 @@
package space.kscience.kmath.integration
import kotlin.reflect.KClass
public interface IntegrandFeature
public interface Integrand {
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

11
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt Normal file
View File

@ -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 integrate(integrand: I): I
}

27
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/MultivariateIntegrand.kt Normal file
View File

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

64
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt Normal file
View File

@ -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 = integrate(
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 integrate(
UnivariateIntegrand(function, *features.toTypedArray())
).value ?: error("Unexpected: no value after integration.")
}