Fix EJML inversion issue

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: `<Type>.algebra.<nd/etc>`
+- Adjusted commons-math API for linear solvers to match conventions.
 
 ### Deprecated
 - Specialized `DoubleBufferAlgebra`

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

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) {
+    fun doubleDot(blackhole: Blackhole) {
-        LinearSpace.double {
+        with(Double.algebra.linearSpace) {
             blackhole.consume(matrix1 dot matrix2)
         }
     }

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<InverseMatrixFeature<Double>>()?.inverse)
+            blackhole.consume(matrix.toEjml().inverse())
         }
     }
 }

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 autoField = DoubleField.autoNdAlgebra(dim, dim)
-        private val specializedField = DoubleField.nd(dim, dim)
+        private val specializedField = DoubleField.ndAlgebra(dim, dim)
-        private val genericField = DoubleField.nd(Buffer.Companion::boxing, dim, dim)
+        private val genericField = DoubleField.ndAlgebra(Buffer.Companion::boxing, dim, dim)
     }
 }

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.autoNdAlgebra
-import space.kscience.kmath.nd.nd
+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 autoField = DoubleField.autoNdAlgebra(dim, dim)
-        private val realField = DoubleField.nd(dim, dim)
+        private val realField = DoubleField.ndAlgebra(dim, dim)
         private val viktorField = ViktorNDField(dim, dim)
     }
 }

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.autoNdAlgebra
-import space.kscience.kmath.nd.nd
+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 autoField = DoubleField.autoNdAlgebra(dim, dim)
-        private val realNdField = DoubleField.nd(dim, dim)
+        private val realNdField = DoubleField.ndAlgebra(dim, dim)
         private val viktorField = ViktorFieldND(intArrayOf(dim, dim))
     }
 }

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.nd.withNdAlgebra
-import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.algebra
 import space.kscience.kmath.operations.invoke
 
-fun main(): Unit = DoubleField {
+fun main(): Unit = Double.algebra {
-    withNd(2, 2) {
+    withNdAlgebra(2, 2) {
 
         //Produce a diagonal StructureND
         fun diagonal(v: Double) = produce { (i, j) ->

examples/src/main/kotlin/space/kscience/kmath/operations/ComplexDemo.kt

@@ -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<Complex> = 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)
-}

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<Complex> = ndAlgebra(2, 2).produce { (i, j) ->
+        Complex(i - j, i + j)
+    }
+    println(element)
+
+    // 1d element operation
+    val result: StructureND<Complex> = 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)
+}

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 realField = DoubleField.ndAlgebra(dim, dim)
-    val complexField: ComplexFieldND = ComplexField.nd(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)

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.autoNdAlgebra
-import space.kscience.kmath.nd.nd
+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

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 <reified R : Any> MutableBuffer.Companion.same(
     n: Int,
@@ -17,6 +17,6 @@ inline fun <reified R : Any> 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))
     }
 }

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<Double> {
     override operator fun get(i: Int, j: Int): Double = origin.getEntry(i, j)
 }
 
