From 89eebbecb76280586f39646bc08cf6a6975bdae7 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Tue, 21 Sep 2021 21:24:27 +0300 Subject: [PATCH] Fix EJML inversion issue --- CHANGELOG.md | 1 + benchmarks/build.gradle.kts | 10 +++++ .../kscience/kmath/benchmarks/DotBenchmark.kt | 11 +++-- .../benchmarks/MatrixInverseBenchmark.kt | 15 ++++--- .../kmath/benchmarks/NDFieldBenchmark.kt | 10 +++-- .../kmath/benchmarks/ViktorBenchmark.kt | 8 ++-- .../kmath/benchmarks/ViktorLogBenchmark.kt | 8 ++-- .../kmath/functions/matrixIntegration.kt | 8 ++-- .../kscience/kmath/operations/ComplexDemo.kt | 31 -------------- .../kscience/kmath/operations/complexDemo.kt | 41 +++++++++++++++++++ .../kscience/kmath/structures/ComplexND.kt | 8 ++-- .../kscience/kmath/structures/NDField.kt | 10 ++--- .../kscience/kmath/structures/buffers.kt | 4 +- .../kscience/kmath/commons/linear/CMMatrix.kt | 6 ++- .../kscience/kmath/commons/linear/CMSolver.kt | 12 +++--- .../kscience/kmath/complex/ComplexFieldND.kt | 9 +++- .../kmath/linear/BufferedLinearSpace.kt | 19 +++++---- .../kscience/kmath/linear/LinearSpace.kt | 21 ++++++---- .../kscience/kmath/linear/LupDecomposition.kt | 6 +-- .../kscience/kmath/linear/MatrixWrapper.kt | 9 ++-- .../kscience/kmath/nd/BufferAlgebraND.kt | 20 ++++----- .../space/kscience/kmath/nd/DoubleFieldND.kt | 4 +- .../space/kscience/kmath/nd/ShortRingND.kt | 2 +- .../space/kscience/kmath/nd/Structure1D.kt | 6 ++- .../space/kscience/kmath/nd/Structure2D.kt | 2 - .../space/kscience/kmath/nd/StructureND.kt | 23 +++++++++-- .../kmath/operations/BufferAlgebra.kt | 24 +++++++---- .../kscience/kmath/structures/ArrayBuffer.kt | 2 + .../space/kscience/kmath/structures/Buffer.kt | 7 ++++ .../kmath/structures/BufferAccessor2D.kt | 2 + .../kmath/structures/FlaggedBuffer.kt | 8 +++- .../kscience/kmath/structures/ListBuffer.kt | 2 + .../kscience/kmath/structures/MemoryBuffer.kt | 2 + .../kscience/kmath/structures/NDFieldTest.kt | 6 +-- .../kmath/structures/NumberNDFieldTest.kt | 6 +-- .../kscience/kmath/streaming/RingBuffer.kt | 2 + .../kscience/kmath/ejml/EjmlLinearSpace.kt | 6 +++ .../kscience/kmath/ejml/EjmlMatrixTest.kt | 25 +++++++++-- .../space/kscience/kmath/real/RealMatrix.kt | 37 +++++++++-------- .../kotlin/space/kscience/kmath/real/dot.kt | 5 ++- .../kmath/integration/SplineIntegrator.kt | 2 + .../kmath/histogram/DoubleHistogramSpace.kt | 2 +- .../kscience/kmath/viktor/ViktorBuffer.kt | 6 ++- 43 files changed, 284 insertions(+), 164 deletions(-) delete mode 100644 examples/src/main/kotlin/space/kscience/kmath/operations/ComplexDemo.kt create mode 100644 examples/src/main/kotlin/space/kscience/kmath/operations/complexDemo.kt diff --git a/CHANGELOG.md b/CHANGELOG.md index 19b7dae32..c27dafeb4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,6 +40,7 @@ - `FeatureSet` now accepts only `Feature`. It is possible to override keys and use interfaces. - Use `Symbol` factory function instead of `StringSymbol` - New discoverability pattern: `.algebra.` +- Adjusted commons-math API for linear solvers to match conventions. ### Deprecated - Specialized `DoubleBufferAlgebra` diff --git a/benchmarks/build.gradle.kts b/benchmarks/build.gradle.kts index d96c5a8b6..fa72dc8ee 100644 --- a/benchmarks/build.gradle.kts +++ b/benchmarks/build.gradle.kts @@ -105,6 +105,16 @@ benchmark { commonConfiguration() include("JafamaBenchmark") } + + configurations.register("viktor") { + commonConfiguration() + include("ViktorBenchmark") + } + + configurations.register("viktorLog") { + commonConfiguration() + include("ViktorLogBenchmark") + } } // Fix kotlinx-benchmarks bug diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt index 7806cd06c..6373d1ce5 100644 --- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt @@ -13,7 +13,10 @@ import space.kscience.kmath.commons.linear.CMLinearSpace import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM import space.kscience.kmath.linear.LinearSpace import space.kscience.kmath.linear.invoke +import space.kscience.kmath.linear.linearSpace import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.operations.algebra +import space.kscience.kmath.structures.Buffer import kotlin.random.Random @State(Scope.Benchmark) @@ -35,7 +38,7 @@ internal class DotBenchmark { @Benchmark fun cmDot(blackhole: Blackhole) { - CMLinearSpace.run { + CMLinearSpace { blackhole.consume(cmMatrix1 dot cmMatrix2) } } @@ -56,14 +59,14 @@ internal class DotBenchmark { @Benchmark fun bufferedDot(blackhole: Blackhole) { - LinearSpace.auto(DoubleField).invoke { + with(DoubleField.linearSpace(Buffer.Companion::auto)) { blackhole.consume(matrix1 dot matrix2) } } @Benchmark - fun realDot(blackhole: Blackhole) { - LinearSpace.double { + fun doubleDot(blackhole: Blackhole) { + with(Double.algebra.linearSpace) { blackhole.consume(matrix1 dot matrix2) } } diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt index 21d7ca41c..5d331af9a 100644 --- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt @@ -10,13 +10,12 @@ import kotlinx.benchmark.Blackhole import kotlinx.benchmark.Scope import kotlinx.benchmark.State import space.kscience.kmath.commons.linear.CMLinearSpace -import space.kscience.kmath.commons.linear.inverse +import space.kscience.kmath.commons.linear.lupSolver import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM -import space.kscience.kmath.linear.InverseMatrixFeature -import space.kscience.kmath.linear.LinearSpace import space.kscience.kmath.linear.invoke +import space.kscience.kmath.linear.linearSpace import space.kscience.kmath.linear.lupSolver -import space.kscience.kmath.nd.getFeature +import space.kscience.kmath.operations.algebra import kotlin.random.Random @State(Scope.Benchmark) @@ -25,7 +24,7 @@ internal class MatrixInverseBenchmark { private val random = Random(1224) private const val dim = 100 - private val space = LinearSpace.double + private val space = Double.algebra.linearSpace //creating invertible matrix private val u = space.buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } @@ -35,20 +34,20 @@ internal class MatrixInverseBenchmark { @Benchmark fun kmathLupInversion(blackhole: Blackhole) { - blackhole.consume(LinearSpace.double.lupSolver().inverse(matrix)) + blackhole.consume(Double.algebra.linearSpace.lupSolver().inverse(matrix)) } @Benchmark fun cmLUPInversion(blackhole: Blackhole) { CMLinearSpace { - blackhole.consume(inverse(matrix)) + blackhole.consume(lupSolver().inverse(matrix)) } } @Benchmark fun ejmlInverse(blackhole: Blackhole) { EjmlLinearSpaceDDRM { - blackhole.consume(matrix.getFeature>()?.inverse) + blackhole.consume(matrix.toEjml().inverse()) } } } diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/NDFieldBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/NDFieldBenchmark.kt index fbf4bdaba..f72bc3ba0 100644 --- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/NDFieldBenchmark.kt +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/NDFieldBenchmark.kt @@ -9,7 +9,9 @@ import kotlinx.benchmark.Benchmark import kotlinx.benchmark.Blackhole import kotlinx.benchmark.Scope import kotlinx.benchmark.State -import space.kscience.kmath.nd.* +import space.kscience.kmath.nd.StructureND +import space.kscience.kmath.nd.autoNdAlgebra +import space.kscience.kmath.nd.ndAlgebra import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.structures.Buffer @@ -46,8 +48,8 @@ internal class NDFieldBenchmark { private companion object { private const val dim = 1000 private const val n = 100 - private val autoField = DoubleField.autoNd(dim, dim) - private val specializedField = DoubleField.nd(dim, dim) - private val genericField = DoubleField.nd(Buffer.Companion::boxing, dim, dim) + private val autoField = DoubleField.autoNdAlgebra(dim, dim) + private val specializedField = DoubleField.ndAlgebra(dim, dim) + private val genericField = DoubleField.ndAlgebra(Buffer.Companion::boxing, dim, dim) } } diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorBenchmark.kt index 6ea64ba3c..b97a05a52 100644 --- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorBenchmark.kt +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorBenchmark.kt @@ -11,8 +11,8 @@ import kotlinx.benchmark.Scope import kotlinx.benchmark.State import org.jetbrains.bio.viktor.F64Array import space.kscience.kmath.nd.StructureND -import space.kscience.kmath.nd.autoNd -import space.kscience.kmath.nd.nd +import space.kscience.kmath.nd.autoNdAlgebra +import space.kscience.kmath.nd.ndAlgebra import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.viktor.ViktorNDField @@ -58,8 +58,8 @@ internal class ViktorBenchmark { private const val n = 100 // automatically build context most suited for given type. - private val autoField = DoubleField.autoNd(dim, dim) - private val realField = DoubleField.nd(dim, dim) + private val autoField = DoubleField.autoNdAlgebra(dim, dim) + private val realField = DoubleField.ndAlgebra(dim, dim) private val viktorField = ViktorNDField(dim, dim) } } diff --git a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorLogBenchmark.kt b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorLogBenchmark.kt index afe44ea99..91e9dcd76 100644 --- a/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorLogBenchmark.kt +++ b/benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/ViktorLogBenchmark.kt @@ -10,8 +10,8 @@ import kotlinx.benchmark.Blackhole import kotlinx.benchmark.Scope import kotlinx.benchmark.State import org.jetbrains.bio.viktor.F64Array -import space.kscience.kmath.nd.autoNd -import space.kscience.kmath.nd.nd +import space.kscience.kmath.nd.autoNdAlgebra +import space.kscience.kmath.nd.ndAlgebra import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.viktor.ViktorFieldND @@ -50,8 +50,8 @@ internal class ViktorLogBenchmark { private const val n = 100 // automatically build context most suited for given type. - private val autoField = DoubleField.autoNd(dim, dim) - private val realNdField = DoubleField.nd(dim, dim) + private val autoField = DoubleField.autoNdAlgebra(dim, dim) + private val realNdField = DoubleField.ndAlgebra(dim, dim) private val viktorField = ViktorFieldND(intArrayOf(dim, dim)) } } diff --git a/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt b/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt index d932fdb9f..93b5671fe 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/functions/matrixIntegration.kt @@ -9,12 +9,12 @@ import space.kscience.kmath.integration.gaussIntegrator import space.kscience.kmath.integration.integrate import space.kscience.kmath.integration.value import space.kscience.kmath.nd.StructureND -import space.kscience.kmath.nd.withNd -import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.nd.withNdAlgebra +import space.kscience.kmath.operations.algebra import space.kscience.kmath.operations.invoke -fun main(): Unit = DoubleField { - withNd(2, 2) { +fun main(): Unit = Double.algebra { + withNdAlgebra(2, 2) { //Produce a diagonal StructureND fun diagonal(v: Double) = produce { (i, j) -> diff --git a/examples/src/main/kotlin/space/kscience/kmath/operations/ComplexDemo.kt b/examples/src/main/kotlin/space/kscience/kmath/operations/ComplexDemo.kt deleted file mode 100644 index 743f05b13..000000000 --- a/examples/src/main/kotlin/space/kscience/kmath/operations/ComplexDemo.kt +++ /dev/null @@ -1,31 +0,0 @@ -/* - * 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.complex.Complex -import space.kscience.kmath.complex.ComplexField -import space.kscience.kmath.complex.nd -import space.kscience.kmath.complex.withNd -import space.kscience.kmath.nd.StructureND - -fun main() { - // 2d element - val element = ComplexField.nd(2, 2).produce { (i, j) -> - Complex(i - j, i + j) - } - println(element) - - // 1d element operation - val result: StructureND = ComplexField.withNd(8) { - val a = produce { (it) -> i * it - it.toDouble() } - val b = 3 - val c = Complex(1.0, 1.0) - - (a pow b) + c - } - - println(result) -} diff --git a/examples/src/main/kotlin/space/kscience/kmath/operations/complexDemo.kt b/examples/src/main/kotlin/space/kscience/kmath/operations/complexDemo.kt new file mode 100644 index 000000000..319221bcc --- /dev/null +++ b/examples/src/main/kotlin/space/kscience/kmath/operations/complexDemo.kt @@ -0,0 +1,41 @@ +/* + * 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.complex.Complex +import space.kscience.kmath.complex.algebra +import space.kscience.kmath.complex.bufferAlgebra +import space.kscience.kmath.complex.ndAlgebra +import space.kscience.kmath.nd.BufferND +import space.kscience.kmath.nd.StructureND + +fun main() = Complex.algebra { + val complex = 2 + 2 * i + println(complex * 8 - 5 * i) + + //flat buffer + val buffer = bufferAlgebra(8).run { + buffer { Complex(it, -it) }.map { Complex(it.im, it.re) } + } + println(buffer) + + + // 2d element + val element: BufferND = ndAlgebra(2, 2).produce { (i, j) -> + Complex(i - j, i + j) + } + println(element) + + // 1d element operation + val result: StructureND = ndAlgebra(8).run { + val a = produce { (it) -> i * it - it.toDouble() } + val b = 3 + val c = Complex(1.0, 1.0) + + (a pow b) + c + } + println(result) +} diff --git a/examples/src/main/kotlin/space/kscience/kmath/structures/ComplexND.kt b/examples/src/main/kotlin/space/kscience/kmath/structures/ComplexND.kt index c591f5682..d4554b3ba 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/structures/ComplexND.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/structures/ComplexND.kt @@ -11,7 +11,7 @@ import space.kscience.kmath.complex.* import space.kscience.kmath.linear.transpose import space.kscience.kmath.nd.StructureND import space.kscience.kmath.nd.as2D -import space.kscience.kmath.nd.nd +import space.kscience.kmath.nd.ndAlgebra import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.invoke import kotlin.system.measureTimeMillis @@ -20,8 +20,8 @@ fun main() { val dim = 1000 val n = 1000 - val realField = DoubleField.nd(dim, dim) - val complexField: ComplexFieldND = ComplexField.nd(dim, dim) + val realField = DoubleField.ndAlgebra(dim, dim) + val complexField: ComplexFieldND = ComplexField.ndAlgebra(dim, dim) val realTime = measureTimeMillis { realField { @@ -49,7 +49,7 @@ fun main() { fun complexExample() { //Create a context for 2-d structure with complex values ComplexField { - withNd(4, 8) { + withNdAlgebra(4, 8) { //a constant real-valued structure val x = one * 2.5 operator fun Number.plus(other: Complex) = Complex(this.toDouble() + other.re, other.im) diff --git a/examples/src/main/kotlin/space/kscience/kmath/structures/NDField.kt b/examples/src/main/kotlin/space/kscience/kmath/structures/NDField.kt index c9cafbba8..5b0e2eb30 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/structures/NDField.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/structures/NDField.kt @@ -9,8 +9,8 @@ import kotlinx.coroutines.DelicateCoroutinesApi import kotlinx.coroutines.GlobalScope import org.nd4j.linalg.factory.Nd4j import space.kscience.kmath.nd.StructureND -import space.kscience.kmath.nd.autoNd -import space.kscience.kmath.nd.nd +import space.kscience.kmath.nd.autoNdAlgebra +import space.kscience.kmath.nd.ndAlgebra import space.kscience.kmath.nd4j.Nd4jArrayField import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.invoke @@ -33,11 +33,11 @@ fun main() { val n = 1000 // automatically build context most suited for given type. - val autoField = DoubleField.autoNd(dim, dim) + val autoField = DoubleField.autoNdAlgebra(dim, dim) // specialized nd-field for Double. It works as generic Double field as well. - val realField = DoubleField.nd(dim, dim) + val realField = DoubleField.ndAlgebra(dim, dim) //A generic boxing field. It should be used for objects, not primitives. - val boxingField = DoubleField.nd(Buffer.Companion::boxing, dim, dim) + val boxingField = DoubleField.ndAlgebra(Buffer.Companion::boxing, dim, dim) // Nd4j specialized field. val nd4jField = Nd4jArrayField.real(dim, dim) //viktor field diff --git a/examples/src/main/kotlin/space/kscience/kmath/structures/buffers.kt b/examples/src/main/kotlin/space/kscience/kmath/structures/buffers.kt index 0f7fdc46a..d78141507 100644 --- a/examples/src/main/kotlin/space/kscience/kmath/structures/buffers.kt +++ b/examples/src/main/kotlin/space/kscience/kmath/structures/buffers.kt @@ -6,8 +6,8 @@ package space.kscience.kmath.structures import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.operations.buffer import space.kscience.kmath.operations.bufferAlgebra -import space.kscience.kmath.operations.produce inline fun MutableBuffer.Companion.same( n: Int, @@ -17,6 +17,6 @@ inline fun MutableBuffer.Companion.same( fun main() { with(DoubleField.bufferAlgebra(5)) { - println(number(2.0) + produce(1, 2, 3, 4, 5)) + println(number(2.0) + buffer(1, 2, 3, 4, 5)) } } diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt index 7defd6841..14e7fc365 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt @@ -10,6 +10,7 @@ import space.kscience.kmath.linear.* import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.nd.StructureFeature import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.DoubleBuffer import kotlin.reflect.KClass import kotlin.reflect.cast @@ -21,12 +22,15 @@ public class CMMatrix(public val origin: RealMatrix) : Matrix { override operator fun get(i: Int, j: Int): Double = origin.getEntry(i, j) } -public class CMVector(public val origin: RealVector) : Point { +@JvmInline +public value class CMVector(public val origin: RealVector) : Point { override val size: Int get() = origin.dimension override operator fun get(index: Int): Double = origin.getEntry(index) override operator fun iterator(): Iterator = origin.toArray().iterator() + + override fun toString(): String = Buffer.toString(this) } public fun RealVector.toPoint(): CMVector = CMVector(this) diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMSolver.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMSolver.kt index b9fbc07ae..d1fb441b0 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMSolver.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMSolver.kt @@ -18,7 +18,7 @@ public enum class CMDecomposition { CHOLESKY } -public fun CMLinearSpace.solver( +private fun CMLinearSpace.solver( a: Matrix, decomposition: CMDecomposition = CMDecomposition.LUP, ): DecompositionSolver = when (decomposition) { @@ -48,9 +48,11 @@ public fun CMLinearSpace.inverse( public fun CMLinearSpace.solver(decomposition: CMDecomposition): LinearSolver = object : LinearSolver { - override fun solve(a: Matrix, b: Matrix): Matrix = solve(a, b, decomposition) + override fun solve(a: Matrix, b: Matrix): Matrix = solver(a, decomposition).solve(b.toCM().origin).wrap() - override fun solve(a: Matrix, b: Point): Point = solve(a, b, decomposition) + override fun solve(a: Matrix, b: Point): Point = solver(a, decomposition).solve(b.toCM().origin).toPoint() - override fun inverse(matrix: Matrix): Matrix = inverse(matrix, decomposition) -} \ No newline at end of file + override fun inverse(matrix: Matrix): Matrix = solver(matrix, decomposition).inverse.wrap() +} + +public fun CMLinearSpace.lupSolver(): LinearSolver = solver((CMDecomposition.LUP)) \ No newline at end of file diff --git a/kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/ComplexFieldND.kt b/kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/ComplexFieldND.kt index 41bcb83df..29e790d16 100644 --- a/kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/ComplexFieldND.kt +++ b/kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/ComplexFieldND.kt @@ -9,8 +9,10 @@ import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.nd.BufferND import space.kscience.kmath.nd.BufferedFieldND import space.kscience.kmath.nd.StructureND +import space.kscience.kmath.operations.BufferField import space.kscience.kmath.operations.ExtendedField import space.kscience.kmath.operations.NumbersAddOperations +import space.kscience.kmath.operations.bufferAlgebra import space.kscience.kmath.structures.Buffer import kotlin.contracts.InvocationKind import kotlin.contracts.contract @@ -111,13 +113,16 @@ public inline fun BufferedFieldND.produceInline(initializ return BufferND(strides, buffer) } +@UnstableKMathAPI +public fun ComplexField.bufferAlgebra(size: Int): BufferField = + bufferAlgebra(Buffer.Companion::complex, size) -public fun ComplexField.nd(vararg shape: Int): ComplexFieldND = ComplexFieldND(shape) +public fun ComplexField.ndAlgebra(vararg shape: Int): ComplexFieldND = ComplexFieldND(shape) /** * Produce a context for n-dimensional operations inside this real field */ -public inline fun ComplexField.withNd(vararg shape: Int, action: ComplexFieldND.() -> R): R { +public inline fun ComplexField.withNdAlgebra(vararg shape: Int, action: ComplexFieldND.() -> R): R { contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) } return ComplexFieldND(shape).action() } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/BufferedLinearSpace.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/BufferedLinearSpace.kt index 14235515c..5f17de607 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/BufferedLinearSpace.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/BufferedLinearSpace.kt @@ -8,17 +8,15 @@ package space.kscience.kmath.linear import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.nd.BufferedRingND import space.kscience.kmath.nd.as2D -import space.kscience.kmath.nd.nd +import space.kscience.kmath.nd.ndAlgebra import space.kscience.kmath.nd.unwrap +import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.Ring import space.kscience.kmath.operations.invoke -import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.BufferFactory -import space.kscience.kmath.structures.VirtualBuffer -import space.kscience.kmath.structures.indices +import space.kscience.kmath.structures.* -public class BufferedLinearSpace>( +public class BufferedLinearSpace>( override val elementAlgebra: A, private val bufferFactory: BufferFactory, ) : LinearSpace { @@ -26,7 +24,7 @@ public class BufferedLinearSpace>( private fun ndRing( rows: Int, cols: Int, - ): BufferedRingND = elementAlgebra.nd(bufferFactory, rows, cols) + ): BufferedRingND = elementAlgebra.ndAlgebra(bufferFactory, rows, cols) override fun buildMatrix(rows: Int, columns: Int, initializer: A.(i: Int, j: Int) -> T): Matrix = ndRing(rows, columns).produce { (i, j) -> elementAlgebra.initializer(i, j) }.as2D() @@ -92,3 +90,10 @@ public class BufferedLinearSpace>( unwrap().map { it * value }.as2D() } } + + +public fun > A.linearSpace(bufferFactory: BufferFactory): BufferedLinearSpace = + BufferedLinearSpace(this, bufferFactory) + +public val DoubleField.linearSpace: BufferedLinearSpace + get() = BufferedLinearSpace(this, ::DoubleBuffer) \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt index e94e984b3..1d8985b59 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt @@ -6,8 +6,13 @@ package space.kscience.kmath.linear import space.kscience.kmath.misc.UnstableKMathAPI -import space.kscience.kmath.nd.* -import space.kscience.kmath.operations.* +import space.kscience.kmath.nd.MutableStructure2D +import space.kscience.kmath.nd.Structure2D +import space.kscience.kmath.nd.StructureFeature +import space.kscience.kmath.nd.as1D +import space.kscience.kmath.operations.DoubleField +import space.kscience.kmath.operations.Ring +import space.kscience.kmath.operations.invoke import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.BufferFactory import space.kscience.kmath.structures.DoubleBuffer @@ -34,7 +39,7 @@ public typealias Point = Buffer * @param T the type of items in the matrices. * @param A the type of ring over [T]. */ -public interface LinearSpace> { +public interface LinearSpace> { public val elementAlgebra: A /** @@ -172,7 +177,8 @@ public interface LinearSpace> { * @return a feature object or `null` if it isn't present. */ @UnstableKMathAPI - public fun computeFeature(structure: Matrix, type: KClass): F? = structure.getFeature(type) + public fun computeFeature(structure: Matrix, type: KClass): F? = + structure.getFeature(type) public companion object { @@ -184,6 +190,7 @@ public interface LinearSpace> { bufferFactory: BufferFactory = Buffer.Companion::boxing, ): LinearSpace = BufferedLinearSpace(algebra, bufferFactory) + @Deprecated("use DoubleField.linearSpace") public val double: LinearSpace = buffered(DoubleField, ::DoubleBuffer) /** @@ -213,10 +220,8 @@ public inline operator fun , R> LS.invoke(block: LS.() -> * Convert matrix to vector if it is possible. */ public fun Matrix.asVector(): Point = - if (this.colNum == 1) - as1D() - else - error("Can't convert matrix with more than one column to vector") + if (this.colNum == 1) as1D() + else error("Can't convert matrix with more than one column to vector") /** * Creates an n × 1 [VirtualMatrix], where n is the size of the given buffer. diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt index 34fc661ed..95dd6d45c 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt @@ -5,9 +5,7 @@ package space.kscience.kmath.linear -import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.misc.UnstableKMathAPI -import space.kscience.kmath.nd.getFeature import space.kscience.kmath.operations.* import space.kscience.kmath.structures.BufferAccessor2D import space.kscience.kmath.structures.DoubleBuffer @@ -214,7 +212,6 @@ internal fun LupDecomposition.solve( /** * Produce a generic solver based on LUP decomposition */ -@PerformancePitfall() @OptIn(UnstableKMathAPI::class) public fun , F : Field> LinearSpace.lupSolver( bufferFactory: MutableBufferFactory, @@ -222,13 +219,12 @@ public fun , F : Field> LinearSpace.lupSolver( ): LinearSolver = object : LinearSolver { override fun solve(a: Matrix, b: Matrix): Matrix { // Use existing decomposition if it is provided by matrix - val decomposition = a.getFeature() ?: lup(bufferFactory, a, singularityCheck) + val decomposition = computeFeature(a) ?: lup(bufferFactory, a, singularityCheck) return decomposition.solve(bufferFactory, b) } override fun inverse(matrix: Matrix): Matrix = solve(matrix, one(matrix.rowNum, matrix.colNum)) } -@PerformancePitfall public fun LinearSpace.lupSolver(singularityThreshold: Double = 1e-11): LinearSolver = lupSolver(::DoubleBuffer) { it < singularityThreshold } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt index 3ba4e8648..a40c0384c 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt @@ -8,7 +8,6 @@ package space.kscience.kmath.linear import space.kscience.kmath.misc.FeatureSet import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.nd.StructureFeature -import space.kscience.kmath.nd.getFeature import space.kscience.kmath.operations.Ring import kotlin.reflect.KClass @@ -26,7 +25,6 @@ public class MatrixWrapper internal constructor( * Get the first feature matching given class. Does not guarantee that matrix has only one feature matching the * criteria. */ - @UnstableKMathAPI @Suppress("UNCHECKED_CAST") override fun getFeature(type: KClass): F? = features.getFeature(type) ?: origin.getFeature(type) @@ -90,8 +88,7 @@ public class TransposedFeature(public val original: Matrix) : Ma /** * Create a virtual transposed matrix without copying anything. `A.transpose().transpose() === A` */ +@Suppress("UNCHECKED_CAST") @OptIn(UnstableKMathAPI::class) -public fun Matrix.transpose(): Matrix = getFeature>()?.original ?: VirtualMatrix( - colNum, - rowNum, -) { i, j -> get(j, i) }.withFeature(TransposedFeature(this)) +public fun Matrix.transpose(): Matrix = getFeature(TransposedFeature::class)?.original as? Matrix + ?: VirtualMatrix(colNum, rowNum) { i, j -> get(j, i) }.withFeature(TransposedFeature(this)) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt index b01beb02d..ae72f3689 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt @@ -87,39 +87,39 @@ public open class BufferedFieldND>( } // group factories -public fun > A.nd( +public fun > A.ndAlgebra( bufferFactory: BufferFactory, vararg shape: Int, ): BufferedGroupND = BufferedGroupND(shape, this, bufferFactory) @JvmName("withNdGroup") -public inline fun , R> A.withNd( +public inline fun , R> A.withNdAlgebra( noinline bufferFactory: BufferFactory, vararg shape: Int, action: BufferedGroupND.() -> R, ): R { contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) } - return nd(bufferFactory, *shape).run(action) + return ndAlgebra(bufferFactory, *shape).run(action) } //ring factories -public fun > A.nd( +public fun > A.ndAlgebra( bufferFactory: BufferFactory, vararg shape: Int, ): BufferedRingND = BufferedRingND(shape, this, bufferFactory) @JvmName("withNdRing") -public inline fun , R> A.withNd( +public inline fun , R> A.withNdAlgebra( noinline bufferFactory: BufferFactory, vararg shape: Int, action: BufferedRingND.() -> R, ): R { contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) } - return nd(bufferFactory, *shape).run(action) + return ndAlgebra(bufferFactory, *shape).run(action) } //field factories -public fun > A.nd( +public fun > A.ndAlgebra( bufferFactory: BufferFactory, vararg shape: Int, ): BufferedFieldND = BufferedFieldND(shape, this, bufferFactory) @@ -129,7 +129,7 @@ public fun > A.nd( */ @UnstableKMathAPI @Suppress("UNCHECKED_CAST") -public inline fun > A.autoNd( +public inline fun > A.autoNdAlgebra( vararg shape: Int, ): FieldND = when (this) { DoubleField -> DoubleFieldND(shape) as FieldND @@ -137,11 +137,11 @@ public inline fun > A.autoNd( } @JvmName("withNdField") -public inline fun , R> A.withNd( +public inline fun , R> A.withNdAlgebra( noinline bufferFactory: BufferFactory, vararg shape: Int, action: BufferedFieldND.() -> R, ): R { contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) } - return nd(bufferFactory, *shape).run(action) + return ndAlgebra(bufferFactory, *shape).run(action) } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/DoubleFieldND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/DoubleFieldND.kt index fed4aca0b..0c7d4d5d1 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/DoubleFieldND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/DoubleFieldND.kt @@ -103,12 +103,12 @@ public class DoubleFieldND( override fun atanh(arg: StructureND): BufferND = arg.map { atanh(it) } } -public fun DoubleField.nd(vararg shape: Int): DoubleFieldND = DoubleFieldND(shape) +public fun DoubleField.ndAlgebra(vararg shape: Int): DoubleFieldND = DoubleFieldND(shape) /** * Produce a context for n-dimensional operations inside this real field */ -public inline fun DoubleField.withNd(vararg shape: Int, action: DoubleFieldND.() -> R): R { +public inline fun DoubleField.withNdAlgebra(vararg shape: Int, action: DoubleFieldND.() -> R): R { contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) } return DoubleFieldND(shape).run(action) } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/ShortRingND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/ShortRingND.kt index 5b6e2cd22..b56bef230 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/ShortRingND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/ShortRingND.kt @@ -34,7 +34,7 @@ public class ShortRingND( public inline fun BufferedRingND.produceInline(crossinline initializer: ShortRing.(Int) -> Short): BufferND = BufferND(strides, ShortBuffer(ShortArray(strides.linearSize) { offset -> ShortRing.initializer(offset) })) -public inline fun ShortRing.withNd(vararg shape: Int, action: ShortRingND.() -> R): R { +public inline fun ShortRing.withNdAlgebra(vararg shape: Int, action: ShortRingND.() -> R): R { contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) } return ShortRingND(shape).run(action) } 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 d0e2354d2..91cd5abe9 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 @@ -67,12 +67,14 @@ private class MutableStructure1DWrapper(val structure: MutableStructureND) structure[intArrayOf(index)] = value } - @PerformancePitfall + @OptIn(PerformancePitfall::class) override fun copy(): MutableBuffer = structure .elements() .map(Pair::second) .toMutableList() .asMutableBuffer() + + override fun toString(): String = Buffer.toString(this) } @@ -107,6 +109,8 @@ internal class MutableBuffer1DWrapper(val buffer: MutableBuffer) : Mutable } override fun copy(): MutableBuffer = buffer.copy() + + override fun toString(): String = Buffer.toString(this) } /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt index 7fb8ea251..57836a9ef 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt @@ -6,7 +6,6 @@ package space.kscience.kmath.nd import space.kscience.kmath.misc.PerformancePitfall -import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.MutableListBuffer import space.kscience.kmath.structures.VirtualBuffer @@ -108,7 +107,6 @@ private value class Structure2DWrapper(val structure: StructureND) : S override operator fun get(i: Int, j: Int): T = structure[i, j] - @UnstableKMathAPI override fun getFeature(type: KClass): F? = structure.getFeature(type) @PerformancePitfall diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt index 0675f5baf..6123336ba 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt @@ -9,12 +9,12 @@ import space.kscience.kmath.linear.LinearSpace import space.kscience.kmath.misc.Feature import space.kscience.kmath.misc.Featured import space.kscience.kmath.misc.PerformancePitfall -import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.operations.Ring import space.kscience.kmath.operations.invoke import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.BufferFactory import kotlin.jvm.JvmName +import kotlin.math.abs import kotlin.native.concurrent.ThreadLocal import kotlin.reflect.KClass @@ -61,7 +61,6 @@ public interface StructureND : Featured { * Feature is some additional structure information that allows to access it special properties or hints. * If the feature is not present, `null` is returned. */ - @UnstableKMathAPI override fun getFeature(type: KClass): F? = null public companion object { @@ -80,6 +79,22 @@ public interface StructureND : Featured { return st1.elements().all { (index, value) -> value == st2[index] } } + @PerformancePitfall + public fun contentEquals( + st1: StructureND, + st2: StructureND, + tolerance: Double = 1e-11 + ): Boolean { + if (st1 === st2) return true + + // fast comparison of buffers if possible + if (st1 is BufferND && st2 is BufferND && st1.strides == st2.strides) + return Buffer.contentEquals(st1.buffer, st2.buffer) + + //element by element comparison if it could not be avoided + return st1.elements().all { (index, value) -> abs(value - st2[index]) < tolerance } + } + /** * Debug output to string */ @@ -196,8 +211,8 @@ public fun > LinearSpace>.contentEquals( */ public operator fun StructureND.get(vararg index: Int): T = get(index) -@UnstableKMathAPI -public inline fun StructureND<*>.getFeature(): T? = getFeature(T::class) +//@UnstableKMathAPI +//public inline fun StructureND<*>.getFeature(): T? = getFeature(T::class) /** * Represents mutable [StructureND]. diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt index 9f4d741fd..e82b62c1b 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt @@ -17,6 +17,12 @@ import space.kscience.kmath.structures.DoubleBuffer public interface BufferAlgebra> : Algebra> { public val bufferFactory: BufferFactory public val elementAlgebra: A + public val size: Int + + public fun buffer(vararg elements: T): Buffer { + require(elements.size == size) { "Expected $size elements but found ${elements.size}" } + return bufferFactory(size) { elements[it] } + } //TODO move to multi-receiver inline extension public fun Buffer.map(block: (T) -> T): Buffer = bufferFactory(size) { block(get(it)) } @@ -39,6 +45,11 @@ public interface BufferAlgebra> : Algebra> { } } +@UnstableKMathAPI +public fun BufferField.buffer(initializer: (Int) -> T): Buffer { + return bufferFactory(size, initializer) +} + @UnstableKMathAPI public fun > BufferAlgebra.sin(arg: Buffer): Buffer = arg.map(elementAlgebra::sin) @@ -104,14 +115,9 @@ public fun > BufferAlgebra.pow(arg: Buffer, p public class BufferField>( override val bufferFactory: BufferFactory, override val elementAlgebra: A, - public val size: Int + override val size: Int ) : BufferAlgebra, Field> { - public fun produce(vararg elements: T): Buffer { - require(elements.size == size) { "Expected $size elements but found ${elements.size}" } - return bufferFactory(size) { elements[it] } - } - override val zero: Buffer = bufferFactory(size) { elementAlgebra.zero } override val one: Buffer = bufferFactory(size) { elementAlgebra.one } @@ -135,11 +141,15 @@ public class BufferField>( //Double buffer specialization @UnstableKMathAPI -public fun BufferField.produce(vararg elements: Number): Buffer { +public fun BufferField.buffer(vararg elements: Number): Buffer { require(elements.size == size) { "Expected $size elements but found ${elements.size}" } return bufferFactory(size) { elements[it].toDouble() } } +@UnstableKMathAPI +public fun > A.bufferAlgebra(bufferFactory: BufferFactory, size: Int): BufferField = + BufferField(bufferFactory, this, size) + @UnstableKMathAPI public fun DoubleField.bufferAlgebra(size: Int): BufferField = BufferField(::DoubleBuffer, DoubleField, size) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/ArrayBuffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/ArrayBuffer.kt index d49f70355..393ee99d6 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/ArrayBuffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/ArrayBuffer.kt @@ -23,6 +23,8 @@ public class ArrayBuffer(internal val array: Array) : MutableBuffer { override operator fun iterator(): Iterator = array.iterator() override fun copy(): MutableBuffer = ArrayBuffer(array.copyOf()) + + override fun toString(): String = Buffer.toString(this) } 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 c22a4ba39..fc23169ca 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 @@ -45,7 +45,12 @@ public interface Buffer { */ public operator fun iterator(): Iterator + override fun toString(): String + public companion object { + + public fun toString(buffer: Buffer<*>): String = buffer.asSequence().joinToString(prefix = "[", separator = ", ", postfix = "]") + /** * Check the element-by-element match of content of two buffers. */ @@ -126,6 +131,8 @@ public class VirtualBuffer(override val size: Int, private val generator: } override operator fun iterator(): Iterator = (0 until size).asSequence().map(generator).iterator() + + override fun toString(): String = Buffer.toString(this) } /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt index 005a693eb..d6a48f42d 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt @@ -49,6 +49,8 @@ internal class BufferAccessor2D( override fun copy(): MutableBuffer = factory(colNum) { get(it) } override operator fun iterator(): Iterator = (0 until colNum).map(::get).iterator() + override fun toString(): String = Buffer.toString(this) + } /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/FlaggedBuffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/FlaggedBuffer.kt index 665558829..700a4f17f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/FlaggedBuffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/FlaggedBuffer.kt @@ -53,8 +53,10 @@ public fun FlaggedBuffer<*>.isMissing(index: Int): Boolean = hasFlag(index, Valu /** * A [Double] buffer that supports flags for each value like `NaN` or Missing. */ -public class FlaggedDoubleBuffer(public val values: DoubleArray, public val flags: ByteArray) : FlaggedBuffer, - Buffer { +public class FlaggedDoubleBuffer( + public val values: DoubleArray, + public val flags: ByteArray +) : FlaggedBuffer, Buffer { init { require(values.size == flags.size) { "Values and flags must have the same dimensions" } } @@ -68,6 +70,8 @@ public class FlaggedDoubleBuffer(public val values: DoubleArray, public val flag override operator fun iterator(): Iterator = values.indices.asSequence().map { if (isValid(it)) values[it] else null }.iterator() + + override fun toString(): String = Buffer.toString(this) } public inline fun FlaggedDoubleBuffer.forEachValid(block: (Double) -> Unit) { diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/ListBuffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/ListBuffer.kt index fdba68d19..666722177 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/ListBuffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/ListBuffer.kt @@ -21,6 +21,8 @@ public class ListBuffer(public val list: List) : Buffer { override operator fun get(index: Int): T = list[index] override operator fun iterator(): Iterator = list.iterator() + + override fun toString(): String = Buffer.toString(this) } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MemoryBuffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MemoryBuffer.kt index 996785570..3e08dbbb1 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MemoryBuffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MemoryBuffer.kt @@ -22,6 +22,8 @@ public open class MemoryBuffer(protected val memory: Memory, protected override operator fun get(index: Int): T = reader.read(spec, spec.objectSize * index) override operator fun iterator(): Iterator = (0 until size).asSequence().map { get(it) }.iterator() + override fun toString(): String = Buffer.toString(this) + public companion object { public fun create(spec: MemorySpec, size: Int): MemoryBuffer = MemoryBuffer(Memory.allocate(size * spec.objectSize), spec) diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NDFieldTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NDFieldTest.kt index 6d4531f41..05a67ab09 100644 --- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NDFieldTest.kt +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NDFieldTest.kt @@ -6,7 +6,7 @@ package space.kscience.kmath.structures import space.kscience.kmath.nd.get -import space.kscience.kmath.nd.nd +import space.kscience.kmath.nd.ndAlgebra import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.invoke import space.kscience.kmath.testutils.FieldVerifier @@ -16,12 +16,12 @@ import kotlin.test.assertEquals internal class NDFieldTest { @Test fun verify() { - (DoubleField.nd(12, 32)) { FieldVerifier(this, one + 3, one - 23, one * 12, 6.66) } + (DoubleField.ndAlgebra(12, 32)) { FieldVerifier(this, one + 3, one - 23, one * 12, 6.66) } } @Test fun testStrides() { - val ndArray = DoubleField.nd(10, 10).produce { (it[0] + it[1]).toDouble() } + val ndArray = DoubleField.ndAlgebra(10, 10).produce { (it[0] + it[1]).toDouble() } assertEquals(ndArray[5, 5], 10.0) } } diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt index 1045933b7..b2982b335 100644 --- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt @@ -10,7 +10,7 @@ import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.nd.StructureND import space.kscience.kmath.nd.combine import space.kscience.kmath.nd.get -import space.kscience.kmath.nd.nd +import space.kscience.kmath.nd.ndAlgebra import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.Norm import space.kscience.kmath.operations.invoke @@ -21,7 +21,7 @@ import kotlin.test.assertEquals @Suppress("UNUSED_VARIABLE") class NumberNDFieldTest { - val algebra = DoubleField.nd(3, 3) + val algebra = DoubleField.ndAlgebra(3, 3) val array1 = algebra.produce { (i, j) -> (i + j).toDouble() } val array2 = algebra.produce { (i, j) -> (i - j).toDouble() } @@ -87,7 +87,7 @@ class NumberNDFieldTest { @Test fun testInternalContext() { algebra { - (DoubleField.nd(*array1.shape)) { with(L2Norm) { 1 + norm(array1) + exp(array2) } } + (DoubleField.ndAlgebra(*array1.shape)) { with(L2Norm) { 1 + norm(array1) + exp(array2) } } } } } diff --git a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/RingBuffer.kt b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/RingBuffer.kt index 4b12b031d..573b406e2 100644 --- a/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/RingBuffer.kt +++ b/kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/RingBuffer.kt @@ -69,6 +69,8 @@ public class RingBuffer( @Suppress("NOTHING_TO_INLINE") private inline fun Int.forward(n: Int): Int = (this + n) % (buffer.size) + override fun toString(): String = Buffer.toString(this) + public companion object { public inline fun build(size: Int, empty: T): RingBuffer { val buffer = MutableBuffer.auto(size) { empty } as MutableBuffer diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt index 022a7874e..4a53af60d 100644 --- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt +++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt @@ -5,9 +5,12 @@ package space.kscience.kmath.ejml +import space.kscience.kmath.linear.InverseMatrixFeature import space.kscience.kmath.linear.LinearSpace import space.kscience.kmath.linear.Matrix import space.kscience.kmath.linear.Point +import space.kscience.kmath.misc.UnstableKMathAPI +import space.kscience.kmath.nd.Structure2D import space.kscience.kmath.operations.Ring /** @@ -36,4 +39,7 @@ public abstract class EjmlLinearSpace, out M : org.ejml ): EjmlMatrix public abstract override fun buildVector(size: Int, initializer: A.(Int) -> T): EjmlVector + + @UnstableKMathAPI + public fun EjmlMatrix.inverse(): Structure2D = computeFeature(this, InverseMatrixFeature::class)?.inverse as Structure2D } diff --git a/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt b/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt index 5667815ac..5b8b2af98 100644 --- a/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt +++ b/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt @@ -9,12 +9,11 @@ import org.ejml.data.DMatrixRMaj import org.ejml.dense.row.CommonOps_DDRM import org.ejml.dense.row.RandomMatrices_DDRM import org.ejml.dense.row.factory.DecompositionFactory_DDRM -import space.kscience.kmath.linear.DeterminantFeature -import space.kscience.kmath.linear.LupDecompositionFeature -import space.kscience.kmath.linear.computeFeature +import space.kscience.kmath.linear.* import space.kscience.kmath.misc.PerformancePitfall import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.nd.StructureND +import space.kscience.kmath.operations.algebra import kotlin.random.Random import kotlin.random.asJavaRandom import kotlin.test.* @@ -82,4 +81,24 @@ internal class EjmlMatrixTest { val m = randomMatrix assertSame(m, EjmlDoubleMatrix(m).origin) } + + @Test + fun inverse() = EjmlLinearSpaceDDRM { + val random = Random(1224) + val dim = 20 + + val space = Double.algebra.linearSpace + + //creating invertible matrix + val u = space.buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } + val l = space.buildMatrix(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 } + val matrix = space { l dot u } + val inverted = matrix.toEjml().inverse() + + val res = matrix dot inverted + + println(StructureND.toString(res)) + + assertTrue { StructureND.contentEquals(one(dim, dim), res, 1e-3) } + } } 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 098beb2cf..88932d59b 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 @@ -12,6 +12,7 @@ import space.kscience.kmath.linear.* 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.structures.Buffer import space.kscience.kmath.structures.DoubleBuffer import space.kscience.kmath.structures.asIterable @@ -32,18 +33,18 @@ import kotlin.math.pow public typealias RealMatrix = Matrix public fun realMatrix(rowNum: Int, colNum: Int, initializer: DoubleField.(i: Int, j: Int) -> Double): RealMatrix = - LinearSpace.double.buildMatrix(rowNum, colNum, initializer) + Double.algebra.linearSpace.buildMatrix(rowNum, colNum, initializer) @OptIn(UnstableKMathAPI::class) public fun realMatrix(rowNum: Int, colNum: Int): MatrixBuilder = - LinearSpace.double.matrix(rowNum, colNum) + Double.algebra.linearSpace.matrix(rowNum, colNum) public fun Array.toMatrix(): RealMatrix { - return LinearSpace.double.buildMatrix(size, this[0].size) { row, col -> this@toMatrix[row][col] } + return Double.algebra.linearSpace.buildMatrix(size, this[0].size) { row, col -> this@toMatrix[row][col] } } public fun Sequence.toMatrix(): RealMatrix = toList().let { - LinearSpace.double.buildMatrix(it.size, it[0].size) { row, col -> it[row][col] } + Double.algebra.linearSpace.buildMatrix(it.size, it[0].size) { row, col -> it[row][col] } } public fun RealMatrix.repeatStackVertical(n: Int): RealMatrix = @@ -56,37 +57,37 @@ public fun RealMatrix.repeatStackVertical(n: Int): RealMatrix = */ public operator fun RealMatrix.times(double: Double): RealMatrix = - LinearSpace.double.buildMatrix(rowNum, colNum) { row, col -> + Double.algebra.linearSpace.buildMatrix(rowNum, colNum) { row, col -> get(row, col) * double } public operator fun RealMatrix.plus(double: Double): RealMatrix = - LinearSpace.double.buildMatrix(rowNum, colNum) { row, col -> + Double.algebra.linearSpace.buildMatrix(rowNum, colNum) { row, col -> get(row, col) + double } public operator fun RealMatrix.minus(double: Double): RealMatrix = - LinearSpace.double.buildMatrix(rowNum, colNum) { row, col -> + Double.algebra.linearSpace.buildMatrix(rowNum, colNum) { row, col -> get(row, col) - double } public operator fun RealMatrix.div(double: Double): RealMatrix = - LinearSpace.double.buildMatrix(rowNum, colNum) { row, col -> + Double.algebra.linearSpace.buildMatrix(rowNum, colNum) { row, col -> get(row, col) / double } public operator fun Double.times(matrix: RealMatrix): RealMatrix = - LinearSpace.double.buildMatrix(matrix.rowNum, matrix.colNum) { row, col -> + Double.algebra.linearSpace.buildMatrix(matrix.rowNum, matrix.colNum) { row, col -> this@times * matrix[row, col] } public operator fun Double.plus(matrix: RealMatrix): RealMatrix = - LinearSpace.double.buildMatrix(matrix.rowNum, matrix.colNum) { row, col -> + Double.algebra.linearSpace.buildMatrix(matrix.rowNum, matrix.colNum) { row, col -> this@plus + matrix[row, col] } public operator fun Double.minus(matrix: RealMatrix): RealMatrix = - LinearSpace.double.buildMatrix(matrix.rowNum, matrix.colNum) { row, col -> + Double.algebra.linearSpace.buildMatrix(matrix.rowNum, matrix.colNum) { row, col -> this@minus - matrix[row, col] } @@ -101,20 +102,20 @@ public operator fun Double.minus(matrix: RealMatrix): RealMatrix = @UnstableKMathAPI public operator fun RealMatrix.times(other: RealMatrix): RealMatrix = - LinearSpace.double.buildMatrix(rowNum, colNum) { row, col -> this@times[row, col] * other[row, col] } + Double.algebra.linearSpace.buildMatrix(rowNum, colNum) { row, col -> this@times[row, col] * other[row, col] } public operator fun RealMatrix.plus(other: RealMatrix): RealMatrix = - LinearSpace.double.run { this@plus + other } + Double.algebra.linearSpace.run { this@plus + other } public operator fun RealMatrix.minus(other: RealMatrix): RealMatrix = - LinearSpace.double.buildMatrix(rowNum, colNum) { row, col -> this@minus[row, col] - other[row, col] } + Double.algebra.linearSpace.buildMatrix(rowNum, colNum) { row, col -> this@minus[row, col] - other[row, col] } /* * Operations on columns */ public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer) -> Double): RealMatrix = - LinearSpace.double.buildMatrix(rowNum, colNum + 1) { row, col -> + Double.algebra.linearSpace.buildMatrix(rowNum, colNum + 1) { row, col -> if (col < colNum) get(row, col) else @@ -122,7 +123,7 @@ public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer) - } public fun RealMatrix.extractColumns(columnRange: IntRange): RealMatrix = - LinearSpace.double.buildMatrix(rowNum, columnRange.count()) { row, col -> + Double.algebra.linearSpace.buildMatrix(rowNum, columnRange.count()) { row, col -> this@extractColumns[row, columnRange.first + col] } @@ -155,14 +156,14 @@ public fun RealMatrix.max(): Double? = elements().map { (_, value) -> value }.ma public fun RealMatrix.average(): Double = elements().map { (_, value) -> value }.average() public inline fun RealMatrix.map(crossinline transform: (Double) -> Double): RealMatrix = - LinearSpace.double.buildMatrix(rowNum, colNum) { i, j -> + Double.algebra.linearSpace.buildMatrix(rowNum, colNum) { i, j -> transform(get(i, j)) } /** * Inverse a square real matrix using LUP decomposition */ -public fun RealMatrix.inverseWithLup(): RealMatrix = LinearSpace.double.lupSolver().inverse(this) +public fun RealMatrix.inverseWithLup(): RealMatrix = Double.algebra.linearSpace.lupSolver().inverse(this) //extended operations diff --git a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/dot.kt b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/dot.kt index ca2db8131..883a63f46 100644 --- a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/dot.kt +++ b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/dot.kt @@ -5,13 +5,14 @@ package space.kscience.kmath.real -import space.kscience.kmath.linear.LinearSpace import space.kscience.kmath.linear.Matrix +import space.kscience.kmath.linear.linearSpace +import space.kscience.kmath.operations.algebra /** * Optimized dot product for real matrices */ -public infix fun Matrix.dot(other: Matrix): Matrix = LinearSpace.double.run { +public infix fun Matrix.dot(other: Matrix): Matrix = Double.algebra.linearSpace.run { this@dot dot other } \ No newline at end of file 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 c99f3ff99..662cdf3d7 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 @@ -10,6 +10,7 @@ import space.kscience.kmath.functions.integrate import space.kscience.kmath.interpolation.PolynomialInterpolator 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 @@ -23,6 +24,7 @@ import space.kscience.kmath.structures.map /** * Compute analytical indefinite integral of this [PiecewisePolynomial], keeping all intervals intact */ +@OptIn(PerformancePitfall::class) @UnstableKMathAPI public fun > PiecewisePolynomial.integrate(algebra: Field): PiecewisePolynomial = PiecewisePolynomial(pieces.map { it.first to it.second.integrate(algebra) }) diff --git a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/DoubleHistogramSpace.kt b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/DoubleHistogramSpace.kt index 37ab8a1b2..5446f05f8 100644 --- a/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/DoubleHistogramSpace.kt +++ b/kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/DoubleHistogramSpace.kt @@ -29,7 +29,7 @@ public class DoubleHistogramSpace( public val dimension: Int get() = lower.size private val shape = IntArray(binNums.size) { binNums[it] + 2 } - override val histogramValueSpace: DoubleFieldND = DoubleField.nd(*shape) + override val histogramValueSpace: DoubleFieldND = DoubleField.ndAlgebra(*shape) override val strides: Strides get() = histogramValueSpace.strides private val binSize = DoubleBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] } diff --git a/kmath-viktor/src/main/kotlin/space/kscience/kmath/viktor/ViktorBuffer.kt b/kmath-viktor/src/main/kotlin/space/kscience/kmath/viktor/ViktorBuffer.kt index caebd9783..32fb65b8a 100644 --- a/kmath-viktor/src/main/kotlin/space/kscience/kmath/viktor/ViktorBuffer.kt +++ b/kmath-viktor/src/main/kotlin/space/kscience/kmath/viktor/ViktorBuffer.kt @@ -6,10 +6,12 @@ package space.kscience.kmath.viktor import org.jetbrains.bio.viktor.F64FlatArray +import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.MutableBuffer @Suppress("NOTHING_TO_INLINE", "OVERRIDE_BY_INLINE") -public class ViktorBuffer(public val flatArray: F64FlatArray) : MutableBuffer { +@JvmInline +public value class ViktorBuffer(public val flatArray: F64FlatArray) : MutableBuffer { override val size: Int get() = flatArray.size @@ -21,4 +23,6 @@ public class ViktorBuffer(public val flatArray: F64FlatArray) : MutableBuffer = ViktorBuffer(flatArray.copy().flatten()) override operator fun iterator(): Iterator = flatArray.data.iterator() + + override fun toString(): String = Buffer.toString(this) }