This commit is contained in:
Alexander Nozik 2023-11-01 08:55:47 +03:00
parent 544b8610e1
commit ea887b8c72
8 changed files with 45 additions and 66 deletions

View File

@ -19,13 +19,19 @@ import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC
import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC
import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC
import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC
import space.kscience.kmath.UnstableKMathAPI
import space.kscience.kmath.linear.* import space.kscience.kmath.linear.*
import space.kscience.kmath.linear.Matrix import space.kscience.kmath.linear.Matrix
import space.kscience.kmath.UnstableKMathAPI
import space.kscience.kmath.nd.StructureFeature import space.kscience.kmath.nd.StructureFeature
import space.kscience.kmath.operations.Float32Field import space.kscience.kmath.structures.Float64
import space.kscience.kmath.structures.Float32
import space.kscience.kmath.operations.Float64Field import space.kscience.kmath.operations.Float64Field
import space.kscience.kmath.operations.Float32Field
import space.kscience.kmath.operations.DoubleField
import space.kscience.kmath.operations.FloatField
import space.kscience.kmath.operations.invoke import space.kscience.kmath.operations.invoke
import space.kscience.kmath.structures.Float64Buffer
import space.kscience.kmath.structures.Float32Buffer
import space.kscience.kmath.structures.DoubleBuffer import space.kscience.kmath.structures.DoubleBuffer
import space.kscience.kmath.structures.FloatBuffer import space.kscience.kmath.structures.FloatBuffer
import kotlin.reflect.KClass import kotlin.reflect.KClass

View File

