From abae29bbede98231c17ebe00708443ac39fbb854 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sat, 2 Oct 2021 18:54:45 +0300 Subject: [PATCH] DoubleBuffer algebra refactoring --- .../space/kscience/kmath/fit/chiSquared.kt | 11 +- .../kotlin/space/kscience/kmath/fit/qowFit.kt | 4 +- .../kmath/functions/interpolateSquare.kt | 2 +- .../commons/transform/Transformations.kt | 7 +- .../commons/optimization/OptimizeTest.kt | 2 +- .../kmath/expressions/specialExpressions.kt | 2 +- .../kmath/linear/DoubleLinearSpace.kt | 17 +- .../space/kscience/kmath/nd/Structure1D.kt | 2 +- .../kmath/operations/DoubleBufferField.kt | 130 +++ .../operations/DoubleBufferOperations.kt | 150 +-- .../bufferOperation.kt | 3 +- .../kscience/kmath/operations/numbers.kt | 5 +- .../space/kscience/kmath/structures/Buffer.kt | 4 +- .../kmath/streaming/RingBufferTest.kt | 2 +- .../space/kscience/kmath/ejml/_generated.kt | 995 ++++++++++++++++++ .../space/kscience/kmath/real/RealMatrix.kt | 2 +- .../integration/GaussIntegratorRuleFactory.kt | 2 +- .../kmath/integration/SplineIntegrator.kt | 6 +- .../kscience/kmath/geometry/Vector2DTest.kt | 5 +- .../kscience/kmath/geometry/Vector3DTest.kt | 2 +- .../kmath/histogram/UnivariateHistogram.kt | 4 +- .../kscience/kmath/jupyter/KMathJupyter.kt | 2 +- .../space/kscience/kmath/stat/Median.kt | 2 +- .../kmath/tensors/core/internal/linUtils.kt | 2 +- .../kmath/tensors/core/internal/utils.kt | 1 + 25 files changed, 1203 insertions(+), 161 deletions(-) create mode 100644 kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferField.kt rename kmath-core/src/commonMain/kotlin/space/kscience/kmath/{structures => operations}/bufferOperation.kt (97%) create mode 100644 kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt diff --git a/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt b/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt index 8db97d672..dbe0b8454 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt @@ -10,20 +10,21 @@ import kotlinx.html.h3 import space.kscience.kmath.commons.expressions.DSProcessor import space.kscience.kmath.commons.optimization.CMOptimizer import space.kscience.kmath.distributions.NormalDistribution -import space.kscience.kmath.expressions.binding import space.kscience.kmath.expressions.chiSquaredExpression import space.kscience.kmath.expressions.symbol -import space.kscience.kmath.optimization.* +import space.kscience.kmath.operations.asIterable +import space.kscience.kmath.operations.toList +import space.kscience.kmath.optimization.FunctionOptimizationTarget +import space.kscience.kmath.optimization.optimizeWith +import space.kscience.kmath.optimization.resultPoint +import space.kscience.kmath.optimization.resultValue import space.kscience.kmath.real.DoubleVector import space.kscience.kmath.real.map import space.kscience.kmath.real.step import space.kscience.kmath.stat.RandomGenerator -import space.kscience.kmath.structures.asIterable -import space.kscience.kmath.structures.toList import space.kscience.plotly.* import space.kscience.plotly.models.ScatterMode import space.kscience.plotly.models.TraceValues -import kotlin.math.abs import kotlin.math.pow import kotlin.math.sqrt diff --git a/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt b/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt index 04764d763..d52976671 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt @@ -13,6 +13,8 @@ import space.kscience.kmath.distributions.NormalDistribution import space.kscience.kmath.expressions.Symbol import space.kscience.kmath.expressions.binding import space.kscience.kmath.expressions.symbol +import space.kscience.kmath.operations.asIterable +import space.kscience.kmath.operations.toList import space.kscience.kmath.optimization.QowOptimizer import space.kscience.kmath.optimization.chiSquaredOrNull import space.kscience.kmath.optimization.fitWith @@ -20,8 +22,6 @@ import space.kscience.kmath.optimization.resultPoint import space.kscience.kmath.real.map import space.kscience.kmath.real.step import space.kscience.kmath.stat.RandomGenerator -import space.kscience.kmath.structures.asIterable -import space.kscience.kmath.structures.toList import space.kscience.plotly.* import space.kscience.plotly.models.ScatterMode import kotlin.math.abs diff --git a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt b/examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt index 76422d658..3f958b3b0 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt @@ -8,8 +8,8 @@ package space.kscience.kmath.functions import space.kscience.kmath.interpolation.SplineInterpolator import space.kscience.kmath.interpolation.interpolatePolynomials import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.real.map import space.kscience.kmath.real.step -import space.kscience.kmath.structures.map import space.kscience.plotly.Plotly import space.kscience.plotly.UnstablePlotlyAPI import space.kscience.plotly.makeFile diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/transform/Transformations.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/transform/Transformations.kt index 1a99e9fc6..73ab91542 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/transform/Transformations.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/transform/Transformations.kt @@ -10,10 +10,13 @@ import kotlinx.coroutines.flow.Flow import kotlinx.coroutines.flow.map import org.apache.commons.math3.transform.* import space.kscience.kmath.complex.Complex +import space.kscience.kmath.operations.SuspendBufferTransform import space.kscience.kmath.streaming.chunked import space.kscience.kmath.streaming.spread -import space.kscience.kmath.structures.* - +import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.DoubleBuffer +import space.kscience.kmath.structures.VirtualBuffer +import space.kscience.kmath.structures.asBuffer /** diff --git a/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt b/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt index 681ec9ebc..c670ceead 100644 --- a/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt +++ b/kmath-commons/src/test/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt @@ -13,11 +13,11 @@ import space.kscience.kmath.expressions.Symbol.Companion.x import space.kscience.kmath.expressions.Symbol.Companion.y import space.kscience.kmath.expressions.chiSquaredExpression import space.kscience.kmath.expressions.symbol +import space.kscience.kmath.operations.map import space.kscience.kmath.optimization.* import space.kscience.kmath.stat.RandomGenerator import space.kscience.kmath.structures.DoubleBuffer import space.kscience.kmath.structures.asBuffer -import space.kscience.kmath.structures.map import kotlin.math.pow import kotlin.test.Test diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt index 6b17dfca5..907ce4004 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt @@ -6,8 +6,8 @@ package space.kscience.kmath.expressions import space.kscience.kmath.operations.ExtendedField +import space.kscience.kmath.operations.asIterable import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.asIterable import space.kscience.kmath.structures.indices import kotlin.jvm.JvmName diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/DoubleLinearSpace.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/DoubleLinearSpace.kt index a3e031902..4a1311b54 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/DoubleLinearSpace.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/DoubleLinearSpace.kt @@ -9,6 +9,7 @@ import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.nd.DoubleFieldND import space.kscience.kmath.nd.as2D import space.kscience.kmath.nd.asND +import space.kscience.kmath.operations.DoubleBufferOperations import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.DoubleBuffer @@ -89,17 +90,21 @@ public object DoubleLinearSpace : LinearSpace { asND().map { it * value }.as2D() } - public override fun Point.plus(other: Point): DoubleBuffer = DoubleBuffer(size) { - get(it) + other[it] + public override fun Point.plus(other: Point): DoubleBuffer = DoubleBufferOperations.run { + this@plus + other } - public override fun Point.minus(other: Point): DoubleBuffer = DoubleBuffer(size) { - get(it) - other[it] + public override fun Point.minus(other: Point): DoubleBuffer = DoubleBufferOperations.run { + this@minus - other } - public override fun Point.times(value: Double): DoubleBuffer = DoubleBuffer(size) { i -> get(i) * value } + public override fun Point.times(value: Double): DoubleBuffer = DoubleBufferOperations.run { + scale(this@times, value) + } - public operator fun Point.div(value: Double): DoubleBuffer = DoubleBuffer(size) { i -> get(i) / value } + public operator fun Point.div(value: Double): DoubleBuffer = DoubleBufferOperations.run { + scale(this@div, 1.0 / value) + } public override fun Double.times(v: Point): DoubleBuffer = v * this diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt index f2a4336eb..3dcc77334 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt @@ -6,10 +6,10 @@ package space.kscience.kmath.nd import space.kscience.kmath.misc.PerformancePitfall +import space.kscience.kmath.operations.asSequence import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.MutableBuffer import space.kscience.kmath.structures.asMutableBuffer -import space.kscience.kmath.structures.asSequence import kotlin.jvm.JvmInline /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferField.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferField.kt new file mode 100644 index 000000000..acc2c2dc0 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferField.kt @@ -0,0 +1,130 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file. + */ + +package space.kscience.kmath.operations + +import space.kscience.kmath.structures.Buffer +import space.kscience.kmath.structures.DoubleBuffer + +/** + * [ExtendedField] over [DoubleBuffer]. + * + * @property size the size of buffers to operate on. + */ +public class DoubleBufferField(public val size: Int) : ExtendedField>, DoubleBufferOperations() { + override val zero: Buffer by lazy { DoubleBuffer(size) { 0.0 } } + override val one: Buffer by lazy { DoubleBuffer(size) { 1.0 } } + + override fun sinh(arg: Buffer): DoubleBuffer = super.sinh(arg) + + override fun cosh(arg: Buffer): DoubleBuffer = super.cosh(arg) + + override fun tanh(arg: Buffer): DoubleBuffer = super.tanh(arg) + + override fun asinh(arg: Buffer): DoubleBuffer = super.asinh(arg) + + override fun acosh(arg: Buffer): DoubleBuffer = super.acosh(arg) + + override fun atanh(arg: Buffer): DoubleBuffer= super.atanh(arg) + + // override fun number(value: Number): Buffer = DoubleBuffer(size) { value.toDouble() } +// +// override fun Buffer.unaryMinus(): Buffer = DoubleBufferOperations.run { +// -this@unaryMinus +// } +// +// override fun add(a: Buffer, b: Buffer): DoubleBuffer { +// require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } +// return DoubleBufferOperations.add(a, b) +// } +// + +// +// override fun multiply(a: Buffer, b: Buffer): DoubleBuffer { +// require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } +// return DoubleBufferOperations.multiply(a, b) +// } +// +// override fun divide(a: Buffer, b: Buffer): DoubleBuffer { +// require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } +// return DoubleBufferOperations.divide(a, b) +// } +// +// override fun sin(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.sin(arg) +// } +// +// override fun cos(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.cos(arg) +// } +// +// override fun tan(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.tan(arg) +// } +// +// override fun asin(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.asin(arg) +// } +// +// override fun acos(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.acos(arg) +// } +// +// override fun atan(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.atan(arg) +// } +// +// override fun sinh(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.sinh(arg) +// } +// +// override fun cosh(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.cosh(arg) +// } +// +// override fun tanh(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.tanh(arg) +// } +// +// override fun asinh(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.asinh(arg) +// } +// +// override fun acosh(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.acosh(arg) +// } +// +// override fun atanh(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.atanh(arg) +// } +// +// override fun power(arg: Buffer, pow: Number): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.power(arg, pow) +// } +// +// override fun exp(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.exp(arg) +// } +// +// override fun ln(arg: Buffer): DoubleBuffer { +// require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } +// return DoubleBufferOperations.ln(arg) +// } + +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferOperations.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferOperations.kt index a95bd2541..50b821962 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferOperations.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferOperations.kt @@ -8,13 +8,13 @@ package space.kscience.kmath.operations import space.kscience.kmath.linear.Point import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.DoubleBuffer -import space.kscience.kmath.structures.fold + import kotlin.math.* /** * [ExtendedFieldOperations] over [DoubleBuffer]. */ -public object DoubleBufferOperations : ExtendedFieldOperations> { +public abstract class DoubleBufferOperations : ExtendedFieldOperations>, Norm, Double> { override fun Buffer.unaryMinus(): DoubleBuffer = if (this is DoubleBuffer) { DoubleBuffer(size) { -array[it] } } else { @@ -32,7 +32,22 @@ public object DoubleBufferOperations : ExtendedFieldOperations> { DoubleBuffer(DoubleArray(a.size) { aArray[it] + bArray[it] }) } else DoubleBuffer(DoubleArray(a.size) { a[it] + b[it] }) } -// + + override fun Buffer.plus(b: Buffer): DoubleBuffer = add(this, b) + + override fun Buffer.minus(b: Buffer): DoubleBuffer { + require(b.size == this.size) { + "The size of the first buffer ${this.size} should be the same as for second one: ${b.size} " + } + + return if (this is DoubleBuffer && b is DoubleBuffer) { + val aArray = this.array + val bArray = b.array + DoubleBuffer(DoubleArray(this.size) { aArray[it] - bArray[it] }) + } else DoubleBuffer(DoubleArray(this.size) { this[it] - b[it] }) + } + + // // override fun multiply(a: Buffer, k: Number): RealBuffer { // val kValue = k.toDouble() // @@ -159,128 +174,21 @@ public object DoubleBufferOperations : ExtendedFieldOperations> { override fun ln(arg: Buffer): DoubleBuffer = if (arg is DoubleBuffer) { val array = arg.array DoubleBuffer(DoubleArray(arg.size) { ln(array[it]) }) - } else + } else { DoubleBuffer(DoubleArray(arg.size) { ln(arg[it]) }) + } + + override fun norm(arg: Buffer): Double = DoubleL2Norm.norm(arg) + + override fun scale(a: Buffer, value: Double): DoubleBuffer = if (a is DoubleBuffer) { + val aArray = a.array + DoubleBuffer(DoubleArray(a.size) { aArray[it] * value }) + } else DoubleBuffer(DoubleArray(a.size) { a[it] * value }) + + public companion object : DoubleBufferOperations() } public object DoubleL2Norm : Norm, Double> { override fun norm(arg: Point): Double = sqrt(arg.fold(0.0) { acc: Double, d: Double -> acc + d.pow(2) }) } -/** - * [ExtendedField] over [DoubleBuffer]. - * - * @property size the size of buffers to operate on. - */ -@Deprecated("To be replaced by generic BufferAlgebra") -public class DoubleBufferField(public val size: Int) : ExtendedField>, Norm, Double> { - override val zero: Buffer by lazy { DoubleBuffer(size) { 0.0 } } - override val one: Buffer by lazy { DoubleBuffer(size) { 1.0 } } - - override fun number(value: Number): Buffer = DoubleBuffer(size) { value.toDouble() } - - override fun Buffer.unaryMinus(): Buffer = DoubleBufferOperations.run { - -this@unaryMinus - } - - override fun add(a: Buffer, b: Buffer): DoubleBuffer { - require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } - return DoubleBufferOperations.add(a, b) - } - - override fun scale(a: Buffer, value: Double): DoubleBuffer { - require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } - - return if (a is DoubleBuffer) { - val aArray = a.array - DoubleBuffer(DoubleArray(a.size) { aArray[it] * value }) - } else DoubleBuffer(DoubleArray(a.size) { a[it] * value }) - } - - override fun multiply(a: Buffer, b: Buffer): DoubleBuffer { - require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } - return DoubleBufferOperations.multiply(a, b) - } - - override fun divide(a: Buffer, b: Buffer): DoubleBuffer { - require(a.size == size) { "The buffer size ${a.size} does not match context size $size" } - return DoubleBufferOperations.divide(a, b) - } - - override fun sin(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.sin(arg) - } - - override fun cos(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.cos(arg) - } - - override fun tan(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.tan(arg) - } - - override fun asin(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.asin(arg) - } - - override fun acos(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.acos(arg) - } - - override fun atan(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.atan(arg) - } - - override fun sinh(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.sinh(arg) - } - - override fun cosh(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.cosh(arg) - } - - override fun tanh(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.tanh(arg) - } - - override fun asinh(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.asinh(arg) - } - - override fun acosh(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.acosh(arg) - } - - override fun atanh(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.atanh(arg) - } - - override fun power(arg: Buffer, pow: Number): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.power(arg, pow) - } - - override fun exp(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.exp(arg) - } - - override fun ln(arg: Buffer): DoubleBuffer { - require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" } - return DoubleBufferOperations.ln(arg) - } - - override fun norm(arg: Buffer): Double = DoubleL2Norm.norm(arg) -} diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferOperation.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/bufferOperation.kt similarity index 97% rename from kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferOperation.kt rename to kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/bufferOperation.kt index bc3f5b64b..6bf3266e3 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferOperation.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/bufferOperation.kt @@ -3,9 +3,10 @@ * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file. */ -package space.kscience.kmath.structures +package space.kscience.kmath.operations import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.structures.* /** * Typealias for buffer transformations. diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt index c90553aaa..4c0010bf9 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/numbers.kt @@ -14,7 +14,8 @@ public interface ExtendedFieldOperations : FieldOperations, TrigonometricOperations, PowerOperations, - ExponentialOperations { + ExponentialOperations, + ScaleOperations { override fun tan(arg: T): T = sin(arg) / cos(arg) override fun tanh(arg: T): T = sinh(arg) / cosh(arg) @@ -41,7 +42,7 @@ public interface ExtendedFieldOperations : /** * Advanced Number-like field that implements basic operations. */ -public interface ExtendedField : ExtendedFieldOperations, Field, NumericAlgebra, ScaleOperations { +public interface ExtendedField : ExtendedFieldOperations, Field, NumericAlgebra{ override fun sinh(arg: T): T = (exp(arg) - exp(-arg)) / 2.0 override fun cosh(arg: T): T = (exp(arg) + exp(-arg)) / 2.0 override fun tanh(arg: T): T = (exp(arg) - exp(-arg)) / (exp(-arg) + exp(arg)) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt index fc23169ca..c68bca2d9 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt @@ -5,6 +5,7 @@ package space.kscience.kmath.structures +import space.kscience.kmath.operations.asSequence import kotlin.jvm.JvmInline import kotlin.reflect.KClass @@ -49,7 +50,8 @@ public interface Buffer { public companion object { - public fun toString(buffer: Buffer<*>): String = buffer.asSequence().joinToString(prefix = "[", separator = ", ", postfix = "]") + public fun toString(buffer: Buffer<*>): String = + buffer.asSequence().joinToString(prefix = "[", separator = ", ", postfix = "]") /** * Check the element-by-element match of content of two buffers. diff --git a/kmath-coroutines/src/jvmTest/kotlin/space/kscience/kmath/streaming/RingBufferTest.kt b/kmath-coroutines/src/jvmTest/kotlin/space/kscience/kmath/streaming/RingBufferTest.kt index 14081b0f5..a3143a1ac 100644 --- a/kmath-coroutines/src/jvmTest/kotlin/space/kscience/kmath/streaming/RingBufferTest.kt +++ b/kmath-coroutines/src/jvmTest/kotlin/space/kscience/kmath/streaming/RingBufferTest.kt @@ -7,7 +7,7 @@ package space.kscience.kmath.streaming import kotlinx.coroutines.flow.* import kotlinx.coroutines.runBlocking -import space.kscience.kmath.structures.asSequence +import space.kscience.kmath.operations.asSequence import kotlin.test.Test import kotlin.test.assertEquals 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 new file mode 100644 index 000000000..dce739dc2 --- /dev/null +++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt @@ -0,0 +1,995 @@ +/* + * Copyright 2018-2021 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file. + */ + +/* This file is generated with buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt */ + +package space.kscience.kmath.ejml + +import org.ejml.data.* +import org.ejml.dense.row.CommonOps_DDRM +import org.ejml.dense.row.CommonOps_FDRM +import org.ejml.dense.row.factory.DecompositionFactory_DDRM +import org.ejml.dense.row.factory.DecompositionFactory_FDRM +import org.ejml.sparse.FillReducing +import org.ejml.sparse.csc.CommonOps_DSCC +import org.ejml.sparse.csc.CommonOps_FSCC +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.linear.* +import space.kscience.kmath.linear.Matrix +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.nd.StructureFeature +import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.operations.FloatField +import space.kscience.kmath.operations.invoke +import space.kscience.kmath.structures.DoubleBuffer +import space.kscience.kmath.structures.FloatBuffer +import kotlin.reflect.KClass +import kotlin.reflect.cast + +/** + * [EjmlVector] specialization for [Double]. + */ +public class EjmlDoubleVector(override val origin: M) : EjmlVector(origin) { + init { + require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" } + } + + override operator fun get(index: Int): Double = origin[0, index] +} + +/** + * [EjmlVector] specialization for [Float]. + */ +public class EjmlFloatVector(override val origin: M) : EjmlVector(origin) { + init { + require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" } + } + + override operator fun get(index: Int): Float = origin[0, index] +} + +/** + * [EjmlMatrix] specialization for [Double]. + */ +public class EjmlDoubleMatrix(override val origin: M) : EjmlMatrix(origin) { + override operator fun get(i: Int, j: Int): Double = origin[i, j] +} + +/** + * [EjmlMatrix] specialization for [Float]. + */ +public class EjmlFloatMatrix(override val origin: M) : EjmlMatrix(origin) { + override operator fun get(i: Int, j: Int): Float = origin[i, j] +} + +/** + * [EjmlLinearSpace] implementation based on [CommonOps_DDRM], [DecompositionFactory_DDRM] operations and + * [DMatrixRMaj] matrices. + */ +public object EjmlLinearSpaceDDRM : EjmlLinearSpace() { + /** + * The [DoubleField] reference. + */ + override val elementAlgebra: DoubleField get() = DoubleField + + @Suppress("UNCHECKED_CAST") + override fun Matrix.toEjml(): EjmlDoubleMatrix = when { + this is EjmlDoubleMatrix<*> && origin is DMatrixRMaj -> this as EjmlDoubleMatrix + else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) } + } + + @Suppress("UNCHECKED_CAST") + override fun Point.toEjml(): EjmlDoubleVector = when { + this is EjmlDoubleVector<*> && origin is DMatrixRMaj -> this as EjmlDoubleVector + else -> EjmlDoubleVector(DMatrixRMaj(size, 1).also { + (0 until it.numRows).forEach { row -> it[row, 0] = get(row) } + }) + } + + override fun buildMatrix( + rows: Int, + columns: Int, + initializer: DoubleField.(i: Int, j: Int) -> Double, + ): EjmlDoubleMatrix = DMatrixRMaj(rows, columns).also { + (0 until rows).forEach { row -> + (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) } + } + }.wrapMatrix() + + override fun buildVector( + size: Int, + initializer: DoubleField.(Int) -> Double, + ): EjmlDoubleVector = EjmlDoubleVector(DMatrixRMaj(size, 1).also { + (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) } + }) + + private fun T.wrapMatrix() = EjmlDoubleMatrix(this) + private fun T.wrapVector() = EjmlDoubleVector(this) + + override fun Matrix.unaryMinus(): Matrix = this * elementAlgebra { -one } + + override fun Matrix.dot(other: Matrix): EjmlDoubleMatrix { + val out = DMatrixRMaj(1, 1) + CommonOps_DDRM.mult(toEjml().origin, other.toEjml().origin, out) + return out.wrapMatrix() + } + + override fun Matrix.dot(vector: Point): EjmlDoubleVector { + val out = DMatrixRMaj(1, 1) + CommonOps_DDRM.mult(toEjml().origin, vector.toEjml().origin, out) + return out.wrapVector() + } + + override operator fun Matrix.minus(other: Matrix): EjmlDoubleMatrix { + val out = DMatrixRMaj(1, 1) + + CommonOps_DDRM.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra { -one }, + other.toEjml().origin, + out, + ) + + return out.wrapMatrix() + } + + override operator fun Matrix.times(value: Double): EjmlDoubleMatrix { + val res = DMatrixRMaj(1, 1) + CommonOps_DDRM.scale(value, toEjml().origin, res) + return res.wrapMatrix() + } + + override fun Point.unaryMinus(): EjmlDoubleVector { + val res = DMatrixRMaj(1, 1) + CommonOps_DDRM.changeSign(toEjml().origin, res) + return res.wrapVector() + } + + override fun Matrix.plus(other: Matrix): EjmlDoubleMatrix { + val out = DMatrixRMaj(1, 1) + + CommonOps_DDRM.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra.one, + other.toEjml().origin, + out, + ) + + return out.wrapMatrix() + } + + override fun Point.plus(other: Point): EjmlDoubleVector { + val out = DMatrixRMaj(1, 1) + + CommonOps_DDRM.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra.one, + other.toEjml().origin, + out, + ) + + return out.wrapVector() + } + + override fun Point.minus(other: Point): EjmlDoubleVector { + val out = DMatrixRMaj(1, 1) + + CommonOps_DDRM.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra { -one }, + other.toEjml().origin, + out, + ) + + return out.wrapVector() + } + + override fun Double.times(m: Matrix): EjmlDoubleMatrix = m * this + + override fun Point.times(value: Double): EjmlDoubleVector { + val res = DMatrixRMaj(1, 1) + CommonOps_DDRM.scale(value, toEjml().origin, res) + return res.wrapVector() + } + + override fun Double.times(v: Point): EjmlDoubleVector = v * this + + @UnstableKMathAPI + override fun computeFeature(structure: Matrix, type: KClass): F? { + structure.getFeature(type)?.let { return it } + val origin = structure.toEjml().origin + + return when (type) { + InverseMatrixFeature::class -> object : InverseMatrixFeature { + override val inverse: Matrix by lazy { + val res = origin.copy() + CommonOps_DDRM.invert(res) + res.wrapMatrix() + } + } + + DeterminantFeature::class -> object : DeterminantFeature { + override val determinant: Double by lazy { CommonOps_DDRM.det(origin) } + } + + SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature { + private val svd by lazy { + DecompositionFactory_DDRM.svd(origin.numRows, origin.numCols, true, true, false) + .apply { decompose(origin.copy()) } + } + + override val u: Matrix by lazy { svd.getU(null, false).wrapMatrix() } + override val s: Matrix by lazy { svd.getW(null).wrapMatrix() } + override val v: Matrix by lazy { svd.getV(null, false).wrapMatrix() } + override val singularValues: Point by lazy { DoubleBuffer(svd.singularValues) } + } + + QRDecompositionFeature::class -> object : QRDecompositionFeature { + private val qr by lazy { + DecompositionFactory_DDRM.qr().apply { decompose(origin.copy()) } + } + + override val q: Matrix by lazy { + qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature) + } + + override val r: Matrix by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) } + } + + CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature { + override val l: Matrix by lazy { + val cholesky = + DecompositionFactory_DDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) } + + cholesky.getT(null).wrapMatrix().withFeature(LFeature) + } + } + + LupDecompositionFeature::class -> object : LupDecompositionFeature { + private val lup by lazy { + DecompositionFactory_DDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) } + } + + override val l: Matrix by lazy { + lup.getLower(null).wrapMatrix().withFeature(LFeature) + } + + override val u: Matrix by lazy { + lup.getUpper(null).wrapMatrix().withFeature(UFeature) + } + + override val p: Matrix by lazy { lup.getRowPivot(null).wrapMatrix() } + } + + else -> null + }?.let(type::cast) + } + + /** + * Solves for *x* in the following equation: *x = [a] -1 · [b]*. + * + * @param a the base matrix. + * @param b n by p matrix. + * @return the solution for *x* that is n by p. + */ + public fun solve(a: Matrix, b: Matrix): EjmlDoubleMatrix { + val res = DMatrixRMaj(1, 1) + CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res) + return res.wrapMatrix() + } + + /** + * Solves for *x* in the following equation: *x = [a] -1 · [b]*. + * + * @param a the base matrix. + * @param b n by p vector. + * @return the solution for *x* that is n by p. + */ + public fun solve(a: Matrix, b: Point): EjmlDoubleVector { + val res = DMatrixRMaj(1, 1) + CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res) + return EjmlDoubleVector(res) + } +} + +/** + * [EjmlLinearSpace] implementation based on [CommonOps_FDRM], [DecompositionFactory_FDRM] operations and + * [FMatrixRMaj] matrices. + */ +public object EjmlLinearSpaceFDRM : EjmlLinearSpace() { + /** + * The [FloatField] reference. + */ + override val elementAlgebra: FloatField get() = FloatField + + @Suppress("UNCHECKED_CAST") + override fun Matrix.toEjml(): EjmlFloatMatrix = when { + this is EjmlFloatMatrix<*> && origin is FMatrixRMaj -> this as EjmlFloatMatrix + else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) } + } + + @Suppress("UNCHECKED_CAST") + override fun Point.toEjml(): EjmlFloatVector = when { + this is EjmlFloatVector<*> && origin is FMatrixRMaj -> this as EjmlFloatVector + else -> EjmlFloatVector(FMatrixRMaj(size, 1).also { + (0 until it.numRows).forEach { row -> it[row, 0] = get(row) } + }) + } + + override fun buildMatrix( + rows: Int, + columns: Int, + initializer: FloatField.(i: Int, j: Int) -> Float, + ): EjmlFloatMatrix = FMatrixRMaj(rows, columns).also { + (0 until rows).forEach { row -> + (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) } + } + }.wrapMatrix() + + override fun buildVector( + size: Int, + initializer: FloatField.(Int) -> Float, + ): EjmlFloatVector = EjmlFloatVector(FMatrixRMaj(size, 1).also { + (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) } + }) + + private fun T.wrapMatrix() = EjmlFloatMatrix(this) + private fun T.wrapVector() = EjmlFloatVector(this) + + override fun Matrix.unaryMinus(): Matrix = this * elementAlgebra { -one } + + override fun Matrix.dot(other: Matrix): EjmlFloatMatrix { + val out = FMatrixRMaj(1, 1) + CommonOps_FDRM.mult(toEjml().origin, other.toEjml().origin, out) + return out.wrapMatrix() + } + + override fun Matrix.dot(vector: Point): EjmlFloatVector { + val out = FMatrixRMaj(1, 1) + CommonOps_FDRM.mult(toEjml().origin, vector.toEjml().origin, out) + return out.wrapVector() + } + + override operator fun Matrix.minus(other: Matrix): EjmlFloatMatrix { + val out = FMatrixRMaj(1, 1) + + CommonOps_FDRM.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra { -one }, + other.toEjml().origin, + out, + ) + + return out.wrapMatrix() + } + + override operator fun Matrix.times(value: Float): EjmlFloatMatrix { + val res = FMatrixRMaj(1, 1) + CommonOps_FDRM.scale(value, toEjml().origin, res) + return res.wrapMatrix() + } + + override fun Point.unaryMinus(): EjmlFloatVector { + val res = FMatrixRMaj(1, 1) + CommonOps_FDRM.changeSign(toEjml().origin, res) + return res.wrapVector() + } + + override fun Matrix.plus(other: Matrix): EjmlFloatMatrix { + val out = FMatrixRMaj(1, 1) + + CommonOps_FDRM.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra.one, + other.toEjml().origin, + out, + ) + + return out.wrapMatrix() + } + + override fun Point.plus(other: Point): EjmlFloatVector { + val out = FMatrixRMaj(1, 1) + + CommonOps_FDRM.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra.one, + other.toEjml().origin, + out, + ) + + return out.wrapVector() + } + + override fun Point.minus(other: Point): EjmlFloatVector { + val out = FMatrixRMaj(1, 1) + + CommonOps_FDRM.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra { -one }, + other.toEjml().origin, + out, + ) + + return out.wrapVector() + } + + override fun Float.times(m: Matrix): EjmlFloatMatrix = m * this + + override fun Point.times(value: Float): EjmlFloatVector { + val res = FMatrixRMaj(1, 1) + CommonOps_FDRM.scale(value, toEjml().origin, res) + return res.wrapVector() + } + + override fun Float.times(v: Point): EjmlFloatVector = v * this + + @UnstableKMathAPI + override fun computeFeature(structure: Matrix, type: KClass): F? { + structure.getFeature(type)?.let { return it } + val origin = structure.toEjml().origin + + return when (type) { + InverseMatrixFeature::class -> object : InverseMatrixFeature { + override val inverse: Matrix by lazy { + val res = origin.copy() + CommonOps_FDRM.invert(res) + res.wrapMatrix() + } + } + + DeterminantFeature::class -> object : DeterminantFeature { + override val determinant: Float by lazy { CommonOps_FDRM.det(origin) } + } + + SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature { + private val svd by lazy { + DecompositionFactory_FDRM.svd(origin.numRows, origin.numCols, true, true, false) + .apply { decompose(origin.copy()) } + } + + override val u: Matrix by lazy { svd.getU(null, false).wrapMatrix() } + override val s: Matrix by lazy { svd.getW(null).wrapMatrix() } + override val v: Matrix by lazy { svd.getV(null, false).wrapMatrix() } + override val singularValues: Point by lazy { FloatBuffer(svd.singularValues) } + } + + QRDecompositionFeature::class -> object : QRDecompositionFeature { + private val qr by lazy { + DecompositionFactory_FDRM.qr().apply { decompose(origin.copy()) } + } + + override val q: Matrix by lazy { + qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature) + } + + override val r: Matrix by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) } + } + + CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature { + override val l: Matrix by lazy { + val cholesky = + DecompositionFactory_FDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) } + + cholesky.getT(null).wrapMatrix().withFeature(LFeature) + } + } + + LupDecompositionFeature::class -> object : LupDecompositionFeature { + private val lup by lazy { + DecompositionFactory_FDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) } + } + + override val l: Matrix by lazy { + lup.getLower(null).wrapMatrix().withFeature(LFeature) + } + + override val u: Matrix by lazy { + lup.getUpper(null).wrapMatrix().withFeature(UFeature) + } + + override val p: Matrix by lazy { lup.getRowPivot(null).wrapMatrix() } + } + + else -> null + }?.let(type::cast) + } + + /** + * Solves for *x* in the following equation: *x = [a] -1 · [b]*. + * + * @param a the base matrix. + * @param b n by p matrix. + * @return the solution for *x* that is n by p. + */ + public fun solve(a: Matrix, b: Matrix): EjmlFloatMatrix { + val res = FMatrixRMaj(1, 1) + CommonOps_FDRM.solve(FMatrixRMaj(a.toEjml().origin), FMatrixRMaj(b.toEjml().origin), res) + return res.wrapMatrix() + } + + /** + * Solves for *x* in the following equation: *x = [a] -1 · [b]*. + * + * @param a the base matrix. + * @param b n by p vector. + * @return the solution for *x* that is n by p. + */ + public fun solve(a: Matrix, b: Point): EjmlFloatVector { + val res = FMatrixRMaj(1, 1) + CommonOps_FDRM.solve(FMatrixRMaj(a.toEjml().origin), FMatrixRMaj(b.toEjml().origin), res) + return EjmlFloatVector(res) + } +} + +/** + * [EjmlLinearSpace] implementation based on [CommonOps_DSCC], [DecompositionFactory_DSCC] operations and + * [DMatrixSparseCSC] matrices. + */ +public object EjmlLinearSpaceDSCC : EjmlLinearSpace() { + /** + * The [DoubleField] reference. + */ + override val elementAlgebra: DoubleField get() = DoubleField + + @Suppress("UNCHECKED_CAST") + override fun Matrix.toEjml(): EjmlDoubleMatrix = when { + this is EjmlDoubleMatrix<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleMatrix + else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) } + } + + @Suppress("UNCHECKED_CAST") + override fun Point.toEjml(): EjmlDoubleVector = when { + this is EjmlDoubleVector<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleVector + else -> EjmlDoubleVector(DMatrixSparseCSC(size, 1).also { + (0 until it.numRows).forEach { row -> it[row, 0] = get(row) } + }) + } + + override fun buildMatrix( + rows: Int, + columns: Int, + initializer: DoubleField.(i: Int, j: Int) -> Double, + ): EjmlDoubleMatrix = DMatrixSparseCSC(rows, columns).also { + (0 until rows).forEach { row -> + (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) } + } + }.wrapMatrix() + + override fun buildVector( + size: Int, + initializer: DoubleField.(Int) -> Double, + ): EjmlDoubleVector = EjmlDoubleVector(DMatrixSparseCSC(size, 1).also { + (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) } + }) + + private fun T.wrapMatrix() = EjmlDoubleMatrix(this) + private fun T.wrapVector() = EjmlDoubleVector(this) + + override fun Matrix.unaryMinus(): Matrix = this * elementAlgebra { -one } + + override fun Matrix.dot(other: Matrix): EjmlDoubleMatrix { + val out = DMatrixSparseCSC(1, 1) + CommonOps_DSCC.mult(toEjml().origin, other.toEjml().origin, out) + return out.wrapMatrix() + } + + override fun Matrix.dot(vector: Point): EjmlDoubleVector { + val out = DMatrixSparseCSC(1, 1) + CommonOps_DSCC.mult(toEjml().origin, vector.toEjml().origin, out) + return out.wrapVector() + } + + override operator fun Matrix.minus(other: Matrix): EjmlDoubleMatrix { + val out = DMatrixSparseCSC(1, 1) + + CommonOps_DSCC.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra { -one }, + other.toEjml().origin, + out, + null, + null, + ) + + return out.wrapMatrix() + } + + override operator fun Matrix.times(value: Double): EjmlDoubleMatrix { + val res = DMatrixSparseCSC(1, 1) + CommonOps_DSCC.scale(value, toEjml().origin, res) + return res.wrapMatrix() + } + + override fun Point.unaryMinus(): EjmlDoubleVector { + val res = DMatrixSparseCSC(1, 1) + CommonOps_DSCC.changeSign(toEjml().origin, res) + return res.wrapVector() + } + + override fun Matrix.plus(other: Matrix): EjmlDoubleMatrix { + val out = DMatrixSparseCSC(1, 1) + + CommonOps_DSCC.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra.one, + other.toEjml().origin, + out, + null, + null, + ) + + return out.wrapMatrix() + } + + override fun Point.plus(other: Point): EjmlDoubleVector { + val out = DMatrixSparseCSC(1, 1) + + CommonOps_DSCC.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra.one, + other.toEjml().origin, + out, + null, + null, + ) + + return out.wrapVector() + } + + override fun Point.minus(other: Point): EjmlDoubleVector { + val out = DMatrixSparseCSC(1, 1) + + CommonOps_DSCC.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra { -one }, + other.toEjml().origin, + out, + null, + null, + ) + + return out.wrapVector() + } + + override fun Double.times(m: Matrix): EjmlDoubleMatrix = m * this + + override fun Point.times(value: Double): EjmlDoubleVector { + val res = DMatrixSparseCSC(1, 1) + CommonOps_DSCC.scale(value, toEjml().origin, res) + return res.wrapVector() + } + + override fun Double.times(v: Point): EjmlDoubleVector = v * this + + @UnstableKMathAPI + override fun computeFeature(structure: Matrix, type: KClass): F? { + structure.getFeature(type)?.let { return it } + val origin = structure.toEjml().origin + + return when (type) { + QRDecompositionFeature::class -> object : QRDecompositionFeature { + private val qr by lazy { + DecompositionFactory_DSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) } + } + + override val q: Matrix by lazy { + qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature) + } + + override val r: Matrix by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) } + } + + CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature { + override val l: Matrix by lazy { + val cholesky = + DecompositionFactory_DSCC.cholesky().apply { decompose(origin.copy()) } + + (cholesky.getT(null) as DMatrix).wrapMatrix().withFeature(LFeature) + } + } + + LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object : + LUDecompositionFeature, DeterminantFeature, InverseMatrixFeature { + private val lu by lazy { + DecompositionFactory_DSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) } + } + + override val l: Matrix by lazy { + lu.getLower(null).wrapMatrix().withFeature(LFeature) + } + + override val u: Matrix by lazy { + lu.getUpper(null).wrapMatrix().withFeature(UFeature) + } + + override val inverse: Matrix by lazy { + var a = origin + val inverse = DMatrixRMaj(1, 1) + val solver = LinearSolverFactory_DSCC.lu(FillReducing.NONE) + if (solver.modifiesA()) a = a.copy() + val i = CommonOps_DDRM.identity(a.numRows) + solver.solve(i, inverse) + inverse.wrapMatrix() + } + + override val determinant: Double by lazy { elementAlgebra.number(lu.computeDeterminant().real) } + } + + else -> null + }?.let(type::cast) + } + + /** + * Solves for *x* in the following equation: *x = [a] -1 · [b]*. + * + * @param a the base matrix. + * @param b n by p matrix. + * @return the solution for *x* that is n by p. + */ + public fun solve(a: Matrix, b: Matrix): EjmlDoubleMatrix { + val res = DMatrixSparseCSC(1, 1) + CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res) + return res.wrapMatrix() + } + + /** + * Solves for *x* in the following equation: *x = [a] -1 · [b]*. + * + * @param a the base matrix. + * @param b n by p vector. + * @return the solution for *x* that is n by p. + */ + public fun solve(a: Matrix, b: Point): EjmlDoubleVector { + val res = DMatrixSparseCSC(1, 1) + CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res) + return EjmlDoubleVector(res) + } +} + +/** + * [EjmlLinearSpace] implementation based on [CommonOps_FSCC], [DecompositionFactory_FSCC] operations and + * [FMatrixSparseCSC] matrices. + */ +public object EjmlLinearSpaceFSCC : EjmlLinearSpace() { + /** + * The [FloatField] reference. + */ + override val elementAlgebra: FloatField get() = FloatField + + @Suppress("UNCHECKED_CAST") + override fun Matrix.toEjml(): EjmlFloatMatrix = when { + this is EjmlFloatMatrix<*> && origin is FMatrixSparseCSC -> this as EjmlFloatMatrix + else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) } + } + + @Suppress("UNCHECKED_CAST") + override fun Point.toEjml(): EjmlFloatVector = when { + this is EjmlFloatVector<*> && origin is FMatrixSparseCSC -> this as EjmlFloatVector + else -> EjmlFloatVector(FMatrixSparseCSC(size, 1).also { + (0 until it.numRows).forEach { row -> it[row, 0] = get(row) } + }) + } + + override fun buildMatrix( + rows: Int, + columns: Int, + initializer: FloatField.(i: Int, j: Int) -> Float, + ): EjmlFloatMatrix = FMatrixSparseCSC(rows, columns).also { + (0 until rows).forEach { row -> + (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) } + } + }.wrapMatrix() + + override fun buildVector( + size: Int, + initializer: FloatField.(Int) -> Float, + ): EjmlFloatVector = EjmlFloatVector(FMatrixSparseCSC(size, 1).also { + (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) } + }) + + private fun T.wrapMatrix() = EjmlFloatMatrix(this) + private fun T.wrapVector() = EjmlFloatVector(this) + + override fun Matrix.unaryMinus(): Matrix = this * elementAlgebra { -one } + + override fun Matrix.dot(other: Matrix): EjmlFloatMatrix { + val out = FMatrixSparseCSC(1, 1) + CommonOps_FSCC.mult(toEjml().origin, other.toEjml().origin, out) + return out.wrapMatrix() + } + + override fun Matrix.dot(vector: Point): EjmlFloatVector { + val out = FMatrixSparseCSC(1, 1) + CommonOps_FSCC.mult(toEjml().origin, vector.toEjml().origin, out) + return out.wrapVector() + } + + override operator fun Matrix.minus(other: Matrix): EjmlFloatMatrix { + val out = FMatrixSparseCSC(1, 1) + + CommonOps_FSCC.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra { -one }, + other.toEjml().origin, + out, + null, + null, + ) + + return out.wrapMatrix() + } + + override operator fun Matrix.times(value: Float): EjmlFloatMatrix { + val res = FMatrixSparseCSC(1, 1) + CommonOps_FSCC.scale(value, toEjml().origin, res) + return res.wrapMatrix() + } + + override fun Point.unaryMinus(): EjmlFloatVector { + val res = FMatrixSparseCSC(1, 1) + CommonOps_FSCC.changeSign(toEjml().origin, res) + return res.wrapVector() + } + + override fun Matrix.plus(other: Matrix): EjmlFloatMatrix { + val out = FMatrixSparseCSC(1, 1) + + CommonOps_FSCC.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra.one, + other.toEjml().origin, + out, + null, + null, + ) + + return out.wrapMatrix() + } + + override fun Point.plus(other: Point): EjmlFloatVector { + val out = FMatrixSparseCSC(1, 1) + + CommonOps_FSCC.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra.one, + other.toEjml().origin, + out, + null, + null, + ) + + return out.wrapVector() + } + + override fun Point.minus(other: Point): EjmlFloatVector { + val out = FMatrixSparseCSC(1, 1) + + CommonOps_FSCC.add( + elementAlgebra.one, + toEjml().origin, + elementAlgebra { -one }, + other.toEjml().origin, + out, + null, + null, + ) + + return out.wrapVector() + } + + override fun Float.times(m: Matrix): EjmlFloatMatrix = m * this + + override fun Point.times(value: Float): EjmlFloatVector { + val res = FMatrixSparseCSC(1, 1) + CommonOps_FSCC.scale(value, toEjml().origin, res) + return res.wrapVector() + } + + override fun Float.times(v: Point): EjmlFloatVector = v * this + + @UnstableKMathAPI + override fun computeFeature(structure: Matrix, type: KClass): F? { + structure.getFeature(type)?.let { return it } + val origin = structure.toEjml().origin + + return when (type) { + QRDecompositionFeature::class -> object : QRDecompositionFeature { + private val qr by lazy { + DecompositionFactory_FSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) } + } + + override val q: Matrix by lazy { + qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature) + } + + override val r: Matrix by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) } + } + + CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature { + override val l: Matrix by lazy { + val cholesky = + DecompositionFactory_FSCC.cholesky().apply { decompose(origin.copy()) } + + (cholesky.getT(null) as FMatrix).wrapMatrix().withFeature(LFeature) + } + } + + LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object : + LUDecompositionFeature, DeterminantFeature, InverseMatrixFeature { + private val lu by lazy { + DecompositionFactory_FSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) } + } + + override val l: Matrix by lazy { + lu.getLower(null).wrapMatrix().withFeature(LFeature) + } + + override val u: Matrix by lazy { + lu.getUpper(null).wrapMatrix().withFeature(UFeature) + } + + override val inverse: Matrix by lazy { + var a = origin + val inverse = FMatrixRMaj(1, 1) + val solver = LinearSolverFactory_FSCC.lu(FillReducing.NONE) + if (solver.modifiesA()) a = a.copy() + val i = CommonOps_FDRM.identity(a.numRows) + solver.solve(i, inverse) + inverse.wrapMatrix() + } + + override val determinant: Float by lazy { elementAlgebra.number(lu.computeDeterminant().real) } + } + + else -> null + }?.let(type::cast) + } + + /** + * Solves for *x* in the following equation: *x = [a] -1 · [b]*. + * + * @param a the base matrix. + * @param b n by p matrix. + * @return the solution for *x* that is n by p. + */ + public fun solve(a: Matrix, b: Matrix): EjmlFloatMatrix { + val res = FMatrixSparseCSC(1, 1) + CommonOps_FSCC.solve(FMatrixSparseCSC(a.toEjml().origin), FMatrixSparseCSC(b.toEjml().origin), res) + return res.wrapMatrix() + } + + /** + * Solves for *x* in the following equation: *x = [a] -1 · [b]*. + * + * @param a the base matrix. + * @param b n by p vector. + * @return the solution for *x* that is n by p. + */ + public fun solve(a: Matrix, b: Point): EjmlFloatVector { + val res = FMatrixSparseCSC(1, 1) + CommonOps_FSCC.solve(FMatrixSparseCSC(a.toEjml().origin), FMatrixSparseCSC(b.toEjml().origin), res) + return EjmlFloatVector(res) + } +} + diff --git a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt index 88932d59b..c1ee8b48f 100644 --- a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt +++ b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt @@ -13,9 +13,9 @@ import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.algebra +import space.kscience.kmath.operations.asIterable import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.DoubleBuffer -import space.kscience.kmath.structures.asIterable import kotlin.math.pow /* diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegratorRuleFactory.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegratorRuleFactory.kt index b09129626..94c73832b 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegratorRuleFactory.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegratorRuleFactory.kt @@ -5,10 +5,10 @@ package space.kscience.kmath.integration +import space.kscience.kmath.operations.map import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.DoubleBuffer import space.kscience.kmath.structures.asBuffer -import space.kscience.kmath.structures.map import kotlin.jvm.Synchronized import kotlin.math.ulp import kotlin.native.concurrent.ThreadLocal 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 662cdf3d7..6abe89aad 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 @@ -12,14 +12,10 @@ import space.kscience.kmath.interpolation.SplineInterpolator import space.kscience.kmath.interpolation.interpolatePolynomials import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.misc.UnstableKMathAPI -import space.kscience.kmath.operations.DoubleField -import space.kscience.kmath.operations.Field -import space.kscience.kmath.operations.invoke -import space.kscience.kmath.operations.sum +import space.kscience.kmath.operations.* import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.DoubleBuffer import space.kscience.kmath.structures.MutableBufferFactory -import space.kscience.kmath.structures.map /** * Compute analytical indefinite integral of this [PiecewisePolynomial], keeping all intervals intact diff --git a/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector2DTest.kt b/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector2DTest.kt index 84b1f4fd6..89ee23354 100644 --- a/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector2DTest.kt +++ b/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector2DTest.kt @@ -1,9 +1,8 @@ package space.kscience.kmath.geometry -import space.kscience.kmath.structures.asSequence -import space.kscience.kmath.structures.toList -import kotlin.test.assertEquals +import space.kscience.kmath.operations.toList import kotlin.test.Test +import kotlin.test.assertEquals internal class Vector2DTest { private val vector = Vector2D(1.0, -7.999) diff --git a/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector3DTest.kt b/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector3DTest.kt index 717e78871..70f8f4ebd 100644 --- a/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector3DTest.kt +++ b/kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector3DTest.kt @@ -1,6 +1,6 @@ package space.kscience.kmath.geometry -import space.kscience.kmath.structures.toList +import space.kscience.kmath.operations.toList import kotlin.test.Test import kotlin.test.assertEquals diff --git a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt b/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt index 0841fcb4c..d5b74fb9b 100644 --- a/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt +++ b/kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt @@ -7,8 +7,8 @@ package space.kscience.kmath.histogram import space.kscience.kmath.domains.UnivariateDomain import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.operations.asSequence import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.asSequence @UnstableKMathAPI @@ -34,7 +34,7 @@ public class UnivariateBin( } @OptIn(UnstableKMathAPI::class) -public interface UnivariateHistogram : Histogram{ +public interface UnivariateHistogram : Histogram { public operator fun get(value: Double): UnivariateBin? override operator fun get(point: Buffer): UnivariateBin? = get(point[0]) diff --git a/kmath-jupyter/src/main/kotlin/space/kscience/kmath/jupyter/KMathJupyter.kt b/kmath-jupyter/src/main/kotlin/space/kscience/kmath/jupyter/KMathJupyter.kt index e646d2bd0..9731908b3 100644 --- a/kmath-jupyter/src/main/kotlin/space/kscience/kmath/jupyter/KMathJupyter.kt +++ b/kmath-jupyter/src/main/kotlin/space/kscience/kmath/jupyter/KMathJupyter.kt @@ -21,9 +21,9 @@ import space.kscience.kmath.expressions.MST import space.kscience.kmath.expressions.MstRing import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.nd.Structure2D +import space.kscience.kmath.operations.asSequence import space.kscience.kmath.operations.invoke import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.asSequence /** * A function for conversion of number to MST for pretty print diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Median.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Median.kt index 54b2e42b3..664e4e8e7 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Median.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Median.kt @@ -5,8 +5,8 @@ package space.kscience.kmath.stat +import space.kscience.kmath.operations.asSequence import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.asSequence /** * Non-composable median 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 a0f5d8080..d31e02677 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 @@ -9,9 +9,9 @@ import space.kscience.kmath.nd.MutableStructure1D import space.kscience.kmath.nd.MutableStructure2D import space.kscience.kmath.nd.as1D import space.kscience.kmath.nd.as2D +import space.kscience.kmath.operations.asSequence import space.kscience.kmath.operations.invoke import space.kscience.kmath.structures.VirtualBuffer -import space.kscience.kmath.structures.asSequence import space.kscience.kmath.tensors.core.BufferedTensor import space.kscience.kmath.tensors.core.DoubleTensor import space.kscience.kmath.tensors.core.DoubleTensorAlgebra diff --git a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/utils.kt b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/utils.kt index 6088c32e4..8428dae5c 100644 --- a/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/utils.kt +++ b/kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/utils.kt @@ -6,6 +6,7 @@ package space.kscience.kmath.tensors.core.internal import space.kscience.kmath.nd.as1D +import space.kscience.kmath.operations.toMutableList import space.kscience.kmath.samplers.GaussianSampler import space.kscience.kmath.stat.RandomGenerator import space.kscience.kmath.structures.*