DoubleBuffer algebra refactoring

examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt

@@ -10,20 +10,21 @@ import kotlinx.html.h3
 import space.kscience.kmath.commons.expressions.DSProcessor
 import space.kscience.kmath.commons.optimization.CMOptimizer
 import space.kscience.kmath.distributions.NormalDistribution
-import space.kscience.kmath.expressions.binding
 import space.kscience.kmath.expressions.chiSquaredExpression
 import space.kscience.kmath.expressions.symbol
-import space.kscience.kmath.optimization.*
+import space.kscience.kmath.operations.asIterable
+import space.kscience.kmath.operations.toList
+import space.kscience.kmath.optimization.FunctionOptimizationTarget
+import space.kscience.kmath.optimization.optimizeWith
+import space.kscience.kmath.optimization.resultPoint
+import space.kscience.kmath.optimization.resultValue
 import space.kscience.kmath.real.DoubleVector
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
 import space.kscience.kmath.stat.RandomGenerator
-import space.kscience.kmath.structures.asIterable
-import space.kscience.kmath.structures.toList
 import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
 import space.kscience.plotly.models.TraceValues
-import kotlin.math.abs
 import kotlin.math.pow
 import kotlin.math.sqrt
 

examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt

@@ -13,6 +13,8 @@ import space.kscience.kmath.distributions.NormalDistribution
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.expressions.binding
 import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.operations.asIterable
+import space.kscience.kmath.operations.toList
 import space.kscience.kmath.optimization.QowOptimizer
 import space.kscience.kmath.optimization.chiSquaredOrNull
 import space.kscience.kmath.optimization.fitWith
@@ -20,8 +22,6 @@ import space.kscience.kmath.optimization.resultPoint
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
 import space.kscience.kmath.stat.RandomGenerator
-import space.kscience.kmath.structures.asIterable
-import space.kscience.kmath.structures.toList
 import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
 import kotlin.math.abs

examples/src/main/kotlin/space/kscience/kmath/functions/interpolateSquare.kt

@@ -8,8 +8,8 @@ package space.kscience.kmath.functions
 import space.kscience.kmath.interpolation.SplineInterpolator
 import space.kscience.kmath.interpolation.interpolatePolynomials
 import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
-import space.kscience.kmath.structures.map
 import space.kscience.plotly.Plotly
 import space.kscience.plotly.UnstablePlotlyAPI
 import space.kscience.plotly.makeFile

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/transform/Transformations.kt

@@ -10,10 +10,13 @@ import kotlinx.coroutines.flow.Flow
 import kotlinx.coroutines.flow.map
 import org.apache.commons.math3.transform.*
 import space.kscience.kmath.complex.Complex
+import space.kscience.kmath.operations.SuspendBufferTransform
 import space.kscience.kmath.streaming.chunked
 import space.kscience.kmath.streaming.spread