@ -54,7 +54,7 @@ public class GaussIntegrator<T : Any>(
} }
} }
override fun process(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = with(algebra) { override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = with(algebra) {
val f = integrand.function val f = integrand.function
val (points, weights) = buildRule(integrand) val (points, weights) = buildRule(integrand)
var res: T = zero var res: T = zero
@ -68,7 +68,7 @@ public class GaussIntegrator<T : Any>(
res = t res = t
} }
return integrand.withAttributes { return integrand.withAttributes {
value(res) IntegrandValue(res)
IntegrandCallsPerformed(integrand.calls + points.size) IntegrandCallsPerformed(integrand.calls + points.size)
} }
} }
@ -102,7 +102,7 @@ public inline fun <reified T : Any> GaussIntegrator<T>.integrate(
val ranges = UnivariateIntegrandRanges( val ranges = UnivariateIntegrandRanges(
(0 until intervals).map { i -> (range.start + rangeSize * i)..(range.start + rangeSize * (i + 1)) to order } (0 until intervals).map { i -> (range.start + rangeSize * i)..(range.start + rangeSize * (i + 1)) to order }
) )
return process( return integrate(
UnivariateIntegrand<T>( UnivariateIntegrand<T>(
attributeBuilder = { attributeBuilder = {
IntegrationRange(range) IntegrationRange(range)

View File

@ -8,9 +8,9 @@ package space.kscience.kmath.integration
/** /**
* A general interface for all integrators. * A general interface for all integrators.
*/ */
public interface Integrator<I : Integrand> { public interface Integrator<T, I : Integrand<T>> {
/** /**
* Runs one integration pass and return a new [Integrand] with a new set of features. * Runs one integration pass and return a new [Integrand] with a new set of features.
*/ */
public fun process(integrand: I): I public fun integrate(integrand: I): I
} }

View File

@ -44,12 +44,12 @@ public class SimpsonIntegrator<T : Any>(
return res return res
} }
override fun process(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> { override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> {
val ranges = integrand[UnivariateIntegrandRanges] val ranges = integrand[UnivariateIntegrandRanges]
return if (ranges != null) { return if (ranges != null) {
val res = algebra.sum(ranges.ranges.map { integrateRange(integrand, it.first, it.second) }) val res = algebra.sum(ranges.ranges.map { integrateRange(integrand, it.first, it.second) })
integrand.withAttributes { integrand.withAttributes {
value(res) IntegrandValue(res)
IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second }) IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second })
} }
} else { } else {
@ -58,7 +58,7 @@ public class SimpsonIntegrator<T : Any>(
val range = integrand[IntegrationRange] ?: 0.0..1.0 val range = integrand[IntegrationRange] ?: 0.0..1.0
val res = integrateRange(integrand, range, numPoints) val res = integrateRange(integrand, range, numPoints)
integrand.withAttributes { integrand.withAttributes {
value(res) IntegrandValue(res)
IntegrandCallsPerformed(integrand.calls + numPoints) IntegrandCallsPerformed(integrand.calls + numPoints)
} }
} }
@ -96,12 +96,12 @@ public object DoubleSimpsonIntegrator : UnivariateIntegrator<Double> {
return res return res
} }
override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> { override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
val ranges = integrand[UnivariateIntegrandRanges] val ranges = integrand[UnivariateIntegrandRanges]
return if (ranges != null) { return if (ranges != null) {
val res = ranges.ranges.sumOf { integrateRange(integrand, it.first, it.second) } val res = ranges.ranges.sumOf { integrateRange(integrand, it.first, it.second) }
integrand.withAttributes { integrand.withAttributes {
value(res) IntegrandValue(res)
IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second }) IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second })
} }
} else { } else {
@ -110,7 +110,7 @@ public object DoubleSimpsonIntegrator : UnivariateIntegrator<Double> {
val range = integrand[IntegrationRange] ?: 0.0..1.0 val range = integrand[IntegrationRange] ?: 0.0..1.0
val res = integrateRange(integrand, range, numPoints) val res = integrateRange(integrand, range, numPoints)
integrand.withAttributes { integrand.withAttributes {
value(res) IntegrandValue(res)
IntegrandCallsPerformed(integrand.calls + numPoints) IntegrandCallsPerformed(integrand.calls + numPoints)
} }
} }

View File

@ -54,7 +54,7 @@ public class SplineIntegrator<T : Comparable<T>>(
public val algebra: Field<T>, public val algebra: Field<T>,
public val bufferFactory: MutableBufferFactory<T>, public val bufferFactory: MutableBufferFactory<T>,
) : UnivariateIntegrator<T> { ) : UnivariateIntegrator<T> {
override fun process(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = algebra { override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = algebra {
val range = integrand[IntegrationRange] ?: 0.0..1.0 val range = integrand[IntegrationRange] ?: 0.0..1.0
val interpolator: PolynomialInterpolator<T> = SplineInterpolator(algebra, bufferFactory) val interpolator: PolynomialInterpolator<T> = SplineInterpolator(algebra, bufferFactory)
@ -72,7 +72,7 @@ public class SplineIntegrator<T : Comparable<T>>(
) )
val res = polynomials.integrate(algebra, number(range.start)..number(range.endInclusive)) val res = polynomials.integrate(algebra, number(range.start)..number(range.endInclusive))
integrand.withAttributes { integrand.withAttributes {
value(res) IntegrandValue(res)
IntegrandCallsPerformed(integrand.calls + nodes.size) IntegrandCallsPerformed(integrand.calls + nodes.size)
} }
} }
@ -86,7 +86,7 @@ public class SplineIntegrator<T : Comparable<T>>(
*/ */
@UnstableKMathAPI @UnstableKMathAPI
public object DoubleSplineIntegrator : UnivariateIntegrator<Double> { public object DoubleSplineIntegrator : UnivariateIntegrator<Double> {
override fun process(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> { override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
val range = integrand[IntegrationRange] ?: 0.0..1.0 val range = integrand[IntegrationRange] ?: 0.0..1.0
val interpolator: PolynomialInterpolator<Double> = SplineInterpolator(Float64Field, ::Float64Buffer) val interpolator: PolynomialInterpolator<Double> = SplineInterpolator(Float64Field, ::Float64Buffer)
@ -100,7 +100,7 @@ public object DoubleSplineIntegrator : UnivariateIntegrator<Double> {
val polynomials = interpolator.interpolatePolynomials(nodes, values) val polynomials = interpolator.interpolatePolynomials(nodes, values)
val res = polynomials.integrate(Float64Field, range) val res = polynomials.integrate(Float64Field, range)
return integrand.withAttributes { return integrand.withAttributes {
value(res) IntegrandValue(res)
IntegrandCallsPerformed(integrand.calls + nodes.size) IntegrandCallsPerformed(integrand.calls + nodes.size)
} }
} }

View File

@ -34,7 +34,7 @@ public inline fun <reified T : Any> UnivariateIntegrand(
noinline function: (Double) -> T, noinline function: (Double) -> T,
): UnivariateIntegrand<T> = UnivariateIntegrand(safeTypeOf(), Attributes(attributeBuilder), function) ): UnivariateIntegrand<T> = UnivariateIntegrand(safeTypeOf(), Attributes(attributeBuilder), function)
public typealias UnivariateIntegrator<T> = Integrator<UnivariateIntegrand<T>> public typealias UnivariateIntegrator<T> = Integrator<T, UnivariateIntegrand<T>>
public object IntegrationRange : IntegrandAttribute<ClosedRange<Double>> public object IntegrationRange : IntegrandAttribute<ClosedRange<Double>>
@ -70,7 +70,7 @@ public fun TypedAttributesBuilder<UnivariateIntegrand<*>>.integrationNodes(varar
public inline fun <reified T : Any> UnivariateIntegrator<T>.integrate( public inline fun <reified T : Any> UnivariateIntegrator<T>.integrate(
attributesBuilder: TypedAttributesBuilder<UnivariateIntegrand<T>>.() -> Unit, attributesBuilder: TypedAttributesBuilder<UnivariateIntegrand<T>>.() -> Unit,
noinline function: (Double) -> T, noinline function: (Double) -> T,
): UnivariateIntegrand<T> = process(UnivariateIntegrand(attributesBuilder, function)) ): UnivariateIntegrand<T> = integrate(UnivariateIntegrand(attributesBuilder, function))
/** /**
* A shortcut method to integrate a [function] in [range] with additional features. * A shortcut method to integrate a [function] in [range] with additional features.
@ -83,7 +83,7 @@ public inline fun <reified T : Any> UnivariateIntegrator<T>.integrate(
noinline function: (Double) -> T, noinline function: (Double) -> T,
): UnivariateIntegrand<T> { ): UnivariateIntegrand<T> {
return process( return integrate(
UnivariateIntegrand( UnivariateIntegrand(
attributeBuilder = { attributeBuilder = {
IntegrationRange(range) IntegrationRange(range)

View File

@ -31,24 +31,3 @@ public fun QuaternionAlgebra.angleBetween(q1: Quaternion, q2: Quaternion): Angle
* Euclidean product of two quaternions * Euclidean product of two quaternions
*/ */
public infix fun Quaternion.dot(other: Quaternion): Double = w * other.w + x * other.x + y * other.y + z * other.z public infix fun Quaternion.dot(other: Quaternion): Double = w * other.w + x * other.x + y * other.y + z * other.z
//
///**
// * Convert a quaternion to XYZ Cardan angles assuming it is normalized.
// */
//private fun Quaternion.normalizedToEuler(): Float32Vector3D {
// val roll = atan2(2 * y * w + 2 * x * z, 1 - 2 * y * y - 2 * z * z)
// val pitch = atan2(2 * x * w - 2 * y * z, 1 - 2 * x * x - 2 * z * z)
// val yaw = asin(2 * x * y + 2 * z * w)
//
// return Float32Vector3D(roll, pitch, yaw)
//}
///**
// * Quaternion to XYZ Cardan angles
// */
//public fun Quaternion.toEuler(): Float32Vector3D = if (QuaternionAlgebra.norm(this) == 0.0) {
// Float32Space3D.zero
//} else {
// normalized().normalizedToEuler()
//}

View File

@ -310,7 +310,11 @@ internal fun DoubleTensorAlgebra.svdHelper(
} }
} }
private fun pythag(a: Double, b: Double): Double { internal fun MutableStructure2D<Double>.svdGolubKahanHelper(
u: MutableStructure2D<Double>, w: BufferedTensor<Double>,
v: MutableStructure2D<Double>, iterations: Int, epsilon: Double,
) {
fun pythag(a: Double, b: Double): Double {
val at: Double = abs(a) val at: Double = abs(a)
val bt: Double = abs(b) val bt: Double = abs(b)
val ct: Double val ct: Double
@ -323,19 +327,9 @@ private fun pythag(a: Double, b: Double): Double {
result = bt * sqrt(1.0 + ct * ct) result = bt * sqrt(1.0 + ct * ct)
} else result = 0.0 } else result = 0.0
return result return result
} }
private fun SIGN(a: Double, b: Double): Double {
if (b >= 0.0)
return abs(a)
else
return -abs(a)
}
internal fun MutableStructure2D<Double>.svdGolubKahanHelper(
u: MutableStructure2D<Double>, w: BufferedTensor<Double>,
v: MutableStructure2D<Double>, iterations: Int, epsilon: Double,
) {
val shape = this.shape val shape = this.shape
val m = shape.component1() val m = shape.component1()
val n = shape.component2() val n = shape.component2()
@ -553,7 +547,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(
h = rv1[k] h = rv1[k]
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y) f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y)
g = pythag(f, 1.0) g = pythag(f, 1.0)
f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x f = ((x - z) * (x + z) + h * ((y / (f + if (f >= 0.0) abs(g) else -abs(g) )) - h)) / x
c = 1.0 c = 1.0
s = 1.0 s = 1.0
@ -563,7 +557,7 @@ internal fun MutableStructure2D<Double>.svdGolubKahanHelper(
g = rv1[i] g = rv1[i]
y = wBuffer[wStart + i] y = wBuffer[wStart + i]
h = s * g h = s * g
g = c * g g *= c
z = pythag(f, h) z = pythag(f, h)
rv1[j] = z rv1[j] = z
c = f / z c = f / z