From ea887b8c728b9a902c80d37a6a092cdf7e7b9a1e Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Wed, 1 Nov 2023 08:55:47 +0300 Subject: [PATCH] 0.4 WIP --- .../space/kscience/kmath/ejml/_generated.kt | 10 ++++- .../kmath/integration/GaussIntegrator.kt | 6 +-- .../kscience/kmath/integration/Integrator.kt | 4 +- .../kmath/integration/SimpsonIntegrator.kt | 12 +++--- .../kmath/integration/SplineIntegrator.kt | 8 ++-- .../kmath/integration/UnivariateIntegrand.kt | 6 +-- .../kmath/geometry/quaternionOperations.kt | 23 +--------- .../kmath/tensors/core/internal/linUtils.kt | 42 ++++++++----------- 8 files changed, 45 insertions(+), 66 deletions(-) diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt index 881649d01..0cef96f49 100644 --- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt +++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt @@ -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.LinearSolverFactory_DSCC import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC -import space.kscience.kmath.UnstableKMathAPI import space.kscience.kmath.linear.* import space.kscience.kmath.linear.Matrix +import space.kscience.kmath.UnstableKMathAPI 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.Float32Field +import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.operations.FloatField 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.FloatBuffer import kotlin.reflect.KClass diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt index 4753edcd0..7d37711cd 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt @@ -54,7 +54,7 @@ public class GaussIntegrator( } } - override fun process(integrand: UnivariateIntegrand): UnivariateIntegrand = with(algebra) { + override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand = with(algebra) { val f = integrand.function val (points, weights) = buildRule(integrand) var res: T = zero @@ -68,7 +68,7 @@ public class GaussIntegrator( res = t } return integrand.withAttributes { - value(res) + IntegrandValue(res) IntegrandCallsPerformed(integrand.calls + points.size) } } @@ -102,7 +102,7 @@ public inline fun GaussIntegrator.integrate( val ranges = UnivariateIntegrandRanges( (0 until intervals).map { i -> (range.start + rangeSize * i)..(range.start + rangeSize * (i + 1)) to order } ) - return process( + return integrate( UnivariateIntegrand( attributeBuilder = { IntegrationRange(range) diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt index 18c46b83b..76fcd0e29 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt @@ -8,9 +8,9 @@ package space.kscience.kmath.integration /** * A general interface for all integrators. */ -public interface Integrator { +public interface Integrator> { /** * 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 } diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SimpsonIntegrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SimpsonIntegrator.kt index c0e4fb394..699d3f0ae 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SimpsonIntegrator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SimpsonIntegrator.kt @@ -44,12 +44,12 @@ public class SimpsonIntegrator( return res } - override fun process(integrand: UnivariateIntegrand): UnivariateIntegrand { + override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand { val ranges = integrand[UnivariateIntegrandRanges] return if (ranges != null) { val res = algebra.sum(ranges.ranges.map { integrateRange(integrand, it.first, it.second) }) integrand.withAttributes { - value(res) + IntegrandValue(res) IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second }) } } else { @@ -58,7 +58,7 @@ public class SimpsonIntegrator( val range = integrand[IntegrationRange] ?: 0.0..1.0 val res = integrateRange(integrand, range, numPoints) integrand.withAttributes { - value(res) + IntegrandValue(res) IntegrandCallsPerformed(integrand.calls + numPoints) } } @@ -96,12 +96,12 @@ public object DoubleSimpsonIntegrator : UnivariateIntegrator { return res } - override fun process(integrand: UnivariateIntegrand): UnivariateIntegrand { + override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand { val ranges = integrand[UnivariateIntegrandRanges] return if (ranges != null) { val res = ranges.ranges.sumOf { integrateRange(integrand, it.first, it.second) } integrand.withAttributes { - value(res) + IntegrandValue(res) IntegrandCallsPerformed(integrand.calls + ranges.ranges.sumOf { it.second }) } } else { @@ -110,7 +110,7 @@ public object DoubleSimpsonIntegrator : UnivariateIntegrator { val range = integrand[IntegrationRange] ?: 0.0..1.0 val res = integrateRange(integrand, range, numPoints) integrand.withAttributes { - value(res) + IntegrandValue(res) IntegrandCallsPerformed(integrand.calls + numPoints) } } diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt index 636ca5a74..c03b248d4 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt @@ -54,7 +54,7 @@ public class SplineIntegrator>( public val algebra: Field, public val bufferFactory: MutableBufferFactory, ) : UnivariateIntegrator { - override fun process(integrand: UnivariateIntegrand): UnivariateIntegrand = algebra { + override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand = algebra { val range = integrand[IntegrationRange] ?: 0.0..1.0 val interpolator: PolynomialInterpolator = SplineInterpolator(algebra, bufferFactory) @@ -72,7 +72,7 @@ public class SplineIntegrator>( ) val res = polynomials.integrate(algebra, number(range.start)..number(range.endInclusive)) integrand.withAttributes { - value(res) + IntegrandValue(res) IntegrandCallsPerformed(integrand.calls + nodes.size) } } @@ -86,7 +86,7 @@ public class SplineIntegrator>( */ @UnstableKMathAPI public object DoubleSplineIntegrator : UnivariateIntegrator { - override fun process(integrand: UnivariateIntegrand): UnivariateIntegrand { + override fun integrate(integrand: UnivariateIntegrand): UnivariateIntegrand { val range = integrand[IntegrationRange] ?: 0.0..1.0 val interpolator: PolynomialInterpolator = SplineInterpolator(Float64Field, ::Float64Buffer) @@ -100,7 +100,7 @@ public object DoubleSplineIntegrator : UnivariateIntegrator { val polynomials = interpolator.interpolatePolynomials(nodes, values) val res = polynomials.integrate(Float64Field, range) return integrand.withAttributes { - value(res) + IntegrandValue(res) IntegrandCallsPerformed(integrand.calls + nodes.size) } } diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt index 6a6d47667..9c39e7edc 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/UnivariateIntegrand.kt @@ -34,7 +34,7 @@ public inline fun UnivariateIntegrand( noinline function: (Double) -> T, ): UnivariateIntegrand = UnivariateIntegrand(safeTypeOf(), Attributes(attributeBuilder), function) -public typealias UnivariateIntegrator = Integrator> +public typealias UnivariateIntegrator = Integrator> public object IntegrationRange : IntegrandAttribute> @@ -70,7 +70,7 @@ public fun TypedAttributesBuilder>.integrationNodes(varar public inline fun UnivariateIntegrator.integrate( attributesBuilder: TypedAttributesBuilder>.() -> Unit, noinline function: (Double) -> T, -): UnivariateIntegrand = process(UnivariateIntegrand(attributesBuilder, function)) +): UnivariateIntegrand = integrate(UnivariateIntegrand(attributesBuilder, function)) /** * A shortcut method to integrate a [function] in [range] with additional features. @@ -83,7 +83,7 @@ public inline fun UnivariateIntegrator.integrate( noinline function: (Double) -> T, ): UnivariateIntegrand { - return process( + return integrate( UnivariateIntegrand( attributeBuilder = { IntegrationRange(range) diff --git a/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/quaternionOperations.kt b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/quaternionOperations.kt index 859255004..4fb3a720c 100644 --- a/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/quaternionOperations.kt +++ b/kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/quaternionOperations.kt @@ -30,25 +30,4 @@ public fun QuaternionAlgebra.angleBetween(q1: Quaternion, q2: Quaternion): Angle /** * 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 - -// -///** -// * 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() -//} \ No newline at end of file +public infix fun Quaternion.dot(other: Quaternion): Double = w * other.w + x * other.x + y * other.y + z * other.z \ No newline at end of file diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt index 272af41d5..66092a532 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt @@ -310,32 +310,26 @@ internal fun DoubleTensorAlgebra.svdHelper( } } -private fun pythag(a: Double, b: Double): Double { - val at: Double = abs(a) - val bt: Double = abs(b) - val ct: Double - val result: Double - if (at > bt) { - ct = bt / at - result = at * sqrt(1.0 + ct * ct) - } else if (bt > 0.0) { - ct = at / bt - result = bt * sqrt(1.0 + ct * ct) - } else result = 0.0 - return result -} - -private fun SIGN(a: Double, b: Double): Double { - if (b >= 0.0) - return abs(a) - else - return -abs(a) -} - internal fun MutableStructure2D.svdGolubKahanHelper( u: MutableStructure2D, w: BufferedTensor, v: MutableStructure2D, iterations: Int, epsilon: Double, ) { + fun pythag(a: Double, b: Double): Double { + val at: Double = abs(a) + val bt: Double = abs(b) + val ct: Double + val result: Double + if (at > bt) { + ct = bt / at + result = at * sqrt(1.0 + ct * ct) + } else if (bt > 0.0) { + ct = at / bt + result = bt * sqrt(1.0 + ct * ct) + } else result = 0.0 + return result + } + + val shape = this.shape val m = shape.component1() val n = shape.component2() @@ -553,7 +547,7 @@ internal fun MutableStructure2D.svdGolubKahanHelper( h = rv1[k] f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y) 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 s = 1.0 @@ -563,7 +557,7 @@ internal fun MutableStructure2D.svdGolubKahanHelper( g = rv1[i] y = wBuffer[wStart + i] h = s * g - g = c * g + g *= c z = pythag(f, h) rv1[j] = z c = f / z