-import space.kscience.kmath.structures.*
+import space.kscience.kmath.structures.Buffer
-
+import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.VirtualBuffer
+import space.kscience.kmath.structures.asBuffer
 
 
 /**

kmath-commons/src/test/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt

@@ -13,11 +13,11 @@ import space.kscience.kmath.expressions.Symbol.Companion.x
 import space.kscience.kmath.expressions.Symbol.Companion.y
 import space.kscience.kmath.expressions.chiSquaredExpression
 import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.operations.map
 import space.kscience.kmath.optimization.*
 import space.kscience.kmath.stat.RandomGenerator
 import space.kscience.kmath.structures.DoubleBuffer
 import space.kscience.kmath.structures.asBuffer
-import space.kscience.kmath.structures.map
 import kotlin.math.pow
 import kotlin.test.Test
 

kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/specialExpressions.kt

@@ -6,8 +6,8 @@
 package space.kscience.kmath.expressions
 
 import space.kscience.kmath.operations.ExtendedField
+import space.kscience.kmath.operations.asIterable
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.asIterable
 import space.kscience.kmath.structures.indices
 import kotlin.jvm.JvmName
 

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/DoubleLinearSpace.kt

@@ -9,6 +9,7 @@ import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.nd.DoubleFieldND
 import space.kscience.kmath.nd.as2D
 import space.kscience.kmath.nd.asND
+import space.kscience.kmath.operations.DoubleBufferOperations
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
@@ -89,17 +90,21 @@ public object DoubleLinearSpace : LinearSpace<Double, DoubleField> {
         asND().map { it * value }.as2D()
     }
 
-    public override fun Point<Double>.plus(other: Point<Double>): DoubleBuffer = DoubleBuffer(size) {
+    public override fun Point<Double>.plus(other: Point<Double>): DoubleBuffer = DoubleBufferOperations.run {
-        get(it) + other[it]
+        this@plus + other
     }
 
-    public override fun Point<Double>.minus(other: Point<Double>): DoubleBuffer = DoubleBuffer(size) {
+    public override fun Point<Double>.minus(other: Point<Double>): DoubleBuffer = DoubleBufferOperations.run {
-        get(it) - other[it]
+        this@minus - other
     }
 
-    public override fun Point<Double>.times(value: Double): DoubleBuffer = DoubleBuffer(size) { i -> get(i) * value }
+    public override fun Point<Double>.times(value: Double): DoubleBuffer = DoubleBufferOperations.run {
+        scale(this@times, value)
+    }
 
-    public operator fun Point<Double>.div(value: Double): DoubleBuffer = DoubleBuffer(size) { i -> get(i) / value }
+    public operator fun Point<Double>.div(value: Double): DoubleBuffer = DoubleBufferOperations.run {
+        scale(this@div, 1.0 / value)
+    }
 
     public override fun Double.times(v: Point<Double>): DoubleBuffer = v * this
 

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

@@ -6,10 +6,10 @@
 package space.kscience.kmath.nd
 
 import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.operations.asSequence
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.MutableBuffer
 import space.kscience.kmath.structures.asMutableBuffer
-import space.kscience.kmath.structures.asSequence
 import kotlin.jvm.JvmInline
 
 /**

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

@@ -0,0 +1,130 @@
+/*
+ * Copyright 2018-2021 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
+ */
+
+package space.kscience.kmath.operations
+
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.DoubleBuffer
+
+/**
+ * [ExtendedField] over [DoubleBuffer].
+ *
+ * @property size the size of buffers to operate on.
+ */
+public class DoubleBufferField(public val size: Int) : ExtendedField<Buffer<Double>>, DoubleBufferOperations() {
+    override val zero: Buffer<Double> by lazy { DoubleBuffer(size) { 0.0 } }
+    override val one: Buffer<Double> by lazy { DoubleBuffer(size) { 1.0 } }
+
+    override fun sinh(arg: Buffer<Double>): DoubleBuffer = super<DoubleBufferOperations>.sinh(arg)
+
+    override fun cosh(arg: Buffer<Double>): DoubleBuffer = super<DoubleBufferOperations>.cosh(arg)
+
+    override fun tanh(arg: Buffer<Double>): DoubleBuffer = super<DoubleBufferOperations>.tanh(arg)
+
+    override fun asinh(arg: Buffer<Double>): DoubleBuffer = super<DoubleBufferOperations>.asinh(arg)
+
+    override fun acosh(arg: Buffer<Double>): DoubleBuffer = super<DoubleBufferOperations>.acosh(arg)
+
+    override fun atanh(arg: Buffer<Double>): DoubleBuffer= super<DoubleBufferOperations>.atanh(arg)
+
+    //    override fun number(value: Number): Buffer<Double> = DoubleBuffer(size) { value.toDouble() }
+//
+//    override fun Buffer<Double>.unaryMinus(): Buffer<Double> = DoubleBufferOperations.run {
+//        -this@unaryMinus
+//    }
+//
+//    override fun add(a: Buffer<Double>, b: Buffer<Double>): DoubleBuffer {
+//        require(a.size == size) { "The buffer size ${a.size} does not match context size $size" }
+//        return DoubleBufferOperations.add(a, b)
+//    }
+//
+
+//
+//    override fun multiply(a: Buffer<Double>, b: Buffer<Double>): DoubleBuffer {
+//        require(a.size == size) { "The buffer size ${a.size} does not match context size $size" }
+//        return DoubleBufferOperations.multiply(a, b)
+//    }
+//
+//    override fun divide(a: Buffer<Double>, b: Buffer<Double>): DoubleBuffer {
+//        require(a.size == size) { "The buffer size ${a.size} does not match context size $size" }
+//        return DoubleBufferOperations.divide(a, b)
+//    }
+//
+//    override fun sin(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.sin(arg)
+//    }
+//
+//    override fun cos(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.cos(arg)
+//    }
+//
+//    override fun tan(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.tan(arg)
+//    }
+//
+//    override fun asin(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.asin(arg)
+//    }
+//
+//    override fun acos(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.acos(arg)
+//    }
+//
+//    override fun atan(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.atan(arg)
+//    }
+//
+//    override fun sinh(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.sinh(arg)
+//    }
+//
+//    override fun cosh(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.cosh(arg)
+//    }
+//
+//    override fun tanh(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.tanh(arg)
+//    }
+//
+//    override fun asinh(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.asinh(arg)
+//    }
+//
+//    override fun acosh(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.acosh(arg)
+//    }
+//
+//    override fun atanh(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.atanh(arg)
+//    }
+//
+//    override fun power(arg: Buffer<Double>, pow: Number): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.power(arg, pow)
+//    }
+//
+//    override fun exp(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.exp(arg)
+//    }
+//
+//    override fun ln(arg: Buffer<Double>): DoubleBuffer {
+//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
+//        return DoubleBufferOperations.ln(arg)
+//    }
+
+}

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

@@ -8,13 +8,13 @@ package space.kscience.kmath.operations
 import space.kscience.kmath.linear.Point
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.fold
+
 import kotlin.math.*
 
 /**
  * [ExtendedFieldOperations] over [DoubleBuffer].
  */