-public class CMVector(public val origin: RealVector) : Point<Double> {
+@JvmInline
+public value class CMVector(public val origin: RealVector) : Point<Double> {
     override val size: Int get() = origin.dimension
 
     override operator fun get(index: Int): Double = origin.getEntry(index)
 
     override operator fun iterator(): Iterator<Double> = origin.toArray().iterator()
+
+    override fun toString(): String = Buffer.toString(this)
 }
 
 public fun RealVector.toPoint(): CMVector = CMVector(this)

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<Double>,
     decomposition: CMDecomposition = CMDecomposition.LUP,
 ): DecompositionSolver = when (decomposition) {
@@ -48,9 +48,11 @@ public fun CMLinearSpace.inverse(
 
 
 public fun CMLinearSpace.solver(decomposition: CMDecomposition): LinearSolver<Double> = object : LinearSolver<Double> {
-    override fun solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> = solve(a, b, decomposition)
+    override fun solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> = solver(a, decomposition).solve(b.toCM().origin).wrap()
 
-    override fun solve(a: Matrix<Double>, b: Point<Double>): Point<Double> = solve(a, b, decomposition)
+    override fun solve(a: Matrix<Double>, b: Point<Double>): Point<Double> = solver(a, decomposition).solve(b.toCM().origin).toPoint()
 
-    override fun inverse(matrix: Matrix<Double>): Matrix<Double> = inverse(matrix, decomposition)
+    override fun inverse(matrix: Matrix<Double>): Matrix<Double> = solver(matrix, decomposition).inverse.wrap()
 }
+
+public fun CMLinearSpace.lupSolver(): LinearSolver<Double> = solver((CMDecomposition.LUP))

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<Complex, ComplexField>.produceInline(initializ
     return BufferND(strides, buffer)
 }
 
+@UnstableKMathAPI
+public fun ComplexField.bufferAlgebra(size: Int): BufferField<Complex, ComplexField> =
+    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 <R> ComplexField.withNd(vararg shape: Int, action: ComplexFieldND.() -> R): R {
+public inline fun <R> ComplexField.withNdAlgebra(vararg shape: Int, action: ComplexFieldND.() -> R): R {
     contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
     return ComplexFieldND(shape).action()
 }

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.*
-import space.kscience.kmath.structures.BufferFactory
-import space.kscience.kmath.structures.VirtualBuffer
-import space.kscience.kmath.structures.indices
 
 
-public class BufferedLinearSpace<T : Any, out A : Ring<T>>(
+public class BufferedLinearSpace<T, out A : Ring<T>>(
     override val elementAlgebra: A,
     private val bufferFactory: BufferFactory<T>,
 ) : LinearSpace<T, A> {
@@ -26,7 +24,7 @@ public class BufferedLinearSpace<T : Any, out A : Ring<T>>(
     private fun ndRing(
         rows: Int,
         cols: Int,
-    ): BufferedRingND<T, A> = elementAlgebra.nd(bufferFactory, rows, cols)
+    ): BufferedRingND<T, A> = elementAlgebra.ndAlgebra(bufferFactory, rows, cols)
 
     override fun buildMatrix(rows: Int, columns: Int, initializer: A.(i: Int, j: Int) -> T): Matrix<T> =
         ndRing(rows, columns).produce { (i, j) -> elementAlgebra.initializer(i, j) }.as2D()
@@ -92,3 +90,10 @@ public class BufferedLinearSpace<T : Any, out A : Ring<T>>(
         unwrap().map { it * value }.as2D()
     }
 }
+
+
+public fun <T, A : Ring<T>> A.linearSpace(bufferFactory: BufferFactory<T>): BufferedLinearSpace<T, A> =
+    BufferedLinearSpace(this, bufferFactory)
+
+public val DoubleField.linearSpace: BufferedLinearSpace<Double, DoubleField>
+    get() = BufferedLinearSpace(this, ::DoubleBuffer)

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.nd.MutableStructure2D
-import space.kscience.kmath.operations.*
+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<T> = Buffer<T>
  * @param T the type of items in the matrices.
  * @param A the type of ring over [T].
  */
-public interface LinearSpace<T : Any, out A : Ring<T>> {
+public interface LinearSpace<T, out A : Ring<T>> {
     public val elementAlgebra: A
 
     /**
@@ -172,7 +177,8 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
      * @return a feature object or `null` if it isn't present.
      */
     @UnstableKMathAPI
-    public fun <F : StructureFeature> computeFeature(structure: Matrix<T>, type: KClass<out F>): F? = structure.getFeature(type)
+    public fun <F : StructureFeature> computeFeature(structure: Matrix<T>, type: KClass<out F>): F? =
+        structure.getFeature(type)
 
     public companion object {
 
@@ -184,6 +190,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
             bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
         ): LinearSpace<T, A> = BufferedLinearSpace(algebra, bufferFactory)
 
+        @Deprecated("use DoubleField.linearSpace")
         public val double: LinearSpace<Double, DoubleField> = buffered(DoubleField, ::DoubleBuffer)
 
         /**
@@ -213,10 +220,8 @@ public inline operator fun <LS : LinearSpace<*, *>, R> LS.invoke(block: LS.() ->
  * Convert matrix to vector if it is possible.
  */
 public fun <T : Any> Matrix<T>.asVector(): Point<T> =
-    if (this.colNum == 1)
+    if (this.colNum == 1) as1D()
-        as1D()
+    else error("Can't convert matrix with more than one column to vector")
-    else
-        error("Can't convert matrix with more than one column to vector")
 
 /**
  * Creates an n &times; 1 [VirtualMatrix], where n is the size of the given buffer.

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 <T : Any> LupDecomposition<T>.solve(
 /**
  * Produce a generic solver based on LUP decomposition
  */
-@PerformancePitfall()
 @OptIn(UnstableKMathAPI::class)
 public fun <T : Comparable<T>, F : Field<T>> LinearSpace<T, F>.lupSolver(
     bufferFactory: MutableBufferFactory<T>,
@@ -222,13 +219,12 @@ public fun <T : Comparable<T>, F : Field<T>> LinearSpace<T, F>.lupSolver(
 ): LinearSolver<T> = object : LinearSolver<T> {
     override fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
         // 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<T>): Matrix<T> = solve(matrix, one(matrix.rowNum, matrix.colNum))
 }
 
-@PerformancePitfall
 public fun LinearSpace<Double, DoubleField>.lupSolver(singularityThreshold: Double = 1e-11): LinearSolver<Double> =
     lupSolver(::DoubleBuffer) { it < singularityThreshold }

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<out T : Any> 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 <F : StructureFeature> getFeature(type: KClass<out F>): F? =
         features.getFeature(type) ?: origin.getFeature(type)
@@ -90,8 +88,7 @@ public class TransposedFeature<out T : Any>(public val original: Matrix<T>) : Ma
 /**
  * Create a virtual transposed matrix without copying anything. `A.transpose().transpose() === A`
  */
+@Suppress("UNCHECKED_CAST")
 @OptIn(UnstableKMathAPI::class)
-public fun <T : Any> Matrix<T>.transpose(): Matrix<T> = getFeature<TransposedFeature<T>>()?.original ?: VirtualMatrix(
+public fun <T : Any> Matrix<T>.transpose(): Matrix<T> = getFeature(TransposedFeature::class)?.original as? Matrix<T>
-    colNum,
+    ?: VirtualMatrix(colNum, rowNum) { i, j -> get(j, i) }.withFeature(TransposedFeature(this))
-    rowNum,
-) { i, j -> get(j, i) }.withFeature(TransposedFeature(this))

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt

@@ -87,39 +87,39 @@ public open class BufferedFieldND<T, out R : Field<T>>(
 }
 
 // group factories
-public fun <T, A : Group<T>> A.nd(
+public fun <T, A : Group<T>> A.ndAlgebra(
     bufferFactory: BufferFactory<T>,
     vararg shape: Int,
 ): BufferedGroupND<T, A> = BufferedGroupND(shape, this, bufferFactory)
 
 @JvmName("withNdGroup")
-public inline fun <T, A : Group<T>, R> A.withNd(
+public inline fun <T, A : Group<T>, R> A.withNdAlgebra(
     noinline bufferFactory: BufferFactory<T>,
     vararg shape: Int,
     action: BufferedGroupND<T, A>.() -> R,
 ): R {
     contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
-    return nd(bufferFactory, *shape).run(action)
+    return ndAlgebra(bufferFactory, *shape).run(action)
 }
 
 //ring factories
-public fun <T, A : Ring<T>> A.nd(
+public fun <T, A : Ring<T>> A.ndAlgebra(
     bufferFactory: BufferFactory<T>,
     vararg shape: Int,
 ): BufferedRingND<T, A> = BufferedRingND(shape, this, bufferFactory)
 
 @JvmName("withNdRing")
-public inline fun <T, A : Ring<T>, R> A.withNd(
+public inline fun <T, A : Ring<T>, R> A.withNdAlgebra(
     noinline bufferFactory: BufferFactory<T>,
     vararg shape: Int,
     action: BufferedRingND<T, A>.() -> R,
 ): R {
     contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
-    return nd(bufferFactory, *shape).run(action)
+    return ndAlgebra(bufferFactory, *shape).run(action)
 }
 
 //field factories
-public fun <T, A : Field<T>> A.nd(
+public fun <T, A : Field<T>> A.ndAlgebra(
     bufferFactory: BufferFactory<T>,
     vararg shape: Int,
 ): BufferedFieldND<T, A> = BufferedFieldND(shape, this, bufferFactory)
@@ -129,7 +129,7 @@ public fun <T, A : Field<T>> A.nd(
  */
 @UnstableKMathAPI
 @Suppress("UNCHECKED_CAST")
-public inline fun <reified T : Any, A : Field<T>> A.autoNd(
+public inline fun <reified T : Any, A : Field<T>> A.autoNdAlgebra(
     vararg shape: Int,
 ): FieldND<T, A> = when (this) {
     DoubleField -> DoubleFieldND(shape) as FieldND<T, A>
@@ -137,11 +137,11 @@ public inline fun <reified T : Any, A : Field<T>> A.autoNd(
 }
 
 @JvmName("withNdField")
-public inline fun <T, A : Field<T>, R> A.withNd(
+public inline fun <T, A : Field<T>, R> A.withNdAlgebra(
     noinline bufferFactory: BufferFactory<T>,
     vararg shape: Int,
     action: BufferedFieldND<T, A>.() -> R,
 ): R {
     contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
-    return nd(bufferFactory, *shape).run(action)
+    return ndAlgebra(bufferFactory, *shape).run(action)
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/DoubleFieldND.kt

@@ -103,12 +103,12 @@ public class DoubleFieldND(
     override fun atanh(arg: StructureND<Double>): BufferND<Double> = 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 <R> DoubleField.withNd(vararg shape: Int, action: DoubleFieldND.() -> R): R {
+public inline fun <R> DoubleField.withNdAlgebra(vararg shape: Int, action: DoubleFieldND.() -> R): R {
     contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
     return DoubleFieldND(shape).run(action)
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/ShortRingND.kt

@@ -34,7 +34,7 @@ public class ShortRingND(
 public inline fun BufferedRingND<Short, ShortRing>.produceInline(crossinline initializer: ShortRing.(Int) -> Short): BufferND<Short> =
     BufferND(strides, ShortBuffer(ShortArray(strides.linearSize) { offset -> ShortRing.initializer(offset) }))
 
-public inline fun <R> ShortRing.withNd(vararg shape: Int, action: ShortRingND.() -> R): R {
+public inline fun <R> ShortRing.withNdAlgebra(vararg shape: Int, action: ShortRingND.() -> R): R {
     contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
     return ShortRingND(shape).run(action)
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt

@@ -67,12 +67,14 @@ private class MutableStructure1DWrapper<T>(val structure: MutableStructureND<T>)
         structure[intArrayOf(index)] = value
     }
 
-    @PerformancePitfall
+    @OptIn(PerformancePitfall::class)
     override fun copy(): MutableBuffer<T> = structure
         .elements()
         .map(Pair<IntArray, T>::second)
         .toMutableList()
         .asMutableBuffer()
+
+    override fun toString(): String = Buffer.toString(this)
 }
 
 
@@ -107,6 +109,8 @@ internal class MutableBuffer1DWrapper<T>(val buffer: MutableBuffer<T>) : Mutable
     }
 
     override fun copy(): MutableBuffer<T> = buffer.copy()
+
+    override fun toString(): String = Buffer.toString(this)
 }
 
 /**

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<out T>(val structure: StructureND<T>) : S
 
     override operator fun get(i: Int, j: Int): T = structure[i, j]
 
-    @UnstableKMathAPI
     override fun <F : StructureFeature> getFeature(type: KClass<out F>): F? = structure.getFeature(type)
 
     @PerformancePitfall

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<out T> : Featured<StructureFeature> {
      * 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 <F : StructureFeature> getFeature(type: KClass<out F>): F? = null
 
     public companion object {
@@ -80,6 +79,22 @@ public interface StructureND<out T> : Featured<StructureFeature> {
             return st1.elements().all { (index, value) -> value == st2[index] }
         }
 
+        @PerformancePitfall
+        public fun contentEquals(
+            st1: StructureND<Double>,
+            st2: StructureND<Double>,
+            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 <T : Comparable<T>> LinearSpace<T, Ring<T>>.contentEquals(
  */
 public operator fun <T> StructureND<T>.get(vararg index: Int): T = get(index)
 
-@UnstableKMathAPI
+//@UnstableKMathAPI
-public inline fun <reified T : StructureFeature> StructureND<*>.getFeature(): T? = getFeature(T::class)
+//public inline fun <reified T : StructureFeature> StructureND<*>.getFeature(): T? = getFeature(T::class)
 
 /**
  * Represents mutable [StructureND].

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt

@@ -17,6 +17,12 @@ import space.kscience.kmath.structures.DoubleBuffer
 public interface BufferAlgebra<T, A : Algebra<T>> : Algebra<Buffer<T>> {
     public val bufferFactory: BufferFactory<T>
     public val elementAlgebra: A
+    public val size: Int
+
+    public fun buffer(vararg elements: T): Buffer<T> {
+        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<T>.map(block: (T) -> T): Buffer<T> = bufferFactory(size) { block(get(it)) }
@@ -39,6 +45,11 @@ public interface BufferAlgebra<T, A : Algebra<T>> : Algebra<Buffer<T>> {
     }
 }
 
+@UnstableKMathAPI
+public fun <T> BufferField<T, *>.buffer(initializer: (Int) -> T): Buffer<T> {
+    return bufferFactory(size, initializer)
+}
+
 @UnstableKMathAPI
 public fun <T, A : TrigonometricOperations<T>> BufferAlgebra<T, A>.sin(arg: Buffer<T>): Buffer<T> =
     arg.map(elementAlgebra::sin)
@@ -104,14 +115,9 @@ public fun <T, A : PowerOperations<T>> BufferAlgebra<T, A>.pow(arg: Buffer<T>, p
 public class BufferField<T, A : Field<T>>(
     override val bufferFactory: BufferFactory<T>,
     override val elementAlgebra: A,
-    public val size: Int
+    override val size: Int
 ) : BufferAlgebra<T, A>, Field<Buffer<T>> {
 
-    public fun produce(vararg elements: T): Buffer<T> {
-        require(elements.size == size) { "Expected $size elements but found ${elements.size}" }
-        return bufferFactory(size) { elements[it] }
-    }
-
     override val zero: Buffer<T> = bufferFactory(size) { elementAlgebra.zero }
     override val one: Buffer<T> = bufferFactory(size) { elementAlgebra.one }
 
@@ -135,11 +141,15 @@ public class BufferField<T, A : Field<T>>(
 //Double buffer specialization
 
 @UnstableKMathAPI
-public fun BufferField<Double, *>.produce(vararg elements: Number): Buffer<Double> {
+public fun BufferField<Double, *>.buffer(vararg elements: Number): Buffer<Double> {
     require(elements.size == size) { "Expected $size elements but found ${elements.size}" }
     return bufferFactory(size) { elements[it].toDouble() }
 }
 
+@UnstableKMathAPI
+public fun <T, A : Field<T>> A.bufferAlgebra(bufferFactory: BufferFactory<T>, size: Int): BufferField<T, A> =
+    BufferField(bufferFactory, this, size)
+
 @UnstableKMathAPI
 public fun DoubleField.bufferAlgebra(size: Int): BufferField<Double, DoubleField> =
     BufferField(::DoubleBuffer, DoubleField, size)

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/ArrayBuffer.kt

@@ -23,6 +23,8 @@ public class ArrayBuffer<T>(internal val array: Array<T>) : MutableBuffer<T> {
 
     override operator fun iterator(): Iterator<T> = array.iterator()
     override fun copy(): MutableBuffer<T> = ArrayBuffer(array.copyOf())
+
+    override fun toString(): String = Buffer.toString(this)
 }
 
 

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt

@@ -45,7 +45,12 @@ public interface Buffer<out T> {
      */
     public operator fun iterator(): Iterator<T>
 
+    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<out T>(override val size: Int, private val generator:
     }
 
     override operator fun iterator(): Iterator<T> = (0 until size).asSequence().map(generator).iterator()
+
+    override fun toString(): String = Buffer.toString(this)
 }
 
 /**

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt

@@ -49,6 +49,8 @@ internal class BufferAccessor2D<T>(
         override fun copy(): MutableBuffer<T> = factory(colNum) { get(it) }
         override operator fun iterator(): Iterator<T> = (0 until colNum).map(::get).iterator()
 
+        override fun toString(): String = Buffer.toString(this)
+
     }
 
     /**

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<Double?>,
+public class FlaggedDoubleBuffer(
-    Buffer<Double?> {
+    public val values: DoubleArray,
+    public val flags: ByteArray
+) : FlaggedBuffer<Double?>, Buffer<Double?> {
     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<Double?> = 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) {

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/ListBuffer.kt

@@ -21,6 +21,8 @@ public class ListBuffer<T>(public val list: List<T>) : Buffer<T> {
 
     override operator fun get(index: Int): T = list[index]
     override operator fun iterator(): Iterator<T> = list.iterator()
+
+    override fun toString(): String = Buffer.toString(this)
 }
 
 

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MemoryBuffer.kt

@@ -22,6 +22,8 @@ public open class MemoryBuffer<T : Any>(protected val memory: Memory, protected
     override operator fun get(index: Int): T = reader.read(spec, spec.objectSize * index)
     override operator fun iterator(): Iterator<T> = (0 until size).asSequence().map { get(it) }.iterator()
 
+    override fun toString(): String = Buffer.toString(this)
+
     public companion object {
         public fun <T : Any> create(spec: MemorySpec<T>, size: Int): MemoryBuffer<T> =
             MemoryBuffer(Memory.allocate(size * spec.objectSize), spec)

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)
     }
 }

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) } }
         }
     }
 }

kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/RingBuffer.kt

@@ -69,6 +69,8 @@ public class RingBuffer<T>(
     @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 <reified T : Any> build(size: Int, empty: T): RingBuffer<T> {
             val buffer = MutableBuffer.auto(size) { empty } as MutableBuffer<T?>

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<T : Any, out A : Ring<T>, out M : org.ejml
     ): EjmlMatrix<T, M>
 
     public abstract override fun buildVector(size: Int, initializer: A.(Int) -> T): EjmlVector<T, M>
+
+    @UnstableKMathAPI
+    public fun EjmlMatrix<T, *>.inverse(): Structure2D<Double> = computeFeature(this, InverseMatrixFeature::class)?.inverse as Structure2D<Double>
 }

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.*
-import space.kscience.kmath.linear.LupDecompositionFeature
-import space.kscience.kmath.linear.computeFeature
 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) }
+    }
 }

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<Double>
 
 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<Double, DoubleField> =
-    LinearSpace.double.matrix(rowNum, colNum)
+    Double.algebra.linearSpace.matrix(rowNum, colNum)
 
 public fun Array<DoubleArray>.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<DoubleArray>.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>) -> 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<Double>) -
     }
 
 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
 

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<Double>.dot(other: Matrix<Double>): Matrix<Double> = LinearSpace.double.run {
+public infix fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = Double.algebra.linearSpace.run {
     this@dot dot other
 }

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 <T : Comparable<T>> PiecewisePolynomial<T>.integrate(algebra: Field<T>): PiecewisePolynomial<T> =
     PiecewisePolynomial(pieces.map { it.first to it.second.integrate(algebra) })

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] }

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<Double> {
+@JvmInline
+public value class ViktorBuffer(public val flatArray: F64FlatArray) : MutableBuffer<Double> {
     override val size: Int
         get() = flatArray.size
 
@@ -21,4 +23,6 @@ public class ViktorBuffer(public val flatArray: F64FlatArray) : MutableBuffer<Do
 
     override fun copy(): MutableBuffer<Double> = ViktorBuffer(flatArray.copy().flatten())
     override operator fun iterator(): Iterator<Double> = flatArray.data.iterator()
+
+    override fun toString(): String = Buffer.toString(this)
 }