From 625e624cabacae0281313eff1c9f08baeeb5eaff Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sun, 29 Nov 2020 21:38:12 +0300 Subject: [PATCH 1/4] kmath-for-real refactoring --- CHANGELOG.md | 2 + .../structures/MultiplicationBenchmark.kt | 51 +++++++++++++++++++ .../kmath/commons/fit/fitWithAutoDiff.kt | 1 + .../kmath/linear/LinearAlgebraBenchmark.kt | 4 +- .../kmath/linear/MultiplicationBenchmark.kt | 38 -------------- .../kscience/kmath/linear/BufferMatrix.kt | 37 +++----------- .../kscience/kmath/linear/FeaturedMatrix.kt | 7 +-- .../kscience/kmath/linear/MatrixTest.kt | 3 +- .../kscience/kmath/ejml/EjmlMatrixContext.kt | 20 ++++---- .../kotlin/kscience/kmath/real/RealMatrix.kt | 44 ++++++++++------ .../kotlin/kscience/kmath/real/RealVector.kt | 37 +++++++++----- .../kotlin/kscience/kmath/real/dot.kt | 31 +++++++++++ .../kaceince/kmath/real/RealVectorTest.kt | 15 +++--- .../kscience/kmath/histogram/RealHistogram.kt | 5 +- .../histogram/MultivariateHistogramTest.kt | 2 + .../kmath/histogram/UnivariateHistogram.kt | 4 +- 16 files changed, 172 insertions(+), 129 deletions(-) create mode 100644 examples/src/benchmarks/kotlin/kscience/kmath/structures/MultiplicationBenchmark.kt delete mode 100644 examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt create mode 100644 kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/dot.kt diff --git a/CHANGELOG.md b/CHANGELOG.md index 99511a161..1fd169b10 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,8 @@ - `kmath-prob` renamed to `kmath-stat` - Grid generators moved to `kmath-for-real` - Use `Point` instead of specialized type in `kmath-for-real` +- Optimized dot product for buffer matrices moved to `kmath-for-real` +- EjmlMatrix context is an object ### Deprecated diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/structures/MultiplicationBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/structures/MultiplicationBenchmark.kt new file mode 100644 index 000000000..3c2c12c0b --- /dev/null +++ b/examples/src/benchmarks/kotlin/kscience/kmath/structures/MultiplicationBenchmark.kt @@ -0,0 +1,51 @@ +package kscience.kmath.structures + +import kotlinx.benchmark.Benchmark +import kscience.kmath.commons.linear.CMMatrixContext +import kscience.kmath.commons.linear.CMMatrixContext.dot +import kscience.kmath.commons.linear.toCM +import kscience.kmath.ejml.EjmlMatrixContext +import kscience.kmath.ejml.toEjml +import kscience.kmath.linear.real +import kscience.kmath.operations.RealField +import kscience.kmath.operations.invoke +import org.openjdk.jmh.annotations.Scope +import org.openjdk.jmh.annotations.State +import kotlin.random.Random + +@State(Scope.Benchmark) +class MultiplicationBenchmark { + companion object { + val random = Random(12224) + val dim = 1000 + + //creating invertible matrix + val matrix1 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } + val matrix2 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } + + val cmMatrix1 = matrix1.toCM() + val cmMatrix2 = matrix2.toCM() + + val ejmlMatrix1 = matrix1.toEjml() + val ejmlMatrix2 = matrix2.toEjml() + } + + @Benchmark + fun commonsMathMultiplication() { + CMMatrixContext.invoke { + cmMatrix1 dot cmMatrix2 + } + } + + @Benchmark + fun ejmlMultiplication() { + EjmlMatrixContext.invoke { + ejmlMatrix1 dot ejmlMatrix2 + } + } + + @Benchmark + fun bufferedMultiplication() { + matrix1 dot matrix2 + } +} \ No newline at end of file diff --git a/examples/src/main/kotlin/kscience/kmath/commons/fit/fitWithAutoDiff.kt b/examples/src/main/kotlin/kscience/kmath/commons/fit/fitWithAutoDiff.kt index 6cfa1961b..c0cd9dc5c 100644 --- a/examples/src/main/kotlin/kscience/kmath/commons/fit/fitWithAutoDiff.kt +++ b/examples/src/main/kotlin/kscience/kmath/commons/fit/fitWithAutoDiff.kt @@ -70,6 +70,7 @@ fun main() { //minimize the chi^2 in given starting point. Derivatives are not required, they are already included. val result: OptimizationResult = chi2.minimize(a to 1.5, b to 0.9, c to 1.0) + //display a page with plot and numerical results val page = Plotly.page { plot { scatter { diff --git a/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt b/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt index 3316f3236..5b1fc9bc7 100644 --- a/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt +++ b/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt @@ -1,10 +1,12 @@ package kscience.kmath.linear import kscience.kmath.commons.linear.CMMatrixContext +import kscience.kmath.commons.linear.CMMatrixContext.dot import kscience.kmath.commons.linear.inverse import kscience.kmath.commons.linear.toCM import kscience.kmath.ejml.EjmlMatrixContext import kscience.kmath.ejml.inverse +import kscience.kmath.ejml.toEjml import kscience.kmath.operations.RealField import kscience.kmath.operations.invoke import kscience.kmath.structures.Matrix @@ -40,7 +42,7 @@ fun main() { println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis") val ejmlTime = measureTimeMillis { - (EjmlMatrixContext(RealField)) { + EjmlMatrixContext { val km = matrix.toEjml() //avoid overhead on conversion repeat(n) { inverse(km) } } diff --git a/examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt b/examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt deleted file mode 100644 index d1011e8f5..000000000 --- a/examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt +++ /dev/null @@ -1,38 +0,0 @@ -package kscience.kmath.linear - -import kscience.kmath.commons.linear.CMMatrixContext -import kscience.kmath.commons.linear.toCM -import kscience.kmath.ejml.EjmlMatrixContext -import kscience.kmath.operations.RealField -import kscience.kmath.operations.invoke -import kscience.kmath.structures.Matrix -import kotlin.random.Random -import kotlin.system.measureTimeMillis - -fun main() { - val random = Random(12224) - val dim = 1000 - //creating invertible matrix - val matrix1 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } - val matrix2 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } - -// //warmup -// matrix1 dot matrix2 - - CMMatrixContext { - val cmMatrix1 = matrix1.toCM() - val cmMatrix2 = matrix2.toCM() - val cmTime = measureTimeMillis { cmMatrix1 dot cmMatrix2 } - println("CM implementation time: $cmTime") - } - - (EjmlMatrixContext(RealField)) { - val ejmlMatrix1 = matrix1.toEjml() - val ejmlMatrix2 = matrix2.toEjml() - val ejmlTime = measureTimeMillis { ejmlMatrix1 dot ejmlMatrix2 } - println("EJML implementation time: $ejmlTime") - } - - val genericTime = measureTimeMillis { val res = matrix1 dot matrix2 } - println("Generic implementation time: $genericTime") -} diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt index d51f40890..f6ed34be0 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt @@ -9,7 +9,7 @@ import kscience.kmath.structures.* */ public class BufferMatrixContext>( public override val elementContext: R, - private val bufferFactory: BufferFactory + private val bufferFactory: BufferFactory, ) : GenericMatrixContext { public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): BufferMatrix { val buffer = bufferFactory(rows * columns) { offset -> initializer(offset / columns, offset % columns) } @@ -29,8 +29,8 @@ public object RealMatrixContext : GenericMatrixContext { public override inline fun produce( rows: Int, columns: Int, - initializer: (i: Int, j: Int) -> Double - ): Matrix { + initializer: (i: Int, j: Int) -> Double, + ): BufferMatrix { val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) } return BufferMatrix(rows, columns, buffer) } @@ -43,15 +43,15 @@ public class BufferMatrix( public override val rowNum: Int, public override val colNum: Int, public val buffer: Buffer, - public override val features: Set = emptySet() + public override val features: Set = emptySet(), ) : FeaturedMatrix { - override val shape: IntArray - get() = intArrayOf(rowNum, colNum) init { require(buffer.size == rowNum * colNum) { "Dimension mismatch for matrix structure" } } + override val shape: IntArray get() = intArrayOf(rowNum, colNum) + public override fun suggestFeature(vararg features: MatrixFeature): BufferMatrix = BufferMatrix(rowNum, colNum, buffer, this.features + features) @@ -86,28 +86,3 @@ public class BufferMatrix( else "Matrix(rowsNum = $rowNum, colNum = $colNum, features=$features)" } } - -/** - * Optimized dot product for real matrices - */ -public infix fun BufferMatrix.dot(other: BufferMatrix): BufferMatrix { - require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" } - val array = DoubleArray(this.rowNum * other.colNum) - - //convert to array to insure there is not memory indirection - fun Buffer.unsafeArray() = if (this is RealBuffer) - array - else - DoubleArray(size) { get(it) } - - val a = this.buffer.unsafeArray() - val b = other.buffer.unsafeArray() - - for (i in (0 until rowNum)) - for (j in (0 until other.colNum)) - for (k in (0 until colNum)) - array[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j] - - val buffer = RealBuffer(array) - return BufferMatrix(rowNum, other.colNum, buffer) -} diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt index 5d9af8608..003741979 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt @@ -27,9 +27,8 @@ public interface FeaturedMatrix : Matrix { public inline fun Structure2D.Companion.real( rows: Int, columns: Int, - initializer: (Int, Int) -> Double -): Matrix = - MatrixContext.real.produce(rows, columns, initializer) + initializer: (Int, Int) -> Double, +): BufferMatrix = MatrixContext.real.produce(rows, columns, initializer) /** * Build a square matrix from given elements. @@ -82,5 +81,3 @@ public fun Matrix.transpose(): Matrix { setOf(TransposedFeature(this)) ) { i, j -> get(j, i) } } - -public infix fun Matrix.dot(other: Matrix): Matrix = with(MatrixContext.real) { dot(other) } diff --git a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt index 7cfa25a66..0a582e339 100644 --- a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt +++ b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt @@ -1,5 +1,6 @@ package kscience.kmath.linear +import kscience.kmath.operations.invoke import kscience.kmath.structures.Matrix import kscience.kmath.structures.NDStructure import kscience.kmath.structures.as2D @@ -38,7 +39,7 @@ class MatrixTest { infix fun Matrix.pow(power: Int): Matrix { var res = this repeat(power - 1) { - res = res dot this + res = RealMatrixContext.invoke { res dot this@pow } } return res } diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt index 52826a7b1..1ed57fcb3 100644 --- a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt +++ b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt @@ -7,17 +7,18 @@ import kscience.kmath.operations.Space import kscience.kmath.operations.invoke import kscience.kmath.structures.Matrix +/** + * Converts this matrix to EJML one. + */ +public fun Matrix.toEjml(): EjmlMatrix = + if (this is EjmlMatrix) this else EjmlMatrixContext.produce(rowNum, colNum) { i, j -> get(i, j) } + /** * Represents context of basic operations operating with [EjmlMatrix]. * * @author Iaroslav Postovalov */ -public class EjmlMatrixContext(private val space: Space) : MatrixContext { - /** - * Converts this matrix to EJML one. - */ - public fun Matrix.toEjml(): EjmlMatrix = - if (this is EjmlMatrix) this else produce(rowNum, colNum) { i, j -> get(i, j) } +public object EjmlMatrixContext : MatrixContext { /** * Converts this vector to EJML one. @@ -47,11 +48,10 @@ public class EjmlMatrixContext(private val space: Space) : MatrixContext EjmlMatrix(toEjml().origin - b.toEjml().origin) public override fun multiply(a: Matrix, k: Number): EjmlMatrix = - produce(a.rowNum, a.colNum) { i, j -> space { a[i, j] * k } } + produce(a.rowNum, a.colNum) { i, j -> a[i, j] * k.toDouble() } - public override operator fun Matrix.times(value: Double): EjmlMatrix = EjmlMatrix(toEjml().origin.scale(value)) - - public companion object + public override operator fun Matrix.times(value: Double): EjmlMatrix = + EjmlMatrix(toEjml().origin.scale(value)) } /** diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt index 784b27f9d..c0657c0bf 100644 --- a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt +++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt @@ -6,10 +6,7 @@ import kscience.kmath.linear.VirtualMatrix import kscience.kmath.misc.UnstableKMathAPI import kscience.kmath.operations.invoke import kscience.kmath.operations.sum -import kscience.kmath.structures.Buffer -import kscience.kmath.structures.Matrix -import kscience.kmath.structures.RealBuffer -import kscience.kmath.structures.asIterable +import kscience.kmath.structures.* import kotlin.math.pow /* @@ -86,18 +83,6 @@ public operator fun Double.minus(matrix: RealMatrix): RealMatrix = // row, col -> matrix[row, col] / this //} -/* - * Per-element (!) square and power operations - */ - -public fun RealMatrix.square(): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { row, col -> - this[row, col].pow(2) -} - -public fun RealMatrix.pow(n: Int): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { i, j -> - this[i, j].pow(n) -} - /* * Operations on two matrices (per-element!) */ @@ -157,3 +142,30 @@ public fun RealMatrix.sum(): Double = elements().map { (_, value) -> value }.sum public fun RealMatrix.min(): Double? = elements().map { (_, value) -> value }.minOrNull() public fun RealMatrix.max(): Double? = elements().map { (_, value) -> value }.maxOrNull() public fun RealMatrix.average(): Double = elements().map { (_, value) -> value }.average() + +public inline fun RealMatrix.map(transform: (Double) -> Double): RealMatrix = + MatrixContext.real.produce(rowNum, colNum) { i, j -> + transform(get(i, j)) + } + +//extended operations + +public fun RealMatrix.pow(p: Double): RealMatrix = map { it.pow(p) } + +public fun RealMatrix.pow(p: Int): RealMatrix = map { it.pow(p) } + +public fun exp(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.exp(it) } + +public fun sqrt(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.sqrt(it) } + +public fun RealMatrix.square(): RealMatrix = map { it.pow(2) } + +public fun sin(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.sin(it) } + +public fun cos(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.cos(it) } + +public fun tan(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.tan(it) } + +public fun ln(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.ln(it) } + +public fun log10(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.log10(it) } \ No newline at end of file diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt index 1182d871c..596692782 100644 --- a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt +++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt @@ -3,7 +3,6 @@ package kscience.kmath.real import kscience.kmath.linear.Point import kscience.kmath.operations.Norm import kscience.kmath.structures.Buffer -import kscience.kmath.structures.RealBuffer import kscience.kmath.structures.asBuffer import kscience.kmath.structures.asIterable import kotlin.math.pow @@ -11,17 +10,11 @@ import kotlin.math.sqrt public typealias RealVector = Point -public inline fun RealVector(size: Int, init: (Int) -> Double): RealVector = RealBuffer(size, init) -public fun RealVector(vararg doubles: Double): RealVector = RealBuffer(doubles) - -public fun DoubleArray.asVector(): RealVector = asBuffer() -public fun List.asVector(): RealVector = asBuffer() - - public object VectorL2Norm : Norm, Double> { override fun norm(arg: Point): Double = sqrt(arg.asIterable().sumByDouble(Number::toDouble)) } +public operator fun Buffer.Companion.invoke(vararg doubles: Double): RealVector = doubles.asBuffer() /** * Fill the vector of given [size] with given [value] @@ -36,12 +29,6 @@ public inline fun RealVector.map(transform: (Double) -> Double): RealVector = public inline fun RealVector.mapIndexed(transform: (index: Int, value: Double) -> Double): RealVector = Buffer.real(size) { transform(it, get(it)) } -public fun RealVector.pow(p: Double): RealVector = map { it.pow(p) } - -public fun RealVector.pow(p: Int): RealVector = map { it.pow(p) } - -public fun exp(vector: RealVector): RealVector = vector.map { kotlin.math.exp(it) } - public operator fun RealVector.plus(other: RealVector): RealVector = mapIndexed { index, value -> value + other[index] } @@ -71,3 +58,25 @@ public operator fun RealVector.div(other: RealVector): RealVector = public operator fun RealVector.div(number: Number): RealVector = map { it / number.toDouble() } public operator fun Number.div(vector: RealVector): RealVector = vector.map { toDouble() / it } + +//extended operations + +public fun RealVector.pow(p: Double): RealVector = map { it.pow(p) } + +public fun RealVector.pow(p: Int): RealVector = map { it.pow(p) } + +public fun exp(vector: RealVector): RealVector = vector.map { kotlin.math.exp(it) } + +public fun sqrt(vector: RealVector): RealVector = vector.map { kotlin.math.sqrt(it) } + +public fun RealVector.square(): RealVector = map { it.pow(2) } + +public fun sin(vector: RealVector): RealVector = vector.map { kotlin.math.sin(it) } + +public fun cos(vector: RealVector): RealVector = vector.map { kotlin.math.cos(it) } + +public fun tan(vector: RealVector): RealVector = vector.map { kotlin.math.tan(it) } + +public fun ln(vector: RealVector): RealVector = vector.map { kotlin.math.ln(it) } + +public fun log10(vector: RealVector): RealVector = vector.map { kotlin.math.log10(it) } diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/dot.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/dot.kt new file mode 100644 index 000000000..9beffe6bb --- /dev/null +++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/dot.kt @@ -0,0 +1,31 @@ +package kscience.kmath.real + +import kscience.kmath.linear.BufferMatrix +import kscience.kmath.structures.Buffer +import kscience.kmath.structures.RealBuffer + + +/** + * Optimized dot product for real matrices + */ +public infix fun BufferMatrix.dot(other: BufferMatrix): BufferMatrix { + require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" } + val resultArray = DoubleArray(this.rowNum * other.colNum) + + //convert to array to insure there is no memory indirection + fun Buffer.unsafeArray() = if (this is RealBuffer) + this.array + else + DoubleArray(size) { get(it) } + + val a = this.buffer.unsafeArray() + val b = other.buffer.unsafeArray() + + for (i in (0 until rowNum)) + for (j in (0 until other.colNum)) + for (k in (0 until colNum)) + resultArray[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j] + + val buffer = RealBuffer(resultArray) + return BufferMatrix(rowNum, other.colNum, buffer) +} \ No newline at end of file diff --git a/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt b/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt index 002db935c..6215ba5e8 100644 --- a/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt +++ b/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt @@ -1,34 +1,33 @@ package kaceince.kmath.real -import kscience.kmath.linear.MatrixContext -import kscience.kmath.linear.asMatrix -import kscience.kmath.linear.transpose +import kscience.kmath.linear.* import kscience.kmath.operations.invoke import kscience.kmath.real.RealVector import kscience.kmath.real.plus +import kscience.kmath.structures.Buffer import kotlin.test.Test import kotlin.test.assertEquals internal class RealVectorTest { @Test fun testSum() { - val vector1 = RealVector(5) { it.toDouble() } - val vector2 = RealVector(5) { 5 - it.toDouble() } + val vector1 = Buffer.real(5) { it.toDouble() } + val vector2 = Buffer.real(5) { 5 - it.toDouble() } val sum = vector1 + vector2 assertEquals(5.0, sum[2]) } @Test fun testVectorToMatrix() { - val vector = RealVector(5) { it.toDouble() } + val vector = Buffer.real(5) { it.toDouble() } val matrix = vector.asMatrix() assertEquals(4.0, matrix[4, 0]) } @Test fun testDot() { - val vector1 = RealVector(5) { it.toDouble() } - val vector2 = RealVector(5) { 5 - it.toDouble() } + val vector1 = Buffer.real(5) { it.toDouble() } + val vector2 = Buffer.real(5) { 5 - it.toDouble() } val matrix1 = vector1.asMatrix() val matrix2 = vector2.asMatrix().transpose() val product = MatrixContext.real { matrix1 dot matrix2 } diff --git a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt index c58952ab4..f95264ee1 100644 --- a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt +++ b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt @@ -3,7 +3,6 @@ package kscience.kmath.histogram import kscience.kmath.linear.Point import kscience.kmath.operations.SpaceOperations import kscience.kmath.operations.invoke -import kscience.kmath.real.asVector import kscience.kmath.structures.* import kotlin.math.floor @@ -123,8 +122,8 @@ public class RealHistogram( *``` */ public fun fromRanges(vararg ranges: ClosedFloatingPointRange): RealHistogram = RealHistogram( - ranges.map(ClosedFloatingPointRange::start).asVector(), - ranges.map(ClosedFloatingPointRange::endInclusive).asVector() + ranges.map(ClosedFloatingPointRange::start).asBuffer(), + ranges.map(ClosedFloatingPointRange::endInclusive).asBuffer() ) /** diff --git a/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt b/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt index eebb41019..af22afc6b 100644 --- a/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt +++ b/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt @@ -4,6 +4,8 @@ import kscience.kmath.histogram.RealHistogram import kscience.kmath.histogram.fill import kscience.kmath.histogram.put import kscience.kmath.real.RealVector +import kscience.kmath.real.invoke +import kscience.kmath.structures.Buffer import kotlin.random.Random import kotlin.test.* diff --git a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt b/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt index 5fada1302..2f3855892 100644 --- a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt +++ b/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt @@ -1,8 +1,8 @@ package kscience.kmath.histogram import kscience.kmath.real.RealVector -import kscience.kmath.real.asVector import kscience.kmath.structures.Buffer +import kscience.kmath.structures.asBuffer import java.util.* import kotlin.math.floor @@ -16,7 +16,7 @@ public class UnivariateBin( //TODO add weighting public override val value: Number get() = counter.sum() - public override val center: RealVector get() = doubleArrayOf(position).asVector() + public override val center: RealVector get() = doubleArrayOf(position).asBuffer() public override val dimension: Int get() = 1 public operator fun contains(value: Double): Boolean = value in (position - size / 2)..(position + size / 2) From 89c0d863d2ddc7ae71f60c7b6d4b28cb30b7a660 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sun, 29 Nov 2020 22:26:05 +0300 Subject: [PATCH 2/4] Benchmark refactoring --- examples/build.gradle.kts | 2 +- .../ArrayBenchmark.kt | 2 +- .../BufferBenchmark.kt | 4 +- .../benchmarks/LinearAlgebraBenchmark.kt | 50 ++++++++++++++++++ .../MultiplicationBenchmark.kt | 13 ++++- .../NDFieldBenchmark.kt | 3 +- .../ViktorBenchmark.kt | 7 ++- .../kotlin/kscience/kmath/utils/utils.kt | 11 ---- .../kmath/linear/LinearAlgebraBenchmark.kt | 52 ------------------- .../kscience/kmath/linear/FeaturedMatrix.kt | 2 +- 10 files changed, 74 insertions(+), 72 deletions(-) rename examples/src/benchmarks/kotlin/kscience/kmath/{structures => benchmarks}/ArrayBenchmark.kt (96%) rename examples/src/benchmarks/kotlin/kscience/kmath/{structures => benchmarks}/BufferBenchmark.kt (85%) create mode 100644 examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt rename examples/src/benchmarks/kotlin/kscience/kmath/{structures => benchmarks}/MultiplicationBenchmark.kt (81%) rename examples/src/benchmarks/kotlin/kscience/kmath/{structures => benchmarks}/NDFieldBenchmark.kt (94%) rename examples/src/benchmarks/kotlin/kscience/kmath/{structures => benchmarks}/ViktorBenchmark.kt (86%) delete mode 100644 examples/src/benchmarks/kotlin/kscience/kmath/utils/utils.kt delete mode 100644 examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts index eead68a9d..c079eaa84 100644 --- a/examples/build.gradle.kts +++ b/examples/build.gradle.kts @@ -69,7 +69,7 @@ benchmark { // This one matches sourceSet name above configurations.register("fast") { - warmups = 5 // number of warmup iterations + warmups = 1 // number of warmup iterations iterations = 3 // number of iterations iterationTime = 500 // time in seconds per iteration iterationTimeUnit = "ms" // time unity for iterationTime, default is seconds diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/structures/ArrayBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ArrayBenchmark.kt similarity index 96% rename from examples/src/benchmarks/kotlin/kscience/kmath/structures/ArrayBenchmark.kt rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ArrayBenchmark.kt index a91d02253..8c44135fb 100644 --- a/examples/src/benchmarks/kotlin/kscience/kmath/structures/ArrayBenchmark.kt +++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ArrayBenchmark.kt @@ -1,4 +1,4 @@ -package kscience.kmath.structures +package kscience.kmath.benchmarks import org.openjdk.jmh.annotations.Benchmark import org.openjdk.jmh.annotations.Scope diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/structures/BufferBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/BufferBenchmark.kt similarity index 85% rename from examples/src/benchmarks/kotlin/kscience/kmath/structures/BufferBenchmark.kt rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/BufferBenchmark.kt index 8b6fd4a51..4c64517f1 100644 --- a/examples/src/benchmarks/kotlin/kscience/kmath/structures/BufferBenchmark.kt +++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/BufferBenchmark.kt @@ -1,7 +1,9 @@ -package kscience.kmath.structures +package kscience.kmath.benchmarks import kscience.kmath.operations.Complex import kscience.kmath.operations.complex +import kscience.kmath.structures.MutableBuffer +import kscience.kmath.structures.RealBuffer import org.openjdk.jmh.annotations.Benchmark import org.openjdk.jmh.annotations.Scope import org.openjdk.jmh.annotations.State diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt new file mode 100644 index 000000000..0ac6ea5f0 --- /dev/null +++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt @@ -0,0 +1,50 @@ +package kscience.kmath.linear + + +import kotlinx.benchmark.Benchmark +import kscience.kmath.commons.linear.CMMatrixContext +import kscience.kmath.commons.linear.CMMatrixContext.dot +import kscience.kmath.commons.linear.inverse +import kscience.kmath.commons.linear.toCM +import kscience.kmath.ejml.EjmlMatrixContext +import kscience.kmath.ejml.inverse +import kscience.kmath.ejml.toEjml +import kscience.kmath.operations.invoke +import kscience.kmath.structures.Matrix +import org.openjdk.jmh.annotations.Scope +import org.openjdk.jmh.annotations.State +import kotlin.random.Random + +@State(Scope.Benchmark) +class LinearAlgebraBenchmark { + companion object { + val random = Random(1224) + val dim = 100 + + //creating invertible matrix + val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } + val l = Matrix.real(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 } + val matrix = l dot u + } + + @Benchmark + fun kmathLUPInversion() { + MatrixContext.real.inverse(matrix) + } + + @Benchmark + fun cmLUPInversion() { + CMMatrixContext { + val cm = matrix.toCM() //avoid overhead on conversion + inverse(cm) + } + } + + @Benchmark + fun ejmlInverse() { + EjmlMatrixContext { + val km = matrix.toEjml() //avoid overhead on conversion + inverse(km) + } + } +} \ No newline at end of file diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/structures/MultiplicationBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt similarity index 81% rename from examples/src/benchmarks/kotlin/kscience/kmath/structures/MultiplicationBenchmark.kt rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt index 3c2c12c0b..9d2b02245 100644 --- a/examples/src/benchmarks/kotlin/kscience/kmath/structures/MultiplicationBenchmark.kt +++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt @@ -1,4 +1,4 @@ -package kscience.kmath.structures +package kscience.kmath.benchmarks import kotlinx.benchmark.Benchmark import kscience.kmath.commons.linear.CMMatrixContext @@ -7,8 +7,8 @@ import kscience.kmath.commons.linear.toCM import kscience.kmath.ejml.EjmlMatrixContext import kscience.kmath.ejml.toEjml import kscience.kmath.linear.real -import kscience.kmath.operations.RealField import kscience.kmath.operations.invoke +import kscience.kmath.structures.Matrix import org.openjdk.jmh.annotations.Scope import org.openjdk.jmh.annotations.State import kotlin.random.Random @@ -44,6 +44,15 @@ class MultiplicationBenchmark { } } + @Benchmark + fun ejmlMultiplicationwithConversion() { + val ejmlMatrix1 = matrix1.toEjml() + val ejmlMatrix2 = matrix2.toEjml() + EjmlMatrixContext.invoke { + ejmlMatrix1 dot ejmlMatrix2 + } + } + @Benchmark fun bufferedMultiplication() { matrix1 dot matrix2 diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/structures/NDFieldBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt similarity index 94% rename from examples/src/benchmarks/kotlin/kscience/kmath/structures/NDFieldBenchmark.kt rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt index 8ec47ae81..1be8e7236 100644 --- a/examples/src/benchmarks/kotlin/kscience/kmath/structures/NDFieldBenchmark.kt +++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt @@ -1,7 +1,8 @@ -package kscience.kmath.structures +package kscience.kmath.benchmarks import kscience.kmath.operations.RealField import kscience.kmath.operations.invoke +import kscience.kmath.structures.* import org.openjdk.jmh.annotations.Benchmark import org.openjdk.jmh.annotations.Scope import org.openjdk.jmh.annotations.State diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/structures/ViktorBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt similarity index 86% rename from examples/src/benchmarks/kotlin/kscience/kmath/structures/ViktorBenchmark.kt rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt index 464925ca0..8663e353c 100644 --- a/examples/src/benchmarks/kotlin/kscience/kmath/structures/ViktorBenchmark.kt +++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt @@ -1,7 +1,10 @@ -package kscience.kmath.structures +package kscience.kmath.benchmarks import kscience.kmath.operations.RealField import kscience.kmath.operations.invoke +import kscience.kmath.structures.BufferedNDField +import kscience.kmath.structures.NDField +import kscience.kmath.structures.RealNDField import kscience.kmath.viktor.ViktorNDField import org.jetbrains.bio.viktor.F64Array import org.openjdk.jmh.annotations.Benchmark @@ -36,7 +39,7 @@ internal class ViktorBenchmark { @Benchmark fun rawViktor() { - val one = F64Array.full(init = 1.0, shape = *intArrayOf(dim, dim)) + val one = F64Array.full(init = 1.0, shape = intArrayOf(dim, dim)) var res = one repeat(n) { res = res + one } } diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/utils/utils.kt b/examples/src/benchmarks/kotlin/kscience/kmath/utils/utils.kt deleted file mode 100644 index 329b7b17b..000000000 --- a/examples/src/benchmarks/kotlin/kscience/kmath/utils/utils.kt +++ /dev/null @@ -1,11 +0,0 @@ -package kscience.kmath.utils - -import kotlin.contracts.InvocationKind -import kotlin.contracts.contract -import kotlin.system.measureTimeMillis - -internal inline fun measureAndPrint(title: String, block: () -> Unit) { - contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) } - val time = measureTimeMillis(block) - println("$title completed in $time millis") -} diff --git a/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt b/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt deleted file mode 100644 index 5b1fc9bc7..000000000 --- a/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt +++ /dev/null @@ -1,52 +0,0 @@ -package kscience.kmath.linear - -import kscience.kmath.commons.linear.CMMatrixContext -import kscience.kmath.commons.linear.CMMatrixContext.dot -import kscience.kmath.commons.linear.inverse -import kscience.kmath.commons.linear.toCM -import kscience.kmath.ejml.EjmlMatrixContext -import kscience.kmath.ejml.inverse -import kscience.kmath.ejml.toEjml -import kscience.kmath.operations.RealField -import kscience.kmath.operations.invoke -import kscience.kmath.structures.Matrix -import kotlin.random.Random -import kotlin.system.measureTimeMillis - -fun main() { - val random = Random(1224) - val dim = 100 - //creating invertible matrix - val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } - val l = Matrix.real(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 } - val matrix = l dot u - - val n = 5000 // iterations - - MatrixContext.real { - repeat(50) { inverse(matrix) } - val inverseTime = measureTimeMillis { repeat(n) { inverse(matrix) } } - println("[kmath] Inversion of $n matrices $dim x $dim finished in $inverseTime millis") - } - - //commons-math - - val commonsTime = measureTimeMillis { - CMMatrixContext { - val cm = matrix.toCM() //avoid overhead on conversion - repeat(n) { inverse(cm) } - } - } - - - println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis") - - val ejmlTime = measureTimeMillis { - EjmlMatrixContext { - val km = matrix.toEjml() //avoid overhead on conversion - repeat(n) { inverse(km) } - } - } - - println("[ejml] Inversion of $n matrices $dim x $dim finished in $ejmlTime millis") -} diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt index 003741979..350379717 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt @@ -80,4 +80,4 @@ public fun Matrix.transpose(): Matrix { rowNum, setOf(TransposedFeature(this)) ) { i, j -> get(j, i) } -} +} \ No newline at end of file From 5368bb5d4a19363d4448e9f985ae04a1e15a799a Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Mon, 30 Nov 2020 11:25:21 +0300 Subject: [PATCH 3/4] Readme update --- kmath-ast/build.gradle.kts | 4 ++++ kmath-ast/src/jvmMain/kotlin/kscience/kmath/ast/parser.kt | 4 +--- kmath-core/build.gradle.kts | 4 +--- kmath-dimensions/build.gradle.kts | 4 ++++ 4 files changed, 10 insertions(+), 6 deletions(-) diff --git a/kmath-ast/build.gradle.kts b/kmath-ast/build.gradle.kts index 2f95d427d..3e3c0475f 100644 --- a/kmath-ast/build.gradle.kts +++ b/kmath-ast/build.gradle.kts @@ -17,3 +17,7 @@ kotlin.sourceSets { } } } + +readme{ + maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE +} \ No newline at end of file diff --git a/kmath-ast/src/jvmMain/kotlin/kscience/kmath/ast/parser.kt b/kmath-ast/src/jvmMain/kotlin/kscience/kmath/ast/parser.kt index 94cd3b321..79cf77a6a 100644 --- a/kmath-ast/src/jvmMain/kotlin/kscience/kmath/ast/parser.kt +++ b/kmath-ast/src/jvmMain/kotlin/kscience/kmath/ast/parser.kt @@ -83,11 +83,10 @@ public object ArithmeticsEvaluator : Grammar() { } /** - * Tries to parse the string into [MST]. + * Tries to parse the string into [MST]. Returns [ParseResult] representing expression or error. * * @receiver the string to parse. * @return the [MST] node. - * @author Alexander Nozik */ public fun String.tryParseMath(): ParseResult = ArithmeticsEvaluator.tryParseToEnd(this) @@ -96,6 +95,5 @@ public fun String.tryParseMath(): ParseResult = ArithmeticsEvaluator.tryPar * * @receiver the string to parse. * @return the [MST] node. - * @author Alexander Nozik */ public fun String.parseMath(): MST = ArithmeticsEvaluator.parseToEnd(this) diff --git a/kmath-core/build.gradle.kts b/kmath-core/build.gradle.kts index 7f889d9b4..9ed7e690b 100644 --- a/kmath-core/build.gradle.kts +++ b/kmath-core/build.gradle.kts @@ -1,5 +1,3 @@ -import ru.mipt.npm.gradle.Maturity - plugins { id("ru.mipt.npm.mpp") id("ru.mipt.npm.native") @@ -13,7 +11,7 @@ kotlin.sourceSets.commonMain { readme { description = "Core classes, algebra definitions, basic linear algebra" - maturity = Maturity.DEVELOPMENT + maturity = ru.mipt.npm.gradle.Maturity.DEVELOPMENT propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md")) feature( diff --git a/kmath-dimensions/build.gradle.kts b/kmath-dimensions/build.gradle.kts index 412d7162f..9bf89fc43 100644 --- a/kmath-dimensions/build.gradle.kts +++ b/kmath-dimensions/build.gradle.kts @@ -18,3 +18,7 @@ kotlin.sourceSets { } } } + +readme{ + maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE +} From 712df041704cf6194610eb76009aa028c93dc331 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Tue, 1 Dec 2020 21:21:56 +0300 Subject: [PATCH 4/4] Implemented #131 --- CHANGELOG.md | 7 +- .../benchmarks/LinearAlgebraBenchmark.kt | 2 +- .../kscience/kmath/commons/linear/CMMatrix.kt | 7 +- .../kscience/kmath/linear/BufferMatrix.kt | 4 +- .../kscience/kmath/linear/FeaturedMatrix.kt | 4 +- .../kscience/kmath/linear/LUPDecomposition.kt | 86 ++++++++++--------- .../kscience/kmath/linear/MatrixContext.kt | 37 ++++---- .../kmath/operations/AlgebraExtensions.kt | 5 ++ .../kmath/structures/BufferAccessor2D.kt | 22 +++-- .../kscience/kmath/linear/RealLUSolverTest.kt | 4 +- .../kscience/kmath/dimensions/Wrappers.kt | 4 +- .../kscience/kmath/ejml/EjmlMatrixContext.kt | 6 +- .../kotlin/kscience/kmath/real/RealMatrix.kt | 13 ++- 13 files changed, 112 insertions(+), 89 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1fd169b10..e542d210c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,13 +16,14 @@ - ND4J support module submitting `NDStructure` and `NDAlgebra` over `INDArray`. - Coroutine-deterministic Monte-Carlo scope with a random number generator. - Some minor utilities to `kmath-for-real`. +- Generic operation result parameter to `MatrixContext` ### Changed - Package changed from `scientifik` to `kscience.kmath`. -- Gradle version: 6.6 -> 6.7 +- Gradle version: 6.6 -> 6.7.1 - Minor exceptions refactor (throwing `IllegalArgumentException` by argument checks instead of `IllegalStateException`) - `Polynomial` secondary constructor made function. -- Kotlin version: 1.3.72 -> 1.4.20-M1 +- Kotlin version: 1.3.72 -> 1.4.20 - `kmath-ast` doesn't depend on heavy `kotlin-reflect` library. - Full autodiff refactoring based on `Symbol` - `kmath-prob` renamed to `kmath-stat` @@ -30,6 +31,7 @@ - Use `Point` instead of specialized type in `kmath-for-real` - Optimized dot product for buffer matrices moved to `kmath-for-real` - EjmlMatrix context is an object +- Matrix LUP `inverse` renamed to `inverseWithLUP` ### Deprecated @@ -37,6 +39,7 @@ - `kmath-koma` module because it doesn't support Kotlin 1.4. - Support of `legacy` JS backend (we will support only IR) - `toGrid` method. +- Public visibility of `BufferAccessor2D` ### Fixed - `symbol` method in `MstExtendedField` (https://github.com/mipt-npm/kmath/pull/140) diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt index 0ac6ea5f0..ec8714617 100644 --- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt +++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt @@ -29,7 +29,7 @@ class LinearAlgebraBenchmark { @Benchmark fun kmathLUPInversion() { - MatrixContext.real.inverse(matrix) + MatrixContext.real.inverseWithLUP(matrix) } @Benchmark diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt index 377d167bc..712927400 100644 --- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt +++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt @@ -5,8 +5,7 @@ import kscience.kmath.structures.Matrix import kscience.kmath.structures.NDStructure import org.apache.commons.math3.linear.* -public class CMMatrix(public val origin: RealMatrix, features: Set? = null) : - FeaturedMatrix { +public class CMMatrix(public val origin: RealMatrix, features: Set? = null) : FeaturedMatrix { public override val rowNum: Int get() = origin.rowDimension public override val colNum: Int get() = origin.columnDimension @@ -55,7 +54,7 @@ public fun Point.toCM(): CMVector = if (this is CMVector) this else { public fun RealVector.toPoint(): CMVector = CMVector(this) -public object CMMatrixContext : MatrixContext { +public object CMMatrixContext : MatrixContext { public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): CMMatrix { val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } } return CMMatrix(Array2DRowRealMatrix(array)) @@ -79,7 +78,7 @@ public object CMMatrixContext : MatrixContext { public override fun multiply(a: Matrix, k: Number): CMMatrix = CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble())) - public override operator fun Matrix.times(value: Double): Matrix = + public override operator fun Matrix.times(value: Double): CMMatrix = produce(rowNum, colNum) { i, j -> get(i, j) * value } } diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt index f6ed34be0..8b50bbe33 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt @@ -10,7 +10,7 @@ import kscience.kmath.structures.* public class BufferMatrixContext>( public override val elementContext: R, private val bufferFactory: BufferFactory, -) : GenericMatrixContext { +) : GenericMatrixContext> { public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): BufferMatrix { val buffer = bufferFactory(rows * columns) { offset -> initializer(offset / columns, offset % columns) } return BufferMatrix(rows, columns, buffer) @@ -22,7 +22,7 @@ public class BufferMatrixContext>( } @Suppress("OVERRIDE_BY_INLINE") -public object RealMatrixContext : GenericMatrixContext { +public object RealMatrixContext : GenericMatrixContext> { public override val elementContext: RealField get() = RealField diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt index 350379717..68272203c 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt @@ -57,7 +57,7 @@ public inline fun Matrix<*>.getFeature(): T? = /** * Diagonal matrix of ones. The matrix is virtual no actual matrix is created */ -public fun > GenericMatrixContext.one(rows: Int, columns: Int): FeaturedMatrix = +public fun > GenericMatrixContext.one(rows: Int, columns: Int): FeaturedMatrix = VirtualMatrix(rows, columns, DiagonalFeature) { i, j -> if (i == j) elementContext.one else elementContext.zero } @@ -66,7 +66,7 @@ public fun > GenericMatrixContext.one(rows: Int, colu /** * A virtual matrix of zeroes */ -public fun > GenericMatrixContext.zero(rows: Int, columns: Int): FeaturedMatrix = +public fun > GenericMatrixContext.zero(rows: Int, columns: Int): FeaturedMatrix = VirtualMatrix(rows, columns) { _, _ -> elementContext.zero } public class TransposedFeature(public val original: Matrix) : MatrixFeature diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt index bb80bcafd..099fa1909 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt @@ -1,25 +1,18 @@ package kscience.kmath.linear -import kscience.kmath.operations.Field -import kscience.kmath.operations.RealField -import kscience.kmath.operations.Ring -import kscience.kmath.operations.invoke -import kscience.kmath.structures.BufferAccessor2D -import kscience.kmath.structures.Matrix -import kscience.kmath.structures.Structure2D -import kotlin.reflect.KClass +import kscience.kmath.operations.* +import kscience.kmath.structures.* /** * Common implementation of [LUPDecompositionFeature] */ public class LUPDecomposition( - public val context: GenericMatrixContext>, + public val context: MatrixContext>, + public val elementContext: Field, public val lu: Structure2D, public val pivot: IntArray, - private val even: Boolean + private val even: Boolean, ) : LUPDecompositionFeature, DeterminantFeature { - public val elementContext: Field - get() = context.elementContext /** * Returns the matrix L of the decomposition. @@ -64,23 +57,25 @@ public class LUPDecomposition( } -public fun , F : Field> GenericMatrixContext.abs(value: T): T = +@PublishedApi +internal fun , F : Field> GenericMatrixContext.abs(value: T): T = if (value > elementContext.zero) value else elementContext { -value } /** - * Create a lup decomposition of generic matrix + * Create a lup decomposition of generic matrix. */ -public inline fun , F : Field> GenericMatrixContext.lup( - type: KClass, +public fun > MatrixContext>.lup( + factory: MutableBufferFactory, + elementContext: Field, matrix: Matrix, - checkSingular: (T) -> Boolean + checkSingular: (T) -> Boolean, ): LUPDecomposition { require(matrix.rowNum == matrix.colNum) { "LU decomposition supports only square matrices" } val m = matrix.colNum val pivot = IntArray(matrix.rowNum) //TODO just waits for KEEP-176 - BufferAccessor2D(type, matrix.rowNum, matrix.colNum).run { + BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run { elementContext { val lu = create(matrix) @@ -112,14 +107,14 @@ public inline fun , F : Field> GenericMatrixContext.l luRow[col] = sum // maintain best permutation choice - if (this@lup.abs(sum) > largest) { - largest = this@lup.abs(sum) + if (abs(sum) > largest) { + largest = abs(sum) max = row } } // Singularity check - check(!checkSingular(this@lup.abs(lu[max, col]))) { "The matrix is singular" } + check(!checkSingular(abs(lu[max, col]))) { "The matrix is singular" } // Pivot if necessary if (max != col) { @@ -143,23 +138,23 @@ public inline fun , F : Field> GenericMatrixContext.l for (row in col + 1 until m) lu[row, col] /= luDiag } - return LUPDecomposition(this@lup, lu.collect(), pivot, even) + return LUPDecomposition(this@lup, elementContext, lu.collect(), pivot, even) } } } -public inline fun , F : Field> GenericMatrixContext.lup( +public inline fun , F : Field> GenericMatrixContext>.lup( matrix: Matrix, - checkSingular: (T) -> Boolean -): LUPDecomposition = lup(T::class, matrix, checkSingular) + noinline checkSingular: (T) -> Boolean, +): LUPDecomposition = lup(MutableBuffer.Companion::auto, elementContext, matrix, checkSingular) -public fun GenericMatrixContext.lup(matrix: Matrix): LUPDecomposition = - lup(Double::class, matrix) { it < 1e-11 } +public fun MatrixContext>.lup(matrix: Matrix): LUPDecomposition = + lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 } -public fun LUPDecomposition.solve(type: KClass, matrix: Matrix): Matrix { +public fun LUPDecomposition.solveWithLUP(factory: MutableBufferFactory, matrix: Matrix): FeaturedMatrix { require(matrix.rowNum == pivot.size) { "Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}" } - BufferAccessor2D(type, matrix.rowNum, matrix.colNum).run { + BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run { elementContext { // Apply permutations to b val bp = create { _, _ -> zero } @@ -201,27 +196,34 @@ public fun LUPDecomposition.solve(type: KClass, matrix: Matrix LUPDecomposition.solve(matrix: Matrix): Matrix = solve(T::class, matrix) +public inline fun LUPDecomposition.solveWithLUP(matrix: Matrix): Matrix = + solveWithLUP(MutableBuffer.Companion::auto, matrix) /** - * Solve a linear equation **a*x = b** + * Solve a linear equation **a*x = b** using LUP decomposition */ -public inline fun , F : Field> GenericMatrixContext.solve( +public inline fun , F : Field> GenericMatrixContext>.solveWithLUP( a: Matrix, b: Matrix, - checkSingular: (T) -> Boolean -): Matrix { + noinline bufferFactory: MutableBufferFactory = MutableBuffer.Companion::auto, + noinline checkSingular: (T) -> Boolean, +): FeaturedMatrix { // Use existing decomposition if it is provided by matrix - val decomposition = a.getFeature() ?: lup(T::class, a, checkSingular) - return decomposition.solve(T::class, b) + val decomposition = a.getFeature() ?: lup(bufferFactory, elementContext, a, checkSingular) + return decomposition.solveWithLUP(bufferFactory, b) } -public fun RealMatrixContext.solve(a: Matrix, b: Matrix): Matrix = solve(a, b) { it < 1e-11 } +public fun RealMatrixContext.solveWithLUP(a: Matrix, b: Matrix): FeaturedMatrix = + solveWithLUP(a, b) { it < 1e-11 } -public inline fun , F : Field> GenericMatrixContext.inverse( +public inline fun , F : Field> GenericMatrixContext>.inverseWithLUP( matrix: Matrix, - checkSingular: (T) -> Boolean -): Matrix = solve(matrix, one(matrix.rowNum, matrix.colNum), checkSingular) + noinline bufferFactory: MutableBufferFactory = MutableBuffer.Companion::auto, + noinline checkSingular: (T) -> Boolean, +): FeaturedMatrix = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular) -public fun RealMatrixContext.inverse(matrix: Matrix): Matrix = - solve(matrix, one(matrix.rowNum, matrix.colNum)) { it < 1e-11 } +/** + * Inverses a square matrix using LUP decomposition. Non square matrix will throw a error. + */ +public fun RealMatrixContext.inverseWithLUP(matrix: Matrix): FeaturedMatrix = + solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), Buffer.Companion::real) { it < 1e-11 } diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt index f4dbce89a..d9dc57b0f 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt @@ -12,15 +12,16 @@ import kscience.kmath.structures.asSequence /** * Basic operations on matrices. Operates on [Matrix] */ -public interface MatrixContext : SpaceOperations> { +public interface MatrixContext> : SpaceOperations> { /** * Produce a matrix with this context and given dimensions */ - public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix + public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): M - public override fun binaryOperation(operation: String, left: Matrix, right: Matrix): Matrix = when (operation) { + @Suppress("UNCHECKED_CAST") + public override fun binaryOperation(operation: String, left: Matrix, right: Matrix): M = when (operation) { "dot" -> left dot right - else -> super.binaryOperation(operation, left, right) + else -> super.binaryOperation(operation, left, right) as M } /** @@ -30,7 +31,7 @@ public interface MatrixContext : SpaceOperations> { * @param other the multiplier. * @return the dot product. */ - public infix fun Matrix.dot(other: Matrix): Matrix + public infix fun Matrix.dot(other: Matrix): M /** * Computes the dot product of this matrix and a vector. @@ -48,7 +49,7 @@ public interface MatrixContext : SpaceOperations> { * @param value the multiplier. * @receiver the product. */ - public operator fun Matrix.times(value: T): Matrix + public operator fun Matrix.times(value: T): M /** * Multiplies an element by a matrix of it. @@ -57,7 +58,7 @@ public interface MatrixContext : SpaceOperations> { * @param value the multiplier. * @receiver the product. */ - public operator fun T.times(m: Matrix): Matrix = m * this + public operator fun T.times(m: Matrix): M = m * this public companion object { /** @@ -70,18 +71,18 @@ public interface MatrixContext : SpaceOperations> { */ public fun > buffered( ring: R, - bufferFactory: BufferFactory = Buffer.Companion::boxing - ): GenericMatrixContext = BufferMatrixContext(ring, bufferFactory) + bufferFactory: BufferFactory = Buffer.Companion::boxing, + ): GenericMatrixContext> = BufferMatrixContext(ring, bufferFactory) /** * Automatic buffered matrix, unboxed if it is possible */ - public inline fun > auto(ring: R): GenericMatrixContext = + public inline fun > auto(ring: R): GenericMatrixContext> = buffered(ring, Buffer.Companion::auto) } } -public interface GenericMatrixContext> : MatrixContext { +public interface GenericMatrixContext, out M : Matrix> : MatrixContext { /** * The ring context for matrix elements */ @@ -92,7 +93,7 @@ public interface GenericMatrixContext> : MatrixContext { */ public fun point(size: Int, initializer: (Int) -> T): Point - public override infix fun Matrix.dot(other: Matrix): Matrix { + public override infix fun Matrix.dot(other: Matrix): M { //TODO add typed error require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" } @@ -113,10 +114,10 @@ public interface GenericMatrixContext> : MatrixContext { } } - public override operator fun Matrix.unaryMinus(): Matrix = + public override operator fun Matrix.unaryMinus(): M = produce(rowNum, colNum) { i, j -> elementContext { -get(i, j) } } - public override fun add(a: Matrix, b: Matrix): Matrix { + public override fun add(a: Matrix, b: Matrix): M { require(a.rowNum == b.rowNum && a.colNum == b.colNum) { "Matrix operation dimension mismatch. [${a.rowNum},${a.colNum}] + [${b.rowNum},${b.colNum}]" } @@ -124,7 +125,7 @@ public interface GenericMatrixContext> : MatrixContext { return produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] + b[i, j] } } } - public override operator fun Matrix.minus(b: Matrix): Matrix { + public override operator fun Matrix.minus(b: Matrix): M { require(rowNum == b.rowNum && colNum == b.colNum) { "Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]" } @@ -132,11 +133,11 @@ public interface GenericMatrixContext> : MatrixContext { return produce(rowNum, colNum) { i, j -> elementContext { get(i, j) + b[i, j] } } } - public override fun multiply(a: Matrix, k: Number): Matrix = + public override fun multiply(a: Matrix, k: Number): M = produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] * k } } - public operator fun Number.times(matrix: FeaturedMatrix): Matrix = matrix * this + public operator fun Number.times(matrix: FeaturedMatrix): M = multiply(matrix, this) - public override operator fun Matrix.times(value: T): Matrix = + public override operator fun Matrix.times(value: T): M = produce(rowNum, colNum) { i, j -> elementContext { get(i, j) * value } } } diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt index 657da6b4e..4527a2a42 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt @@ -38,6 +38,11 @@ public fun Space.average(data: Iterable): T = sum(data) / data.count() */ public fun Space.average(data: Sequence): T = sum(data) / data.count() +/** + * Absolute of the comparable [value] + */ +public fun > Space.abs(value: T): T = if (value > zero) value else -value + /** * Returns the sum of all elements in the iterable in provided space. * diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt index 00fc4aef0..5d7ba611f 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt @@ -1,25 +1,31 @@ package kscience.kmath.structures -import kotlin.reflect.KClass - /** * A context that allows to operate on a [MutableBuffer] as on 2d array */ -public class BufferAccessor2D(public val type: KClass, public val rowNum: Int, public val colNum: Int) { +internal class BufferAccessor2D( + public val rowNum: Int, + public val colNum: Int, + val factory: MutableBufferFactory, +) { public operator fun Buffer.get(i: Int, j: Int): T = get(i + colNum * j) public operator fun MutableBuffer.set(i: Int, j: Int, value: T) { set(i + colNum * j, value) } - public inline fun create(init: (i: Int, j: Int) -> T): MutableBuffer = - MutableBuffer.auto(type, rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) } + public inline fun create(crossinline init: (i: Int, j: Int) -> T): MutableBuffer = + factory(rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) } public fun create(mat: Structure2D): MutableBuffer = create { i, j -> mat[i, j] } //TODO optimize wrapper - public fun MutableBuffer.collect(): Structure2D = - NDStructure.auto(type, rowNum, colNum) { (i, j) -> get(i, j) }.as2D() + public fun MutableBuffer.collect(): Structure2D = NDStructure.build( + DefaultStrides(intArrayOf(rowNum, colNum)), + factory + ) { (i, j) -> + get(i, j) + }.as2D() public inner class Row(public val buffer: MutableBuffer, public val rowIndex: Int) : MutableBuffer { override val size: Int get() = colNum @@ -30,7 +36,7 @@ public class BufferAccessor2D(public val type: KClass, public val ro buffer[rowIndex, index] = value } - override fun copy(): MutableBuffer = MutableBuffer.auto(type, colNum) { get(it) } + override fun copy(): MutableBuffer = factory(colNum) { get(it) } override operator fun iterator(): Iterator = (0 until colNum).map(::get).iterator() } diff --git a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt index de07d3639..28dfe46ec 100644 --- a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt +++ b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt @@ -9,7 +9,7 @@ class RealLUSolverTest { @Test fun testInvertOne() { val matrix = MatrixContext.real.one(2, 2) - val inverted = MatrixContext.real.inverse(matrix) + val inverted = MatrixContext.real.inverseWithLUP(matrix) assertEquals(matrix, inverted) } @@ -37,7 +37,7 @@ class RealLUSolverTest { 1.0, 3.0 ) - val inverted = MatrixContext.real.inverse(matrix) + val inverted = MatrixContext.real.inverseWithLUP(matrix) val expected = Matrix.square( 0.375, -0.125, diff --git a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt index 4111cb78d..68a5dc262 100644 --- a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt +++ b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt @@ -42,7 +42,7 @@ public interface DMatrix : Structure2D { * An inline wrapper for a Matrix */ public inline class DMatrixWrapper( - public val structure: Structure2D + private val structure: Structure2D ) : DMatrix { override val shape: IntArray get() = structure.shape override operator fun get(i: Int, j: Int): T = structure[i, j] @@ -81,7 +81,7 @@ public inline class DPointWrapper(public val point: Point) /** * Basic operations on dimension-safe matrices. Operates on [Matrix] */ -public inline class DMatrixContext>(public val context: GenericMatrixContext) { +public inline class DMatrixContext>(public val context: GenericMatrixContext>) { public inline fun Matrix.coerce(): DMatrix { require(rowNum == Dimension.dim().toInt()) { "Row number mismatch: expected ${Dimension.dim()} but found $rowNum" diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt index 1ed57fcb3..31792e39c 100644 --- a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt +++ b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt @@ -1,11 +1,9 @@ package kscience.kmath.ejml -import org.ejml.simple.SimpleMatrix import kscience.kmath.linear.MatrixContext import kscience.kmath.linear.Point -import kscience.kmath.operations.Space -import kscience.kmath.operations.invoke import kscience.kmath.structures.Matrix +import org.ejml.simple.SimpleMatrix /** * Converts this matrix to EJML one. @@ -18,7 +16,7 @@ public fun Matrix.toEjml(): EjmlMatrix = * * @author Iaroslav Postovalov */ -public object EjmlMatrixContext : MatrixContext { +public object EjmlMatrixContext : MatrixContext { /** * Converts this vector to EJML one. diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt index c0657c0bf..e8ad835e5 100644 --- a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt +++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt @@ -1,12 +1,16 @@ package kscience.kmath.real +import kscience.kmath.linear.FeaturedMatrix import kscience.kmath.linear.MatrixContext import kscience.kmath.linear.RealMatrixContext.elementContext import kscience.kmath.linear.VirtualMatrix +import kscience.kmath.linear.inverseWithLUP import kscience.kmath.misc.UnstableKMathAPI import kscience.kmath.operations.invoke import kscience.kmath.operations.sum -import kscience.kmath.structures.* +import kscience.kmath.structures.Buffer +import kscience.kmath.structures.RealBuffer +import kscience.kmath.structures.asIterable import kotlin.math.pow /* @@ -21,7 +25,7 @@ import kotlin.math.pow * Functions that help create a real (Double) matrix */ -public typealias RealMatrix = Matrix +public typealias RealMatrix = FeaturedMatrix public fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double): RealMatrix = MatrixContext.real.produce(rowNum, colNum, initializer) @@ -148,6 +152,11 @@ public inline fun RealMatrix.map(transform: (Double) -> Double): RealMatrix = transform(get(i, j)) } +/** + * Inverse a square real matrix using LUP decomposition + */ +public fun RealMatrix.inverseWithLUP(): RealMatrix = MatrixContext.real.inverseWithLUP(this) + //extended operations public fun RealMatrix.pow(p: Double): RealMatrix = map { it.pow(p) }