-public object DoubleBufferOperations : ExtendedFieldOperations<Buffer<Double>> {
+public abstract class DoubleBufferOperations : ExtendedFieldOperations<Buffer<Double>>, Norm<Buffer<Double>, Double> {
     override fun Buffer<Double>.unaryMinus(): DoubleBuffer = if (this is DoubleBuffer) {
         DoubleBuffer(size) { -array[it] }
     } else {
@@ -32,6 +32,21 @@ public object DoubleBufferOperations : ExtendedFieldOperations<Buffer<Double>> {
             DoubleBuffer(DoubleArray(a.size) { aArray[it] + bArray[it] })
         } else DoubleBuffer(DoubleArray(a.size) { a[it] + b[it] })
     }
+
+    override fun Buffer<Double>.plus(b: Buffer<Double>): DoubleBuffer = add(this, b)
+
+    override fun Buffer<Double>.minus(b: Buffer<Double>): DoubleBuffer {
+        require(b.size == this.size) {
+            "The size of the first buffer ${this.size} should be the same as for second one: ${b.size} "
+        }
+
+        return if (this is DoubleBuffer && b is DoubleBuffer) {
+            val aArray = this.array
+            val bArray = b.array
+            DoubleBuffer(DoubleArray(this.size) { aArray[it] - bArray[it] })
+        } else DoubleBuffer(DoubleArray(this.size) { this[it] - b[it] })
+    }
+
     //
 //    override fun multiply(a: Buffer<Double>, k: Number): RealBuffer {
 //        val kValue = k.toDouble()
@@ -159,128 +174,21 @@ public object DoubleBufferOperations : ExtendedFieldOperations<Buffer<Double>> {
     override fun ln(arg: Buffer<Double>): DoubleBuffer = if (arg is DoubleBuffer) {
         val array = arg.array
         DoubleBuffer(DoubleArray(arg.size) { ln(array[it]) })
-    } else
+    } else {
         DoubleBuffer(DoubleArray(arg.size) { ln(arg[it]) })
     }
 
+    override fun norm(arg: Buffer<Double>): Double = DoubleL2Norm.norm(arg)
+
+    override fun scale(a: Buffer<Double>, value: Double): DoubleBuffer = if (a is DoubleBuffer) {
+        val aArray = a.array
+        DoubleBuffer(DoubleArray(a.size) { aArray[it] * value })
+    } else DoubleBuffer(DoubleArray(a.size) { a[it] * value })
+
+    public companion object : DoubleBufferOperations()
+}
+
 public object DoubleL2Norm : Norm<Point<Double>, Double> {
     override fun norm(arg: Point<Double>): Double = sqrt(arg.fold(0.0) { acc: Double, d: Double -> acc + d.pow(2) })
 }
 
-/**
- * [ExtendedField] over [DoubleBuffer].
- *
- * @property size the size of buffers to operate on.
- */
-@Deprecated("To be replaced by generic BufferAlgebra")
-public class DoubleBufferField(public val size: Int) : ExtendedField<Buffer<Double>>, Norm<Buffer<Double>, Double> {
-    override val zero: Buffer<Double> by lazy { DoubleBuffer(size) { 0.0 } }
-    override val one: Buffer<Double> by lazy { DoubleBuffer(size) { 1.0 } }
-
-    override fun number(value: Number): Buffer<Double> = DoubleBuffer(size) { value.toDouble() }
-
-    override fun Buffer<Double>.unaryMinus(): Buffer<Double> = DoubleBufferOperations.run {
-        -this@unaryMinus
-    }
-
-    override fun add(a: Buffer<Double>, b: Buffer<Double>): DoubleBuffer {
-        require(a.size == size) { "The buffer size ${a.size} does not match context size $size" }
-        return DoubleBufferOperations.add(a, b)
-    }
-
-    override fun scale(a: Buffer<Double>, value: Double): DoubleBuffer {
-        require(a.size == size) { "The buffer size ${a.size} does not match context size $size" }
-
-        return if (a is DoubleBuffer) {
-            val aArray = a.array
-            DoubleBuffer(DoubleArray(a.size) { aArray[it] * value })
-        } else DoubleBuffer(DoubleArray(a.size) { a[it] * value })
-    }
-
-    override fun multiply(a: Buffer<Double>, b: Buffer<Double>): DoubleBuffer {
-        require(a.size == size) { "The buffer size ${a.size} does not match context size $size" }
-        return DoubleBufferOperations.multiply(a, b)
-    }
-
-    override fun divide(a: Buffer<Double>, b: Buffer<Double>): DoubleBuffer {
-        require(a.size == size) { "The buffer size ${a.size} does not match context size $size" }
-        return DoubleBufferOperations.divide(a, b)
-    }
-
-    override fun sin(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.sin(arg)
-    }
-
-    override fun cos(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.cos(arg)
-    }
-
-    override fun tan(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.tan(arg)
-    }
-
-    override fun asin(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.asin(arg)
-    }
-
-    override fun acos(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.acos(arg)
-    }
-
-    override fun atan(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.atan(arg)
-    }
-
-    override fun sinh(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.sinh(arg)
-    }
-
-    override fun cosh(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.cosh(arg)
-    }
-
-    override fun tanh(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.tanh(arg)
-    }
-
-    override fun asinh(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.asinh(arg)
-    }
-
-    override fun acosh(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.acosh(arg)
-    }
-
-    override fun atanh(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.atanh(arg)
-    }
-
-    override fun power(arg: Buffer<Double>, pow: Number): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.power(arg, pow)
-    }
-
-    override fun exp(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.exp(arg)
-    }
-
-    override fun ln(arg: Buffer<Double>): DoubleBuffer {
-        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-        return DoubleBufferOperations.ln(arg)
-    }
-
-    override fun norm(arg: Buffer<Double>): Double  = DoubleL2Norm.norm(arg)
-}

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

@@ -3,9 +3,10 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
  */
 
-package space.kscience.kmath.structures
+package space.kscience.kmath.operations
 
 import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.structures.*
 
 /**
  * Typealias for buffer transformations.

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

@@ -14,7 +14,8 @@ public interface ExtendedFieldOperations<T> :
     FieldOperations<T>,
     TrigonometricOperations<T>,
     PowerOperations<T>,
-    ExponentialOperations<T> {
+    ExponentialOperations<T>,
+    ScaleOperations<T>  {
     override fun tan(arg: T): T = sin(arg) / cos(arg)
     override fun tanh(arg: T): T = sinh(arg) / cosh(arg)
 
@@ -41,7 +42,7 @@ public interface ExtendedFieldOperations<T> :
 /**
  * Advanced Number-like field that implements basic operations.
  */
