teldufalsari/kmath
2
0
Fork 0
forked from kscience/kmath
Code Pull Requests 1 Activity

Update integration to use Attributes

Browse Source
This commit is contained in:
Alexander Nozik 2023-08-13 19:13:39 +03:00
parent eff70eb690
commit 5196322b7a
19 changed files with 173 additions and 160 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attributes.kt
View File

@ -59,6 +59,14 @@ public fun <T : Any, A : Attribute<T>> Attributes.withAttribute(
public fun <A : Attribute<Unit>> Attributes.withAttribute(attribute: A): Attributes =
withAttribute(attribute, Unit)
/**
* Create a new [Attributes] by modifying the current one
*/
public fun Attributes.modify(block: AttributesBuilder.() -> Unit): Attributes = Attributes {
from(this@modify)
block()
}
/**
* Create new [Attributes] by removing [attribute] key
*/

7
attributes-kt/src/commonMain/kotlin/space/kscience/attributes/AttributesBuilder.kt
View File

@ -8,7 +8,9 @@ package space.kscience.attributes
/**
* A safe builder for [Attributes]
*/
public class AttributesBuilder internal constructor(private val map: MutableMap<Attribute<*>, Any> = mutableMapOf()) {
public class AttributesBuilder internal constructor(private val map: MutableMap<Attribute<*>, Any>) {
public constructor() : this(mutableMapOf())
@Suppress("UNCHECKED_CAST")
public operator fun <T> get(attribute: Attribute<T>): T? = map[attribute] as? T
@ -49,4 +51,5 @@ public fun AttributesBuilder(
attributes: Attributes,
): AttributesBuilder = AttributesBuilder(attributes.content.toMutableMap())
public fun Attributes(builder: AttributesBuilder.() -> Unit): Attributes = AttributesBuilder().apply(builder).build()
public inline fun Attributes(builder: AttributesBuilder.() -> Unit): Attributes =
AttributesBuilder().apply(builder).build()

19
kmath-commons/src/main/kotlin/space/kscience/kmath/commons/integration/CMIntegrator.kt
View File

@ -20,10 +20,9 @@ public class CMIntegrator(
override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
val integrator = integratorBuilder(integrand)
val maxCalls = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: defaultMaxCalls
val maxCalls = integrand[IntegrandMaxCalls] ?: defaultMaxCalls
val remainingCalls = maxCalls - integrand.calls
val range = integrand.getFeature<IntegrationRange>()?.range
?: error("Integration range is not provided")
val range = integrand[IntegrationRange] ?: error("Integration range is not provided")
val res = integrator.integrate(remainingCalls, integrand.function, range.start, range.endInclusive)
return integrand +
@ -39,11 +38,9 @@ public class CMIntegrator(
* Create a Simpson integrator based on [SimpsonIntegrator]
*/
public fun simpson(defaultMaxCalls: Int = 200): CMIntegrator = CMIntegrator(defaultMaxCalls) { integrand ->
val absoluteAccuracy = integrand.getFeature<IntegrandAbsoluteAccuracy>()?.accuracy
?: SimpsonIntegrator.DEFAULT_ABSOLUTE_ACCURACY
val relativeAccuracy = integrand.getFeature<IntegrandRelativeAccuracy>()?.accuracy
?: SimpsonIntegrator.DEFAULT_ABSOLUTE_ACCURACY
val iterations = integrand.getFeature<IntegrandIterationsRange>()?.range
val absoluteAccuracy = integrand[IntegrandAbsoluteAccuracy] ?: SimpsonIntegrator.DEFAULT_ABSOLUTE_ACCURACY
val relativeAccuracy = integrand[IntegrandRelativeAccuracy] ?: SimpsonIntegrator.DEFAULT_ABSOLUTE_ACCURACY
val iterations = integrand[IntegrandIterationsRange]
?: SimpsonIntegrator.DEFAULT_MIN_ITERATIONS_COUNT..SimpsonIntegrator.SIMPSON_MAX_ITERATIONS_COUNT
@ -55,11 +52,11 @@ public class CMIntegrator(
*/
public fun legandre(numPoints: Int, defaultMaxCalls: Int = numPoints * 5): CMIntegrator =
CMIntegrator(defaultMaxCalls) { integrand ->
val absoluteAccuracy = integrand.getFeature<IntegrandAbsoluteAccuracy>()?.accuracy
val absoluteAccuracy = integrand[IntegrandAbsoluteAccuracy]
?: IterativeLegendreGaussIntegrator.DEFAULT_ABSOLUTE_ACCURACY
val relativeAccuracy = integrand.getFeature<IntegrandRelativeAccuracy>()?.accuracy
val relativeAccuracy = integrand[IntegrandRelativeAccuracy]
?: IterativeLegendreGaussIntegrator.DEFAULT_ABSOLUTE_ACCURACY
val iterations = integrand.getFeature<IntegrandIterationsRange>()?.range
val iterations = integrand[IntegrandIterationsRange]
?: IterativeLegendreGaussIntegrator.DEFAULT_MIN_ITERATIONS_COUNT..IterativeLegendreGaussIntegrator.DEFAULT_MAX_ITERATIONS_COUNT
IterativeLegendreGaussIntegrator(

2
kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt
View File

@ -135,7 +135,7 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {
override val r: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(qr.r).withAttribute(UpperTriangular) }
}
SingularValueDecompositionAttribute::class -> object : SingularValueDecompositionAttribute<Double> {
SVDAttribute::class -> object : SVDAttribute<Double> {
private val sv by lazy { SingularValueDecomposition(origin) }
override val u: Matrix<Double> by lazy { CMMatrix(sv.u) }
override val s: Matrix<Double> by lazy { CMMatrix(sv.s) }

9
kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt
View File

@ -13,10 +13,6 @@ import space.kscience.attributes.safeTypeOf
import space.kscience.kmath.UnstableKMathAPI
import space.kscience.kmath.operations.*
import space.kscience.kmath.structures.*
import space.kscience.kmath.structures.BufferAccessor2D
import space.kscience.kmath.structures.Float64Buffer
import space.kscience.kmath.structures.MutableBuffer
import space.kscience.kmath.structures.MutableBufferFactory
/**
* Matrices with this feature support LU factorization with partial pivoting: *[p] &middot; a = [l] &middot; [u]* where
@ -197,8 +193,7 @@ public fun <T : Comparable<T>> LinearSpace<T, Field<T>>.lup(
public fun LinearSpace<Double, Float64Field>.lup(
matrix: Matrix<Double>,
singularityThreshold: Double = 1e-11,
): LupDecomposition<Double> =
lup(::Float64Buffer, matrix) { it < singularityThreshold }
): LupDecomposition<Double> = lup(matrix) { it < singularityThreshold }
internal fun <T : Any, A : Field<T>> LinearSpace<T, A>.solve(
lup: LupDecomposition<T>,
@ -265,4 +260,4 @@ public fun <T : Comparable<T>, F : Field<T>> LinearSpace<T, F>.lupSolver(
}
public fun LinearSpace<Double, Float64Field>.lupSolver(singularityThreshold: Double = 1e-11): LinearSolver<Double> =
lupSolver(::Float64Buffer) { it < singularityThreshold }
lupSolver { it < singularityThreshold }

6
kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt → kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/matrixAttributes.kt
View File

@ -163,12 +163,12 @@ public interface SingularValueDecomposition<T> {
*
* @param T the type of matrices' items.
*/
public class SingularValueDecompositionAttribute<T>(type: SafeType<SingularValueDecomposition<T>>) :
public class SVDAttribute<T>(type: SafeType<SingularValueDecomposition<T>>) :
PolymorphicAttribute<SingularValueDecomposition<T>>(type),
MatrixAttribute<SingularValueDecomposition<T>>
public val <T> MatrixOperations<T>.SVD: SingularValueDecompositionAttribute<T>
get() = SingularValueDecompositionAttribute(safeTypeOf())
public val <T> MatrixOperations<T>.SVD: SVDAttribute<T>
get() = SVDAttribute(safeTypeOf())
//TODO add sparse matrix feature

1
kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt
View File

@ -240,7 +240,6 @@ public interface MutableStructureND<T> : StructureND<T> {
* Set value at specified indices
*/
@PerformancePitfall
@Deprecated("")
public operator fun <T> MutableStructureND<T>.set(vararg index: Int, value: T) {
set(index, value)
}

7
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/polynomialConstructors.kt
View File

@ -20,9 +20,4 @@ public fun <C> Polynomial(coefficients: List<C>, reverse: Boolean = false): Poly
*/
@Suppress("FunctionName")
public fun <C> Polynomial(vararg coefficients: C, reverse: Boolean = false): Polynomial<C> =
Polynomial(with(coefficients) { if (reverse) reversed() else toList() })
/**
* Represents [this] constant as a [Polynomial].
*/
public fun <C> C.asPolynomial() : Polynomial<C> = Polynomial(listOf(this))
Polynomial(with(coefficients) { if (reverse) reversed() else toList() })

34
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt
View File

@ -4,6 +4,7 @@
*/
package space.kscience.kmath.integration
import space.kscience.attributes.AttributesBuilder
import space.kscience.kmath.UnstableKMathAPI
import space.kscience.kmath.operations.Field
import space.kscience.kmath.structures.Buffer
@ -11,14 +12,14 @@ import space.kscience.kmath.structures.asBuffer
import space.kscience.kmath.structures.indices
/**
* A simple one-pass integrator based on Gauss rule
* Following integrand features are accepted:
* A simple one-pass integrator based on Gauss rule.
* The following integrand features are accepted:
*
* * [GaussIntegratorRuleFactory]&mdash;a factory for computing the Gauss integration rule. By default, uses
* [GaussLegendreRuleFactory].
* * [IntegrationRange]&mdash;the univariate range of integration. By default, uses `0..1` interval.
* * [IntegrandMaxCalls]&mdash;the maximum number of function calls during integration. For non-iterative rules, always
* uses the maximum number of points. By default, uses 10 points.
* use the maximum number of points. By default, uses 10 points.
* * [UnivariateIntegrandRanges]&mdash;set of ranges and number of points per range. Defaults to given
* [IntegrationRange] and [IntegrandMaxCalls].
*/
@ -27,11 +28,11 @@ public class GaussIntegrator<T : Any>(
) : UnivariateIntegrator<T> {
private fun buildRule(integrand: UnivariateIntegrand<T>): Pair<Buffer<Double>, Buffer<Double>> {
val factory = integrand.getFeature<GaussIntegratorRuleFactory>() ?: GaussLegendreRuleFactory
val predefinedRanges = integrand.getFeature<UnivariateIntegrandRanges>()
val factory = integrand[GaussIntegratorRuleFactory] ?: GaussLegendreRuleFactory
val predefinedRanges = integrand[UnivariateIntegrandRanges]
if (predefinedRanges == null || predefinedRanges.ranges.isEmpty()) {
val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
val numPoints = integrand[IntegrandMaxCalls] ?: 100
val range = integrand[IntegrationRange] ?: 0.0..1.0
return factory.build(numPoints, range)
} else {
val ranges = predefinedRanges.ranges
@ -66,7 +67,10 @@ public class GaussIntegrator<T : Any>(
c = t - res - y
res = t
}
return integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + points.size)
return integrand.modify {
value(res)
IntegrandCallsPerformed(integrand.calls + points.size)
}
}
public companion object
@ -88,7 +92,7 @@ public fun <T : Any> GaussIntegrator<T>.integrate(
range: ClosedRange<Double>,
order: Int = 10,
intervals: Int = 10,
vararg features: IntegrandFeature,
attributesBuilder: AttributesBuilder.() -> Unit,
function: (Double) -> T,
): UnivariateIntegrand<T> {
require(range.endInclusive > range.start) { "The range upper bound should be higher than lower bound" }
@ -100,11 +104,13 @@ public fun <T : Any> GaussIntegrator<T>.integrate(
)
return process(
UnivariateIntegrand(
function,
IntegrationRange(range),
GaussLegendreRuleFactory,
ranges,
*features
attributeBuilder = {
IntegrationRange(range)
GaussIntegratorRuleFactory(GaussLegendreRuleFactory)
UnivariateIntegrandRanges(ranges)
attributesBuilder()
},
function = function,
)
)
}

4
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegratorRuleFactory.kt
View File

@ -12,10 +12,10 @@ import space.kscience.kmath.structures.asBuffer
import kotlin.math.ulp
import kotlin.native.concurrent.ThreadLocal
public interface GaussIntegratorRuleFactory : IntegrandFeature {
public interface GaussIntegratorRuleFactory {
public fun build(numPoints: Int): Pair<Buffer<Double>, Buffer<Double>>
public companion object {
public companion object: IntegrandAttribute<GaussIntegratorRuleFactory>{
public fun double(numPoints: Int, range: ClosedRange<Double>): Pair<Buffer<Double>, Buffer<Double>> =
GaussLegendreRuleFactory.build(numPoints, range)
}

63
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrand.kt
View File

@ -5,44 +5,49 @@
package space.kscience.kmath.integration
import space.kscience.kmath.misc.Feature
import space.kscience.kmath.misc.FeatureSet
import space.kscience.kmath.misc.Featured
import kotlin.reflect.KClass
import space.kscience.attributes.*
import kotlin.reflect.typeOf
public interface IntegrandFeature : Feature<IntegrandFeature> {
override fun toString(): String
public interface IntegrandAttribute<T> : Attribute<T>
public interface Integrand<T> : AttributeContainer {
public fun modify(block: AttributesBuilder.() -> Unit): Integrand<T>
public fun <A : Any> withAttribute(attribute: Attribute<A>, value: A): Integrand<T>
public companion object
}
public interface Integrand : Featured<IntegrandFeature> {
public val features: FeatureSet<IntegrandFeature>
override fun <T : IntegrandFeature> getFeature(type: KClass<out T>): T? = features.getFeature(type)
public operator fun <T> Integrand<*>.get(attribute: Attribute<T>): T? = attributes[attribute]
public class IntegrandValue<T>(type: SafeType<T>) : PolymorphicAttribute<T>(type), IntegrandAttribute<T>
public inline val <reified T : Any> Integrand<T>.Value: IntegrandValue<T> get() = IntegrandValue(safeTypeOf())
public fun <T> AttributesBuilder.value(value: T){
val type: SafeType<T> = typeOf<T>()
IntegrandValue(type).invoke(value)
}
public inline fun <reified T : IntegrandFeature> Integrand.getFeature(): T? = getFeature(T::class)
/**
* Value of the integrand if it is present or null
*/
public inline val <reified T : Any> Integrand<T>.valueOrNull: T? get() = attributes[Value]
public class IntegrandValue<out T : Any>(public val value: T) : IntegrandFeature {
override fun toString(): String = "Value($value)"
}
/**
* Value of the integrand or error
*/
public inline val <reified T : Any> Integrand<T>.value: T get() = valueOrNull ?: error("No value in the integrand")
public class IntegrandRelativeAccuracy(public val accuracy: Double) : IntegrandFeature {
override fun toString(): String = "TargetRelativeAccuracy($accuracy)"
}
public object IntegrandRelativeAccuracy : IntegrandAttribute<Double>
public class IntegrandAbsoluteAccuracy(public val accuracy: Double) : IntegrandFeature {
override fun toString(): String = "TargetAbsoluteAccuracy($accuracy)"
}
public object IntegrandAbsoluteAccuracy : IntegrandAttribute<Double>
public class IntegrandCallsPerformed(public val calls: Int) : IntegrandFeature {
override fun toString(): String = "Calls($calls)"
}
public object IntegrandCallsPerformed : IntegrandAttribute<Int>
public val Integrand.calls: Int get() = getFeature<IntegrandCallsPerformed>()?.calls ?: 0
public val Integrand<*>.calls: Int get() = attributes[IntegrandCallsPerformed] ?: 0
public class IntegrandMaxCalls(public val maxCalls: Int) : IntegrandFeature {
override fun toString(): String = "MaxCalls($maxCalls)"
}
public object IntegrandMaxCalls : IntegrandAttribute<Int>
public class IntegrandIterationsRange(public val range: IntRange) : IntegrandFeature {
override fun toString(): String = "Iterations(${range.first}..${range.last})"
}
public object IntegrandIterationsRange : IntegrandAttribute<IntRange>

22
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/MultivariateIntegrand.kt
View File

@ -5,22 +5,22 @@
package space.kscience.kmath.integration
import space.kscience.attributes.*
import space.kscience.kmath.linear.Point
import space.kscience.kmath.misc.FeatureSet
public class MultivariateIntegrand<T : Any> internal constructor(
override val features: FeatureSet<IntegrandFeature>,
public class MultivariateIntegrand<T> internal constructor(
override val attributes: Attributes,
public val function: (Point<T>) -> T,
) : Integrand {
) : Integrand<T> {
public operator fun <F : IntegrandFeature> plus(feature: F): MultivariateIntegrand<T> =
MultivariateIntegrand(features.with(feature), function)
override fun modify(block: AttributesBuilder.() -> Unit): MultivariateIntegrand<T> =
MultivariateIntegrand(attributes.modify(block), function)
override fun <A : Any> withAttribute(attribute: Attribute<A>, value: A): MultivariateIntegrand<T> =
MultivariateIntegrand(attributes.withAttribute(attribute, value), function)
}
@Suppress("FunctionName")
public fun <T : Any> MultivariateIntegrand(
vararg features: IntegrandFeature,
attributeBuilder: AttributesBuilder.() -> Unit,
function: (Point<T>) -> T,
): MultivariateIntegrand<T> = MultivariateIntegrand(FeatureSet.of(*features), function)
public val <T : Any> MultivariateIntegrand<T>.value: T? get() = getFeature<IntegrandValue<T>>()?.value
): MultivariateIntegrand<T> = MultivariateIntegrand(Attributes(attributeBuilder), function)

32
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SimpsonIntegrator.kt
View File

@ -45,16 +45,22 @@ public class SimpsonIntegrator<T : Any>(
}
override fun process(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> {
val ranges = integrand.getFeature<UnivariateIntegrandRanges>()
val ranges = integrand[UnivariateIntegrandRanges]
return if (ranges != null) {
val res = algebra.sum(ranges.ranges.map { integrateRange(integrand, it.first, it.second) })
integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second })
integrand.modify {
value(res)
IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second })
}
} else {
val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
val numPoints = integrand[IntegrandMaxCalls] ?: 100
require(numPoints >= 4) { "Simpson integrator requires at least 4 nodes" }
val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
val range = integrand[IntegrationRange] ?: 0.0..1.0
val res = integrateRange(integrand, range, numPoints)
integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + numPoints)
integrand.modify {
value(res)
IntegrandCallsPerformed(integrand.calls + numPoints)
}
}
}
}
@ -91,16 +97,22 @@ public object DoubleSimpsonIntegrator : UnivariateIntegrator<Double> {
}
override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
val ranges = integrand.getFeature<UnivariateIntegrandRanges>()
val ranges = integrand[UnivariateIntegrandRanges]
return if (ranges != null) {
val res = ranges.ranges.sumOf { integrateRange(integrand, it.first, it.second) }
integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second })
integrand.modify {
value(res)
IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second })
}
} else {
val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
val numPoints = integrand[IntegrandMaxCalls] ?: 100
require(numPoints >= 4) { "Simpson integrator requires at least 4 nodes" }
val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
val range = integrand[IntegrationRange] ?: 0.0..1.0
val res = integrateRange(integrand, range, numPoints)
integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + numPoints)
integrand.modify {
value(res)
IntegrandCallsPerformed(integrand.calls + numPoints)
}
}
}
}

22
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt
View File

@ -55,12 +55,12 @@ public class SplineIntegrator<T : Comparable<T>>(
public val bufferFactory: MutableBufferFactory<T>,
) : UnivariateIntegrator<T> {
override fun process(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = algebra {
val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
val range = integrand[IntegrationRange] ?: 0.0..1.0
val interpolator: PolynomialInterpolator<T> = SplineInterpolator(algebra, bufferFactory)
val nodes: Buffer<Double> = integrand.getFeature<UnivariateIntegrationNodes>()?.nodes ?: run {
val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
val nodes: Buffer<Double> = integrand[UnivariateIntegrationNodes] ?: run {
val numPoints = integrand[IntegrandMaxCalls] ?: 100
val step = (range.endInclusive - range.start) / (numPoints - 1)
Float64Buffer(numPoints) { i -> range.start + i * step }
}
@ -71,7 +71,10 @@ public class SplineIntegrator<T : Comparable<T>>(
values
)
val res = polynomials.integrate(algebra, number(range.start)..number(range.endInclusive))
integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + nodes.size)
integrand.modify {
value(res)
IntegrandCallsPerformed(integrand.calls + nodes.size)
}
}
}
@ -84,11 +87,11 @@ public class SplineIntegrator<T : Comparable<T>>(
@UnstableKMathAPI
public object DoubleSplineIntegrator : UnivariateIntegrator<Double> {
override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
val range = integrand[IntegrationRange] ?: 0.0..1.0
val interpolator: PolynomialInterpolator<Double> = SplineInterpolator(Float64Field, ::Float64Buffer)
val nodes: Buffer<Double> = integrand.getFeature<UnivariateIntegrationNodes>()?.nodes ?: run {
val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
val nodes: Buffer<Double> = integrand[UnivariateIntegrationNodes] ?: run {
val numPoints = integrand[IntegrandMaxCalls] ?: 100
val step = (range.endInclusive - range.start) / (numPoints - 1)
Float64Buffer(numPoints) { i -> range.start + i * step }
}
@ -96,7 +99,10 @@ public object DoubleSplineIntegrator : UnivariateIntegrator<Double> {
val values = nodes.mapToBuffer(::Float64Buffer) { integrand.function(it) }
val polynomials = interpolator.interpolatePolynomials(nodes, values)
val res = polynomials.integrate(Float64Field, range)
return integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + nodes.size)
return integrand.modify {
value(res)
IntegrandCallsPerformed(integrand.calls + nodes.size)
}
}
}

73
kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt
View File

@ -5,36 +5,38 @@
package space.kscience.kmath.integration
import space.kscience.attributes.*
import space.kscience.kmath.UnstableKMathAPI
import space.kscience.kmath.misc.FeatureSet
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.Float64Buffer
public class UnivariateIntegrand<T> internal constructor(
override val features: FeatureSet<IntegrandFeature>,
override val attributes: Attributes,
public val function: (Double) -> T,
) : Integrand {
public operator fun <F : IntegrandFeature> plus(feature: F): UnivariateIntegrand<T> =
UnivariateIntegrand(features.with(feature), function)
) : Integrand<T> {
override fun <A : Any> withAttribute(attribute: Attribute<A>, value: A): UnivariateIntegrand<T> =
UnivariateIntegrand(attributes.withAttribute(attribute, value), function)
override fun modify(block: AttributesBuilder.() -> Unit): UnivariateIntegrand<T> =
UnivariateIntegrand(attributes.modify(block), function)
}
@Suppress("FunctionName")
public fun <T : Any> UnivariateIntegrand(
attributeBuilder: AttributesBuilder.() -> Unit,
function: (Double) -> T,
vararg features: IntegrandFeature,
): UnivariateIntegrand<T> = UnivariateIntegrand(FeatureSet.of(*features), function)
): UnivariateIntegrand<T> = UnivariateIntegrand(Attributes(attributeBuilder), function)
public typealias UnivariateIntegrator<T> = Integrator<UnivariateIntegrand<T>>
public class IntegrationRange(public val range: ClosedRange<Double>) : IntegrandFeature {
override fun toString(): String = "Range(${range.start}..${range.endInclusive})"
}
public object IntegrationRange : IntegrandAttribute<ClosedRange<Double>>
/**
* Set of univariate integration ranges. First components correspond to the ranges themselves, second components to
* number of integration nodes per range.
* the number of integration nodes per range.
*/
public class UnivariateIntegrandRanges(public val ranges: List<Pair<ClosedRange<Double>, Int>>) : IntegrandFeature {
public class UnivariateIntegrandRanges(public val ranges: List<Pair<ClosedRange<Double>, Int>>) {
public constructor(vararg pairs: Pair<ClosedRange<Double>, Int>) : this(pairs.toList())
override fun toString(): String {
@ -43,45 +45,25 @@ public class UnivariateIntegrandRanges(public val ranges: List<Pair<ClosedRange<
}
return "UnivariateRanges($rangesString)"
}
public companion object : IntegrandAttribute<UnivariateIntegrandRanges>
}
public class UnivariateIntegrationNodes(public val nodes: Buffer<Double>) : IntegrandFeature {
public constructor(vararg nodes: Double) : this(Float64Buffer(nodes))
public object UnivariateIntegrationNodes : IntegrandAttribute<Buffer<Double>>
override fun toString(): String = "UnivariateNodes($nodes)"
public fun AttributesBuilder.integrationNodes(vararg nodes: Double) {
UnivariateIntegrationNodes(Float64Buffer(nodes))
}
/**
* Value of the integrand if it is present or null
*/
public val <T : Any> UnivariateIntegrand<T>.valueOrNull: T? get() = getFeature<IntegrandValue<T>>()?.value
/**
* Value of the integrand or error
*/
public val <T : Any> UnivariateIntegrand<T>.value: T get() = valueOrNull ?: error("No value in the integrand")
/**
* A shortcut method to integrate a [function] with additional [features]. Range must be provided in features.
* The [function] is placed in the end position to allow passing a lambda.
*/
@UnstableKMathAPI
public fun <T : Any> UnivariateIntegrator<T>.integrate(
vararg features: IntegrandFeature,
attributesBuilder: AttributesBuilder.() -> Unit,
function: (Double) -> T,
): UnivariateIntegrand<T> = process(UnivariateIntegrand(function, *features))
/**
* 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 <T : Any> UnivariateIntegrator<T>.integrate(
range: ClosedRange<Double>,
vararg features: IntegrandFeature,
function: (Double) -> T,
): UnivariateIntegrand<T> = process(UnivariateIntegrand(function, IntegrationRange(range), *features))
): UnivariateIntegrand<T> = process(UnivariateIntegrand(attributesBuilder, function))
/**
* A shortcut method to integrate a [function] in [range] with additional features.
@ -90,13 +72,12 @@ public fun <T : Any> UnivariateIntegrator<T>.integrate(
@UnstableKMathAPI
public fun <T : Any> UnivariateIntegrator<T>.integrate(
range: ClosedRange<Double>,
featureBuilder: MutableList<IntegrandFeature>.() -> Unit = {},
attributeBuilder: AttributesBuilder.() -> Unit = {},
function: (Double) -> T,
): UnivariateIntegrand<T> {
//TODO use dedicated feature builder class instead or add extensions to MutableList<IntegrandFeature>
val features = buildList {
featureBuilder()
add(IntegrationRange(range))
val attributes = Attributes {
IntegrationRange(range)
attributeBuilder()
}
return process(UnivariateIntegrand(function, *features.toTypedArray()))
return process(UnivariateIntegrand(attributes, function))
}

10
kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SimpsonIntegralTest.kt
View File

@ -16,7 +16,10 @@ import kotlin.test.assertEquals
class SimpsonIntegralTest {
@Test
fun gaussSin() {
val res = Float64Field.simpsonIntegrator.integrate(0.0..2 * PI, IntegrandMaxCalls(5)) { x ->
val res = Float64Field.simpsonIntegrator.integrate(
0.0..2 * PI,
{ IntegrandMaxCalls(5) }
) { x ->
sin(x)
}
assertEquals(0.0, res.value, 1e-2)
@ -24,7 +27,10 @@ class SimpsonIntegralTest {
@Test
fun gaussUniform() {
val res = Float64Field.simpsonIntegrator.integrate(35.0..100.0, IntegrandMaxCalls(20)) { x ->
val res = Float64Field.simpsonIntegrator.integrate(
35.0..100.0,
{ IntegrandMaxCalls(20) }
) { x ->
if (x in 30.0..50.0) {
1.0
} else {

2
kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/testUtils.kt
View File

@ -27,7 +27,7 @@ fun grid(
return xs.flatMap { x -> ys.map { y -> x to y } }
}
fun assertVectorEquals(expected: DoubleVector2D, actual: DoubleVector2D, absoluteTolerance: Double = 1e-6) {
fun assertVectorEquals(expected: DoubleVector2D, actual: DoubleVector2D, absoluteTolerance: Double = 1e-3) {
assertEquals(expected.x, actual.x, absoluteTolerance)
assertEquals(expected.y, actual.y, absoluteTolerance)
}

2
kmath-optimization/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt
View File

@ -53,7 +53,7 @@ public class OptimizationPrior<T>(type: SafeType<T>):
PolymorphicAttribute<DifferentiableExpression<T>>(safeTypeOf()),
Attribute<DifferentiableExpression<T>>
public val <T> FunctionOptimization.Companion.Optimization get() =
//public val <T> FunctionOptimization.Companion.Optimization get() =
public fun <T> FunctionOptimization<T>.withFeatures(

10
kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/LevenbergMarquardtAlgorithm.kt
View File

@ -6,7 +6,7 @@
package space.kscience.kmath.tensors.core
import space.kscience.kmath.PerformancePitfall
import space.kscience.kmath.linear.transpose
import space.kscience.kmath.linear.transposed
import space.kscience.kmath.nd.*
import kotlin.math.abs
import kotlin.math.max
@ -139,7 +139,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
if (inputData.nargin < 5) {
weight = fromArray(
ShapeND(intArrayOf(1, 1)),
doubleArrayOf((inputData.realValues.transpose().dot(inputData.realValues)).as1D()[0])
doubleArrayOf((inputData.realValues.transposed dot inputData.realValues).as1D()[0])
).as2D()
}
@ -266,12 +266,12 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
settings.funcCalls += 1
// val tmp = deltaY.times(weight)
var X2Try = deltaY.as2D().transpose().dot(deltaY.times(weight)) // Chi-squared error criteria
var X2Try = deltaY.as2D().transposed dot deltaY.times(weight) // Chi-squared error criteria
val alpha = 1.0
if (updateType == 2) { // Quadratic
// One step of quadratic line update in the h direction for minimum X2
val alphaTensor = (jtWdy.transpose() dot h) / ((X2Try - x2) / 2.0 + 2 * (jtWdy.transpose() dot h))
val alphaTensor = (jtWdy.transposed dot h) / ((X2Try - x2) / 2.0 + 2 * (jtWdy.transposed dot h))
h = h dot alphaTensor
pTry = (p + h).as2D() // update only [idx] elements
pTry = smallestElementComparison(
@ -289,7 +289,7 @@ public fun DoubleTensorAlgebra.levenbergMarquardt(inputData: LMInput): LMResultI
) // residual error using p_try
settings.funcCalls += 1
X2Try = deltaY.as2D().transpose() dot deltaY * weight // Chi-squared error criteria
X2Try = deltaY.as2D().transposed dot deltaY * weight // Chi-squared error criteria
}
val rho = when (updateType) { // Nielsen