-public interface ExtendedField<T> : ExtendedFieldOperations<T>, Field<T>, NumericAlgebra<T>, ScaleOperations<T> {
+public interface ExtendedField<T> : ExtendedFieldOperations<T>, Field<T>, NumericAlgebra<T>{
     override fun sinh(arg: T): T = (exp(arg) - exp(-arg)) / 2.0
     override fun cosh(arg: T): T = (exp(arg) + exp(-arg)) / 2.0
     override fun tanh(arg: T): T = (exp(arg) - exp(-arg)) / (exp(-arg) + exp(arg))

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

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.structures
 
+import space.kscience.kmath.operations.asSequence
 import kotlin.jvm.JvmInline
 import kotlin.reflect.KClass
 
@@ -49,7 +50,8 @@ public interface Buffer<out T> {
 
     public companion object {
 
-        public fun toString(buffer: Buffer<*>): String = buffer.asSequence().joinToString(prefix = "[", separator = ", ", postfix = "]")
+        public fun toString(buffer: Buffer<*>): String =
+            buffer.asSequence().joinToString(prefix = "[", separator = ", ", postfix = "]")
 
         /**
          * Check the element-by-element match of content of two buffers.

kmath-coroutines/src/jvmTest/kotlin/space/kscience/kmath/streaming/RingBufferTest.kt

@@ -7,7 +7,7 @@ package space.kscience.kmath.streaming
 
 import kotlinx.coroutines.flow.*
 import kotlinx.coroutines.runBlocking
-import space.kscience.kmath.structures.asSequence
+import space.kscience.kmath.operations.asSequence
 import kotlin.test.Test
 import kotlin.test.assertEquals
 

kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt

@@ -0,0 +1,995 @@
+/*
+ * Copyright 2018-2021 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
+ */
+
+/* This file is generated with buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt */
+
+package space.kscience.kmath.ejml
+
+import org.ejml.data.*
+import org.ejml.dense.row.CommonOps_DDRM
+import org.ejml.dense.row.CommonOps_FDRM
+import org.ejml.dense.row.factory.DecompositionFactory_DDRM
+import org.ejml.dense.row.factory.DecompositionFactory_FDRM
+import org.ejml.sparse.FillReducing
+import org.ejml.sparse.csc.CommonOps_DSCC
+import org.ejml.sparse.csc.CommonOps_FSCC
+import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC
+import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC
+import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC
+import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC
+import space.kscience.kmath.linear.*
+import space.kscience.kmath.linear.Matrix
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.nd.StructureFeature
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.FloatField
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.FloatBuffer
+import kotlin.reflect.KClass
+import kotlin.reflect.cast
+
+/**
+ * [EjmlVector] specialization for [Double].
+ */
+public class EjmlDoubleVector<out M : DMatrix>(override val origin: M) : EjmlVector<Double, M>(origin) {
+    init {
+        require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" }
+    }
+
+    override operator fun get(index: Int): Double = origin[0, index]
+}
+
+/**
+ * [EjmlVector] specialization for [Float].
+ */
+public class EjmlFloatVector<out M : FMatrix>(override val origin: M) : EjmlVector<Float, M>(origin) {
+    init {
+        require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" }
+    }
+
+    override operator fun get(index: Int): Float = origin[0, index]
+}
+
+/**
+ * [EjmlMatrix] specialization for [Double].
+ */
+public class EjmlDoubleMatrix<out M : DMatrix>(override val origin: M) : EjmlMatrix<Double, M>(origin) {
+    override operator fun get(i: Int, j: Int): Double = origin[i, j]
+}
+
+/**
+ * [EjmlMatrix] specialization for [Float].
+ */
+public class EjmlFloatMatrix<out M : FMatrix>(override val origin: M) : EjmlMatrix<Float, M>(origin) {
+    override operator fun get(i: Int, j: Int): Float = origin[i, j]
+}
+
+/**
+ * [EjmlLinearSpace] implementation based on [CommonOps_DDRM], [DecompositionFactory_DDRM] operations and
+ * [DMatrixRMaj] matrices.
+ */
+public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, DoubleField, DMatrixRMaj>() {
+    /**
+     * The [DoubleField] reference.
+     */
+    override val elementAlgebra: DoubleField get() = DoubleField
+
+    @Suppress("UNCHECKED_CAST")
+    override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixRMaj> = when {
+        this is EjmlDoubleMatrix<*> && origin is DMatrixRMaj -> this as EjmlDoubleMatrix<DMatrixRMaj>
+        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixRMaj> = when {
+        this is EjmlDoubleVector<*> && origin is DMatrixRMaj -> this as EjmlDoubleVector<DMatrixRMaj>
+        else -> EjmlDoubleVector(DMatrixRMaj(size, 1).also {
+            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
+        })
+    }
+
+    override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: DoubleField.(i: Int, j: Int) -> Double,
+    ): EjmlDoubleMatrix<DMatrixRMaj> = DMatrixRMaj(rows, columns).also {
+        (0 until rows).forEach { row ->
+            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
+        }
+    }.wrapMatrix()
+
+    override fun buildVector(
+        size: Int,
+        initializer: DoubleField.(Int) -> Double,
+    ): EjmlDoubleVector<DMatrixRMaj> = EjmlDoubleVector(DMatrixRMaj(size, 1).also {
+        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
+    })
+
+    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
+    private fun <T : DMatrix> T.wrapVector() = EjmlDoubleVector(this)
+
+    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
+
+    override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.mult(toEjml().origin, other.toEjml().origin, out)
+        return out.wrapMatrix()
+    }
+
+    override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.mult(toEjml().origin, vector.toEjml().origin, out)
+        return out.wrapVector()
+    }
+
+    override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+
+        CommonOps_DDRM.add(
+            elementAlgebra.one, 
+            toEjml().origin, 
+            elementAlgebra { -one },
+            other.toEjml().origin, 
+            out,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    override operator fun Matrix<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixRMaj> {
+        val res = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.scale(value, toEjml().origin, res)
+        return res.wrapMatrix()
+    }
+
+    override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixRMaj> {
+        val res = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.changeSign(toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+        
+        CommonOps_DDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin, 
+            out,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+
+        CommonOps_DDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin,
+            out,
+        )
+
+        return out.wrapVector()
+    }
+
+    override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+        val out = DMatrixRMaj(1, 1)
+
+        CommonOps_DDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra { -one },
+            other.toEjml().origin,
+            out,
+        )
+
+        return out.wrapVector()
+    }
+
+    override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> = m * this
+
+    override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixRMaj> {
+        val res = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.scale(value, toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixRMaj> = v * this
+
+    @UnstableKMathAPI
+    override fun <F : StructureFeature> computeFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
+        structure.getFeature(type)?.let { return it }
+        val origin = structure.toEjml().origin
+
+        return when (type) {
+            InverseMatrixFeature::class -> object : InverseMatrixFeature<Double> {
+                override val inverse: Matrix<Double> by lazy {
+                    val res = origin.copy()
+                    CommonOps_DDRM.invert(res)
+                    res.wrapMatrix()
+                }
+            }
+
+            DeterminantFeature::class -> object : DeterminantFeature<Double> {
+                override val determinant: Double by lazy { CommonOps_DDRM.det(origin) }
+            }
+
+            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Double> {
+                private val svd by lazy {
+                    DecompositionFactory_DDRM.svd(origin.numRows, origin.numCols, true, true, false)
+                        .apply { decompose(origin.copy()) }
+                }
+
+                override val u: Matrix<Double> by lazy { svd.getU(null, false).wrapMatrix() }
+                override val s: Matrix<Double> by lazy { svd.getW(null).wrapMatrix() }
+                override val v: Matrix<Double> by lazy { svd.getV(null, false).wrapMatrix() }
+                override val singularValues: Point<Double> by lazy { DoubleBuffer(svd.singularValues) }
+            }
+
+            QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
+                private val qr by lazy {
+                    DecompositionFactory_DDRM.qr().apply { decompose(origin.copy()) }
+                }
+
+                override val q: Matrix<Double> by lazy {
+                    qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature)
+                }
+
+                override val r: Matrix<Double> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) }
+            }
+
+            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
+                override val l: Matrix<Double> by lazy {
+                    val cholesky =
+                        DecompositionFactory_DDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
+
+                    cholesky.getT(null).wrapMatrix().withFeature(LFeature)
+                }
+            }
+
+            LupDecompositionFeature::class -> object : LupDecompositionFeature<Double> {
+                private val lup by lazy {
+                    DecompositionFactory_DDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
+                }
+
+                override val l: Matrix<Double> by lazy {
+                    lup.getLower(null).wrapMatrix().withFeature(LFeature)
+                }
+
+                override val u: Matrix<Double> by lazy {
+                    lup.getUpper(null).wrapMatrix().withFeature(UFeature)
+                }
+
+                override val p: Matrix<Double> by lazy { lup.getRowPivot(null).wrapMatrix() }
+            }
+
+            else -> null
+        }?.let(type::cast)
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p matrix.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
+        val res = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res)
+        return res.wrapMatrix()
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p vector.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
+        val res = DMatrixRMaj(1, 1)
+        CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res)
+        return EjmlDoubleVector(res)
+    }
+}
+
+/**
+ * [EjmlLinearSpace] implementation based on [CommonOps_FDRM], [DecompositionFactory_FDRM] operations and
+ * [FMatrixRMaj] matrices.
+ */
+public object EjmlLinearSpaceFDRM : EjmlLinearSpace<Float, FloatField, FMatrixRMaj>() {
+    /**
+     * The [FloatField] reference.
+     */
+    override val elementAlgebra: FloatField get() = FloatField
+
+    @Suppress("UNCHECKED_CAST")
+    override fun Matrix<Float>.toEjml(): EjmlFloatMatrix<FMatrixRMaj> = when {
+        this is EjmlFloatMatrix<*> && origin is FMatrixRMaj -> this as EjmlFloatMatrix<FMatrixRMaj>
+        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    override fun Point<Float>.toEjml(): EjmlFloatVector<FMatrixRMaj> = when {
+        this is EjmlFloatVector<*> && origin is FMatrixRMaj -> this as EjmlFloatVector<FMatrixRMaj>
+        else -> EjmlFloatVector(FMatrixRMaj(size, 1).also {
+            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
+        })
+    }
+
+    override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: FloatField.(i: Int, j: Int) -> Float,
+    ): EjmlFloatMatrix<FMatrixRMaj> = FMatrixRMaj(rows, columns).also {
+        (0 until rows).forEach { row ->
+            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
+        }
+    }.wrapMatrix()
+
+    override fun buildVector(
+        size: Int,
+        initializer: FloatField.(Int) -> Float,
+    ): EjmlFloatVector<FMatrixRMaj> = EjmlFloatVector(FMatrixRMaj(size, 1).also {
+        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
+    })
+
+    private fun <T : FMatrix> T.wrapMatrix() = EjmlFloatMatrix(this)
+    private fun <T : FMatrix> T.wrapVector() = EjmlFloatVector(this)
+
+    override fun Matrix<Float>.unaryMinus(): Matrix<Float> = this * elementAlgebra { -one }
+
+    override fun Matrix<Float>.dot(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.mult(toEjml().origin, other.toEjml().origin, out)
+        return out.wrapMatrix()
+    }
+
+    override fun Matrix<Float>.dot(vector: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.mult(toEjml().origin, vector.toEjml().origin, out)
+        return out.wrapVector()
+    }
+
+    override operator fun Matrix<Float>.minus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+
+        CommonOps_FDRM.add(
+            elementAlgebra.one, 
+            toEjml().origin, 
+            elementAlgebra { -one },
+            other.toEjml().origin, 
+            out,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    override operator fun Matrix<Float>.times(value: Float): EjmlFloatMatrix<FMatrixRMaj> {
+        val res = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.scale(value, toEjml().origin, res)
+        return res.wrapMatrix()
+    }
+
+    override fun Point<Float>.unaryMinus(): EjmlFloatVector<FMatrixRMaj> {
+        val res = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.changeSign(toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    override fun Matrix<Float>.plus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+        
+        CommonOps_FDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin, 
+            out,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    override fun Point<Float>.plus(other: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+
+        CommonOps_FDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin,
+            out,
+        )
+
+        return out.wrapVector()
+    }
+
+    override fun Point<Float>.minus(other: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
+        val out = FMatrixRMaj(1, 1)
+
+        CommonOps_FDRM.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra { -one },
+            other.toEjml().origin,
+            out,
+        )
+
+        return out.wrapVector()
+    }
+
+    override fun Float.times(m: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> = m * this
+
+    override fun Point<Float>.times(value: Float): EjmlFloatVector<FMatrixRMaj> {
+        val res = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.scale(value, toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixRMaj> = v * this
+
+    @UnstableKMathAPI
+    override fun <F : StructureFeature> computeFeature(structure: Matrix<Float>, type: KClass<out F>): F? {
+        structure.getFeature(type)?.let { return it }
+        val origin = structure.toEjml().origin
+
+        return when (type) {
+            InverseMatrixFeature::class -> object : InverseMatrixFeature<Float> {
+                override val inverse: Matrix<Float> by lazy {
+                    val res = origin.copy()
+                    CommonOps_FDRM.invert(res)
+                    res.wrapMatrix()
+                }
+            }
+
+            DeterminantFeature::class -> object : DeterminantFeature<Float> {
+                override val determinant: Float by lazy { CommonOps_FDRM.det(origin) }
+            }
+
+            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Float> {
+                private val svd by lazy {
+                    DecompositionFactory_FDRM.svd(origin.numRows, origin.numCols, true, true, false)
+                        .apply { decompose(origin.copy()) }
+                }
+
+                override val u: Matrix<Float> by lazy { svd.getU(null, false).wrapMatrix() }
+                override val s: Matrix<Float> by lazy { svd.getW(null).wrapMatrix() }
+                override val v: Matrix<Float> by lazy { svd.getV(null, false).wrapMatrix() }
+                override val singularValues: Point<Float> by lazy { FloatBuffer(svd.singularValues) }
+            }
+
+            QRDecompositionFeature::class -> object : QRDecompositionFeature<Float> {
+                private val qr by lazy {
+                    DecompositionFactory_FDRM.qr().apply { decompose(origin.copy()) }
+                }
+
+                override val q: Matrix<Float> by lazy {
+                    qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature)
+                }
+
+                override val r: Matrix<Float> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) }
+            }
+
+            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Float> {
+                override val l: Matrix<Float> by lazy {
+                    val cholesky =
+                        DecompositionFactory_FDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
+
+                    cholesky.getT(null).wrapMatrix().withFeature(LFeature)
+                }
+            }
+
+            LupDecompositionFeature::class -> object : LupDecompositionFeature<Float> {
+                private val lup by lazy {
+                    DecompositionFactory_FDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
+                }
+
+                override val l: Matrix<Float> by lazy {
+                    lup.getLower(null).wrapMatrix().withFeature(LFeature)
+                }
+
+                override val u: Matrix<Float> by lazy {
+                    lup.getUpper(null).wrapMatrix().withFeature(UFeature)
+                }
+
+                override val p: Matrix<Float> by lazy { lup.getRowPivot(null).wrapMatrix() }
+            }
+
+            else -> null
+        }?.let(type::cast)
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p matrix.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Float>, b: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
+        val res = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.solve(FMatrixRMaj(a.toEjml().origin), FMatrixRMaj(b.toEjml().origin), res)
+        return res.wrapMatrix()
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p vector.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Float>, b: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
+        val res = FMatrixRMaj(1, 1)
+        CommonOps_FDRM.solve(FMatrixRMaj(a.toEjml().origin), FMatrixRMaj(b.toEjml().origin), res)
+        return EjmlFloatVector(res)
+    }
+}
+
+/**
+ * [EjmlLinearSpace] implementation based on [CommonOps_DSCC], [DecompositionFactory_DSCC] operations and
+ * [DMatrixSparseCSC] matrices.
+ */
+public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, DoubleField, DMatrixSparseCSC>() {
+    /**
+     * The [DoubleField] reference.
+     */
+    override val elementAlgebra: DoubleField get() = DoubleField
+
+    @Suppress("UNCHECKED_CAST")
+    override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixSparseCSC> = when {
+        this is EjmlDoubleMatrix<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleMatrix<DMatrixSparseCSC>
+        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixSparseCSC> = when {
+        this is EjmlDoubleVector<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleVector<DMatrixSparseCSC>
+        else -> EjmlDoubleVector(DMatrixSparseCSC(size, 1).also {
+            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
+        })
+    }
+
+    override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: DoubleField.(i: Int, j: Int) -> Double,
+    ): EjmlDoubleMatrix<DMatrixSparseCSC> = DMatrixSparseCSC(rows, columns).also {
+        (0 until rows).forEach { row ->
+            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
+        }
+    }.wrapMatrix()
+
+    override fun buildVector(
+        size: Int,
+        initializer: DoubleField.(Int) -> Double,
+    ): EjmlDoubleVector<DMatrixSparseCSC> = EjmlDoubleVector(DMatrixSparseCSC(size, 1).also {
+        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
+    })
+
+    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
+    private fun <T : DMatrix> T.wrapVector() = EjmlDoubleVector(this)
+
+    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
+
+    override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.mult(toEjml().origin, other.toEjml().origin, out)
+        return out.wrapMatrix()
+    }
+
+    override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.mult(toEjml().origin, vector.toEjml().origin, out)
+        return out.wrapVector()
+    }
+
+    override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+
+        CommonOps_DSCC.add(
+            elementAlgebra.one, 
+            toEjml().origin, 
+            elementAlgebra { -one },
+            other.toEjml().origin, 
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    override operator fun Matrix<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixSparseCSC> {
+        val res = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.scale(value, toEjml().origin, res)
+        return res.wrapMatrix()
+    }
+
+    override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixSparseCSC> {
+        val res = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.changeSign(toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+        
+        CommonOps_DSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin, 
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+
+        CommonOps_DSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin,
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapVector()
+    }
+
+    override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+        val out = DMatrixSparseCSC(1, 1)
+
+        CommonOps_DSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra { -one },
+            other.toEjml().origin,
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapVector()
+    }
+
+    override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> = m * this
+
+    override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixSparseCSC> {
+        val res = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.scale(value, toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> = v * this
+
+    @UnstableKMathAPI
+    override fun <F : StructureFeature> computeFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
+        structure.getFeature(type)?.let { return it }
+        val origin = structure.toEjml().origin
+
+        return when (type) {
+            QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
+                private val qr by lazy {
+                    DecompositionFactory_DSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
+                }
+
+                override val q: Matrix<Double> by lazy {
+                    qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature)
+                }
+
+                override val r: Matrix<Double> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) }
+            }
+
+            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
+                override val l: Matrix<Double> by lazy {
+                    val cholesky =
+                        DecompositionFactory_DSCC.cholesky().apply { decompose(origin.copy()) }
+
+                    (cholesky.getT(null) as DMatrix).wrapMatrix().withFeature(LFeature)
+                }
+            }
+
+            LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
+                LUDecompositionFeature<Double>, DeterminantFeature<Double>, InverseMatrixFeature<Double> {
+                private val lu by lazy {
+                    DecompositionFactory_DSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
+                }
+
+                override val l: Matrix<Double> by lazy {
+                    lu.getLower(null).wrapMatrix().withFeature(LFeature)
+                }
+
+                override val u: Matrix<Double> by lazy {
+                    lu.getUpper(null).wrapMatrix().withFeature(UFeature)
+                }
+
+                override val inverse: Matrix<Double> by lazy {
+                    var a = origin
+                    val inverse = DMatrixRMaj(1, 1)
+                    val solver = LinearSolverFactory_DSCC.lu(FillReducing.NONE)
+                    if (solver.modifiesA()) a = a.copy()
+                    val i = CommonOps_DDRM.identity(a.numRows)
+                    solver.solve(i, inverse)
+                    inverse.wrapMatrix()
+                }
+
+                override val determinant: Double by lazy { elementAlgebra.number(lu.computeDeterminant().real) }
+            }
+
+            else -> null
+        }?.let(type::cast)
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p matrix.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
+        val res = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res)
+        return res.wrapMatrix()
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p vector.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
+        val res = DMatrixSparseCSC(1, 1)
+        CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res)
+        return EjmlDoubleVector(res)
+    }
+}
+
+/**
+ * [EjmlLinearSpace] implementation based on [CommonOps_FSCC], [DecompositionFactory_FSCC] operations and
+ * [FMatrixSparseCSC] matrices.
+ */
+public object EjmlLinearSpaceFSCC : EjmlLinearSpace<Float, FloatField, FMatrixSparseCSC>() {
+    /**
+     * The [FloatField] reference.
+     */
+    override val elementAlgebra: FloatField get() = FloatField
+
+    @Suppress("UNCHECKED_CAST")
+    override fun Matrix<Float>.toEjml(): EjmlFloatMatrix<FMatrixSparseCSC> = when {
+        this is EjmlFloatMatrix<*> && origin is FMatrixSparseCSC -> this as EjmlFloatMatrix<FMatrixSparseCSC>
+        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    override fun Point<Float>.toEjml(): EjmlFloatVector<FMatrixSparseCSC> = when {
+        this is EjmlFloatVector<*> && origin is FMatrixSparseCSC -> this as EjmlFloatVector<FMatrixSparseCSC>
+        else -> EjmlFloatVector(FMatrixSparseCSC(size, 1).also {
+            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
+        })
+    }
+
+    override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: FloatField.(i: Int, j: Int) -> Float,
+    ): EjmlFloatMatrix<FMatrixSparseCSC> = FMatrixSparseCSC(rows, columns).also {
+        (0 until rows).forEach { row ->
+            (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
+        }
+    }.wrapMatrix()
+
+    override fun buildVector(
+        size: Int,
+        initializer: FloatField.(Int) -> Float,
+    ): EjmlFloatVector<FMatrixSparseCSC> = EjmlFloatVector(FMatrixSparseCSC(size, 1).also {
+        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
+    })
+
+    private fun <T : FMatrix> T.wrapMatrix() = EjmlFloatMatrix(this)
+    private fun <T : FMatrix> T.wrapVector() = EjmlFloatVector(this)
+
+    override fun Matrix<Float>.unaryMinus(): Matrix<Float> = this * elementAlgebra { -one }
+
+    override fun Matrix<Float>.dot(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.mult(toEjml().origin, other.toEjml().origin, out)
+        return out.wrapMatrix()
+    }
+
+    override fun Matrix<Float>.dot(vector: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.mult(toEjml().origin, vector.toEjml().origin, out)
+        return out.wrapVector()
+    }
+
+    override operator fun Matrix<Float>.minus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+
+        CommonOps_FSCC.add(
+            elementAlgebra.one, 
+            toEjml().origin, 
+            elementAlgebra { -one },
+            other.toEjml().origin, 
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    override operator fun Matrix<Float>.times(value: Float): EjmlFloatMatrix<FMatrixSparseCSC> {
+        val res = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.scale(value, toEjml().origin, res)
+        return res.wrapMatrix()
+    }
+
+    override fun Point<Float>.unaryMinus(): EjmlFloatVector<FMatrixSparseCSC> {
+        val res = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.changeSign(toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    override fun Matrix<Float>.plus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+        
+        CommonOps_FSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin, 
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapMatrix()
+    }
+
+    override fun Point<Float>.plus(other: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+
+        CommonOps_FSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra.one,
+            other.toEjml().origin,
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapVector()
+    }
+
+    override fun Point<Float>.minus(other: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
+        val out = FMatrixSparseCSC(1, 1)
+
+        CommonOps_FSCC.add(
+            elementAlgebra.one,
+            toEjml().origin,
+            elementAlgebra { -one },
+            other.toEjml().origin,
+            out,
+            null, 
+            null,
+        )
+
+        return out.wrapVector()
+    }
+
+    override fun Float.times(m: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> = m * this
+
+    override fun Point<Float>.times(value: Float): EjmlFloatVector<FMatrixSparseCSC> {
+        val res = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.scale(value, toEjml().origin, res)
+        return res.wrapVector()
+    }
+
+    override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> = v * this
+
+    @UnstableKMathAPI
+    override fun <F : StructureFeature> computeFeature(structure: Matrix<Float>, type: KClass<out F>): F? {
+        structure.getFeature(type)?.let { return it }
+        val origin = structure.toEjml().origin
+
+        return when (type) {
+            QRDecompositionFeature::class -> object : QRDecompositionFeature<Float> {
+                private val qr by lazy {
+                    DecompositionFactory_FSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
+                }
+
+                override val q: Matrix<Float> by lazy {
+                    qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature)
+                }
+
+                override val r: Matrix<Float> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) }
+            }
+
+            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Float> {
+                override val l: Matrix<Float> by lazy {
+                    val cholesky =
+                        DecompositionFactory_FSCC.cholesky().apply { decompose(origin.copy()) }
+
+                    (cholesky.getT(null) as FMatrix).wrapMatrix().withFeature(LFeature)
+                }
+            }
+
+            LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
+                LUDecompositionFeature<Float>, DeterminantFeature<Float>, InverseMatrixFeature<Float> {
+                private val lu by lazy {
+                    DecompositionFactory_FSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
+                }
+
+                override val l: Matrix<Float> by lazy {
+                    lu.getLower(null).wrapMatrix().withFeature(LFeature)
+                }
+
+                override val u: Matrix<Float> by lazy {
+                    lu.getUpper(null).wrapMatrix().withFeature(UFeature)
+                }
+
+                override val inverse: Matrix<Float> by lazy {
+                    var a = origin
+                    val inverse = FMatrixRMaj(1, 1)
+                    val solver = LinearSolverFactory_FSCC.lu(FillReducing.NONE)
+                    if (solver.modifiesA()) a = a.copy()
+                    val i = CommonOps_FDRM.identity(a.numRows)
+                    solver.solve(i, inverse)
+                    inverse.wrapMatrix()
+                }
+
+                override val determinant: Float by lazy { elementAlgebra.number(lu.computeDeterminant().real) }
+            }
+
+            else -> null
+        }?.let(type::cast)
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p matrix.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Float>, b: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
+        val res = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.solve(FMatrixSparseCSC(a.toEjml().origin), FMatrixSparseCSC(b.toEjml().origin), res)
+        return res.wrapMatrix()
+    }
+
+    /**
+     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [b]*.
+     *
+     * @param a the base matrix.
+     * @param b n by p vector.
+     * @return the solution for *x* that is n by p.
+     */
+    public fun solve(a: Matrix<Float>, b: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
+        val res = FMatrixSparseCSC(1, 1)
+        CommonOps_FSCC.solve(FMatrixSparseCSC(a.toEjml().origin), FMatrixSparseCSC(b.toEjml().origin), res)
+        return EjmlFloatVector(res)
+    }
+}
+

kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealMatrix.kt

@@ -13,9 +13,9 @@ import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.asIterable
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.asIterable
 import kotlin.math.pow
 
 /*

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegratorRuleFactory.kt

@@ -5,10 +5,10 @@
 
 package space.kscience.kmath.integration
 
+import space.kscience.kmath.operations.map
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
 import space.kscience.kmath.structures.asBuffer
-import space.kscience.kmath.structures.map
 import kotlin.jvm.Synchronized
 import kotlin.math.ulp
 import kotlin.native.concurrent.ThreadLocal

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/SplineIntegrator.kt

@@ -12,14 +12,10 @@ import space.kscience.kmath.interpolation.SplineInterpolator
 import space.kscience.kmath.interpolation.interpolatePolynomials
 import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
-import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.*
-import space.kscience.kmath.operations.Field
-import space.kscience.kmath.operations.invoke
-import space.kscience.kmath.operations.sum
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
 import space.kscience.kmath.structures.MutableBufferFactory
-import space.kscience.kmath.structures.map
 
 /**
  * Compute analytical indefinite integral of this [PiecewisePolynomial], keeping all intervals intact

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector2DTest.kt

@@ -1,9 +1,8 @@
 package space.kscience.kmath.geometry
 
-import space.kscience.kmath.structures.asSequence
+import space.kscience.kmath.operations.toList
-import space.kscience.kmath.structures.toList
-import kotlin.test.assertEquals
 import kotlin.test.Test
+import kotlin.test.assertEquals
 
 internal class Vector2DTest {
     private val vector = Vector2D(1.0, -7.999)

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector3DTest.kt

@@ -1,6 +1,6 @@
 package space.kscience.kmath.geometry
 
-import space.kscience.kmath.structures.toList
+import space.kscience.kmath.operations.toList
 import kotlin.test.Test
 import kotlin.test.assertEquals
 

kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt

@@ -7,8 +7,8 @@ package space.kscience.kmath.histogram
 
 import space.kscience.kmath.domains.UnivariateDomain
 import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.asSequence
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.asSequence
 
 
 @UnstableKMathAPI

kmath-jupyter/src/main/kotlin/space/kscience/kmath/jupyter/KMathJupyter.kt

@@ -21,9 +21,9 @@ import space.kscience.kmath.expressions.MST
 import space.kscience.kmath.expressions.MstRing
 import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.nd.Structure2D
+import space.kscience.kmath.operations.asSequence
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.asSequence
 
 /**
  * A function for conversion of number to MST for pretty print

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/Median.kt

@@ -5,8 +5,8 @@
 
 package space.kscience.kmath.stat
 
+import space.kscience.kmath.operations.asSequence
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.asSequence
 
 /**
  * Non-composable median

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt

@@ -9,9 +9,9 @@ import space.kscience.kmath.nd.MutableStructure1D
 import space.kscience.kmath.nd.MutableStructure2D
 import space.kscience.kmath.nd.as1D
 import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.operations.asSequence
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.VirtualBuffer
-import space.kscience.kmath.structures.asSequence
 import space.kscience.kmath.tensors.core.BufferedTensor
 import space.kscience.kmath.tensors.core.DoubleTensor
 import space.kscience.kmath.tensors.core.DoubleTensorAlgebra

kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/utils.kt

@@ -6,6 +6,7 @@
 package space.kscience.kmath.tensors.core.internal
 
 import space.kscience.kmath.nd.as1D
+import space.kscience.kmath.operations.toMutableList
 import space.kscience.kmath.samplers.GaussianSampler
 import space.kscience.kmath.stat.RandomGenerator
 import space.kscience.kmath.structures.*