WIP vector space refactor

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt

@@ -3,12 +3,13 @@ package space.kscience.kmath.commons.linear
 import org.apache.commons.math3.linear.*
 import space.kscience.kmath.linear.*
 import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.nd.NDStructure
 import space.kscience.kmath.operations.RealField
 import space.kscience.kmath.structures.RealBuffer
 import kotlin.reflect.KClass
 import kotlin.reflect.cast
 
-public inline class CMMatrix(public val origin: RealMatrix) : Matrix<Double> {
+public class CMMatrix(public val origin: RealMatrix) : Matrix<Double> {
     public override val rowNum: Int get() = origin.rowDimension
     public override val colNum: Int get() = origin.columnDimension
 
@@ -51,12 +52,17 @@ public inline class CMMatrix(public val origin: RealMatrix) : Matrix<Double> {
     }?.let(type::cast)
 
     public override operator fun get(i: Int, j: Int): Double = origin.getEntry(i, j)
+
+    override fun equals(other: Any?): Boolean {
+        if (this === other) return true
+        if (other !is NDStructure<*>) return false
+        return NDStructure.contentEquals(this, other)
     }
 
+    override fun hashCode(): Int = origin.hashCode()
+}
 
-public fun RealMatrix.asMatrix(): CMMatrix = CMMatrix(this)
+public inline class CMVector(public val origin: RealVector) : Vector<Double> {
-
-public class CMVector(public val origin: RealVector) : Vector<Double> {
     public override val size: Int get() = origin.dimension
 
     public override operator fun get(index: Int): Double = origin.getEntry(index)
@@ -64,16 +70,17 @@ public class CMVector(public val origin: RealVector) : Vector<Double> {
     public override operator fun iterator(): Iterator<Double> = origin.toArray().iterator()
 }
 
-public fun Point<Double>.toCM(): CMVector = if (this is CMVector) this else {
-    val array = DoubleArray(size) { this[it] }
-    CMVector(ArrayRealVector(array))
-}
-
 public fun RealVector.toPoint(): CMVector = CMVector(this)
 
 public object CMLinearSpace : LinearSpace<Double, RealField> {
-    public override fun buildMatrix(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): CMMatrix {
+    override val elementAlgebra: RealField get() = RealField
-        val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } }
+
+    public override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: RealField.(i: Int, j: Int) -> Double,
+    ): CMMatrix {
+        val array = Array(rows) { i -> DoubleArray(columns) { j -> RealField.initializer(i, j) } }
         return CMMatrix(Array2DRowRealMatrix(array))
     }
 
@@ -83,31 +90,50 @@ public object CMLinearSpace : LinearSpace<Double, RealField> {
         else -> {
             //TODO add feature analysis
             val array = Array(rowNum) { i -> DoubleArray(colNum) { j -> get(i, j) } }
-            CMMatrix(Array2DRowRealMatrix(array))
+            Array2DRowRealMatrix(array).wrap()
         }
     }
 
+    public fun Vector<Double>.toCM(): CMVector = if (this is CMVector) this else {
+        val array = DoubleArray(size) { this[it] }
+        ArrayRealVector(array).wrap()
+    }
+
+    private fun RealMatrix.wrap(): CMMatrix = CMMatrix(this)
+    private fun RealVector.wrap(): CMVector = CMVector(this)
+
+    override fun buildVector(size: Int, initializer: RealField.(Int) -> Double): Vector<Double> =
+        ArrayRealVector(DoubleArray(size) { RealField.initializer(it) }).wrap()
+
+    override fun Matrix<Double>.plus(other: Matrix<Double>): CMMatrix =
+        toCM().origin.add(other.toCM().origin).wrap()
+
+    override fun Vector<Double>.plus(other: Vector<Double>): CMVector =
+        toCM().origin.add(other.toCM().origin).wrap()
+
+    override fun Vector<Double>.minus(other: Vector<Double>): CMVector =
+        toCM().origin.subtract(other.toCM().origin).wrap()
 
     public override fun Matrix<Double>.dot(other: Matrix<Double>): CMMatrix =
-        CMMatrix(toCM().origin.multiply(other.toCM().origin))
+        toCM().origin.multiply(other.toCM().origin).wrap()
 
     public override fun Matrix<Double>.dot(vector: Vector<Double>): CMVector =
-        CMVector(toCM().origin.preMultiply(vector.toCM().origin))
+        toCM().origin.preMultiply(vector.toCM().origin).wrap()
 
-    public override operator fun Matrix<Double>.unaryMinus(): CMMatrix =
+    public override operator fun Matrix<Double>.minus(other: Matrix<Double>): CMMatrix =
-        buildMatrix(rowNum, colNum) { i, j -> -get(i, j) }
+        toCM().origin.subtract(other.toCM().origin).wrap()
-
-    public override fun add(a: Matrix<Double>, b: Matrix<Double>): CMMatrix =
-        CMMatrix(a.toCM().origin.multiply(b.toCM().origin))
-
-    public override operator fun Matrix<Double>.minus(b: Matrix<Double>): CMMatrix =
-        CMMatrix(toCM().origin.subtract(b.toCM().origin))
-
-//    public override fun multiply(a: Matrix<Double>, k: Number): CMMatrix =
-//        CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble()))
 
     public override operator fun Matrix<Double>.times(value: Double): CMMatrix =
-        buildMatrix(rowNum, colNum) { i, j -> get(i, j) * value }
+        toCM().origin.scalarMultiply(value).wrap()
+
+    override fun Double.times(m: Matrix<Double>): CMMatrix =
+        m * this
+
+    override fun Vector<Double>.times(value: Double): CMVector =
+        toCM().origin.mapMultiply(value).wrap()
+
+    override fun Double.times(v: Vector<Double>): CMVector =
+        v * this
 }
 
 public operator fun CMMatrix.plus(other: CMMatrix): CMMatrix =

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

@@ -14,12 +14,7 @@ import kotlin.reflect.KClass
  */
 public typealias Matrix<T> = Structure2D<T>
 
-/**
+public typealias Vector<T> = Point<T>
- * Alias for [Structure1D] with more familiar name.
- *
- * @param T the type of items.
- */
-public typealias Vector<T> = Structure1D<T>
 
 /**
  * Basic operations on matrices and vectors. Operates on [Matrix].
@@ -183,12 +178,13 @@ public interface LinearSpace<T : Any, A : Ring<T>> {
             override fun buildMatrix(
                 rows: Int, columns: Int,
                 initializer: A.(i: Int, j: Int) -> T,
-            ): Matrix<T> = NDStructure.buffered(intArrayOf(rows, columns)) { (i, j) ->
+            ): Matrix<T> = NDStructure.buffered(intArrayOf(rows, columns), bufferFactory) { (i, j) ->
                 algebra.initializer(i, j)
             }.as2D()
-
         }
 
+        public val real: LinearSpace<Double, RealField> = buffered(RealField, Buffer.Companion::real)
+
         /**
          * Automatic buffered matrix, unboxed if it is possible
          */

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

@@ -12,7 +12,7 @@ public fun <T : Any> LinearSpace<T, Ring<T>>.matrix(rows: Int, columns: Int, var
 
 @UnstableKMathAPI
 public fun <T : Any> LinearSpace<T, Ring<T>>.vector(vararg elements: T): Vector<T> {
-    return buildVector(elements.size, elements::get)
+    return buildVector(elements.size) { elements[it] }
 }
 
 public inline fun <T : Any> LinearSpace<T, Ring<T>>.row(

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

@@ -1,88 +1,86 @@
 package space.kscience.kmath.linear
 
-import space.kscience.kmath.operations.RealField
+//public object RealLinearSpace:
-import space.kscience.kmath.operations.ScaleOperations
-import space.kscience.kmath.structures.RealBuffer
-
-public object RealLinearSpace : LinearSpace<Double, RealField>, ScaleOperations<Matrix<Double>> {
-
-    override val elementAlgebra: RealField get() = RealField
-
-    public override fun buildMatrix(
-        rows: Int,
-        columns: Int,
-        initializer: (i: Int, j: Int) -> Double,
-    ): Matrix<Double> {
-        val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
-        return  BufferMatrix(rows, columns, buffer)
-    }
-
-    public fun Matrix<Double>.toBufferMatrix(): BufferMatrix<Double> = if (this is BufferMatrix) this else {
-        buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    public fun one(rows: Int, columns: Int): Matrix<Double> = VirtualMatrix(rows, columns) { i, j ->
-        if (i == j) 1.0 else 0.0
-    } + DiagonalFeature
-
-    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = buildMatrix(rowNum, colNum) { i, j -> -get(i, j) }
-
-    public override infix fun Matrix<Double>.dot(other: Matrix<Double>): BufferMatrix<Double> {
-        require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
-        val bufferMatrix = toBufferMatrix()
-        val otherBufferMatrix = other.toBufferMatrix()
-        return buildMatrix(rowNum, other.colNum) { i, j ->
-            var res = 0.0
-            for (l in 0 until colNum) {
-                res += bufferMatrix[i, l] * otherBufferMatrix[l, j]
-            }
-            res
-        }
-    }
-
-    public override infix fun Matrix<Double>.dot(vector: Point<Double>): Point<Double> {
-        require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
-        val bufferMatrix = toBufferMatrix()
-        return RealBuffer(rowNum) { i ->
-            var res = 0.0
-            for (j in 0 until colNum) {
-                res += bufferMatrix[i, j] * vector[j]
-            }
-            res
-        }
-    }
-
-    override fun add(a: Matrix<Double>, b: Matrix<Double>): BufferMatrix<Double> {
-        require(a.rowNum == b.rowNum) { "Row number mismatch in matrix addition. Left side: ${a.rowNum}, right side: ${b.rowNum}" }
-        require(a.colNum == b.colNum) { "Column number mismatch in matrix addition. Left side: ${a.colNum}, right side: ${b.colNum}" }
-        val aBufferMatrix = a.toBufferMatrix()
-        val bBufferMatrix = b.toBufferMatrix()
-        return buildMatrix(a.rowNum, a.colNum) { i, j ->
-            aBufferMatrix[i, j] + bBufferMatrix[i, j]
-        }
-    }
-
-    override fun scale(a: Matrix<Double>, value: Double): BufferMatrix<Double> {
-        val bufferMatrix = a.toBufferMatrix()
-        return buildMatrix(a.rowNum, a.colNum) { i, j -> bufferMatrix[i, j] * value }
-    }
-
-    override fun Matrix<Double>.times(value: Double): BufferMatrix<Double> = scale(this, value)
 
+//public object RealLinearSpace : LinearSpace<Double, RealField>, ScaleOperations<Matrix<Double>> {
 //
-//    override fun multiply(a: Matrix<Double>, k: Number): BufferMatrix<Double> {
+//    override val elementAlgebra: RealField get() = RealField
-//        val aBufferMatrix = a.toBufferMatrix()
+//
-//        return produce(a.rowNum, a.colNum) { i, j -> aBufferMatrix[i, j] * k.toDouble() }
+//    public override fun buildMatrix(
+//        rows: Int,
+//        columns: Int,
+//        initializer: (i: Int, j: Int) -> Double,
+//    ): Matrix<Double> {
+//        val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
+//        return  BufferMatrix(rows, columns, buffer)
 //    }
 //
-//    override fun divide(a: Matrix<Double>, k: Number): BufferMatrix<Double> {
+//    public fun Matrix<Double>.toBufferMatrix(): BufferMatrix<Double> = if (this is BufferMatrix) this else {
-//        val aBufferMatrix = a.toBufferMatrix()
+//        buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-//        return produce(a.rowNum, a.colNum) { i, j -> aBufferMatrix[i, j] / k.toDouble() }
+//    }
+//
+//    public fun one(rows: Int, columns: Int): Matrix<Double> = VirtualMatrix(rows, columns) { i, j ->
+//        if (i == j) 1.0 else 0.0
+//    } + DiagonalFeature
+//
+//    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = buildMatrix(rowNum, colNum) { i, j -> -get(i, j) }
+//
+//    public override infix fun Matrix<Double>.dot(other: Matrix<Double>): BufferMatrix<Double> {
+//        require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
+//        val bufferMatrix = toBufferMatrix()
+//        val otherBufferMatrix = other.toBufferMatrix()
+//        return buildMatrix(rowNum, other.colNum) { i, j ->
+//            var res = 0.0
+//            for (l in 0 until colNum) {
+//                res += bufferMatrix[i, l] * otherBufferMatrix[l, j]
+//            }
+//            res
+//        }
+//    }
+//
+//    public override infix fun Matrix<Double>.dot(vector: Point<Double>): Point<Double> {
+//        require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
+//        val bufferMatrix = toBufferMatrix()
+//        return RealBuffer(rowNum) { i ->
+//            var res = 0.0
+//            for (j in 0 until colNum) {
+//                res += bufferMatrix[i, j] * vector[j]
+//            }
+//            res
+//        }
+//    }
+//
+//    override fun add(a: Matrix<Double>, b: Matrix<Double>): BufferMatrix<Double> {
+//        require(a.rowNum == b.rowNum) { "Row number mismatch in matrix addition. Left side: ${a.rowNum}, right side: ${b.rowNum}" }
+//        require(a.colNum == b.colNum) { "Column number mismatch in matrix addition. Left side: ${a.colNum}, right side: ${b.colNum}" }
+//        val aBufferMatrix = a.toBufferMatrix()
+//        val bBufferMatrix = b.toBufferMatrix()
+//        return buildMatrix(a.rowNum, a.colNum) { i, j ->
+//            aBufferMatrix[i, j] + bBufferMatrix[i, j]
+//        }
+//    }
+//
+//    override fun scale(a: Matrix<Double>, value: Double): BufferMatrix<Double> {
+//        val bufferMatrix = a.toBufferMatrix()
+//        return buildMatrix(a.rowNum, a.colNum) { i, j -> bufferMatrix[i, j] * value }
+//    }
+//
+//    override fun Matrix<Double>.times(value: Double): BufferMatrix<Double> = scale(this, value)
+//
+////
+////    override fun multiply(a: Matrix<Double>, k: Number): BufferMatrix<Double> {
+////        val aBufferMatrix = a.toBufferMatrix()
+////        return produce(a.rowNum, a.colNum) { i, j -> aBufferMatrix[i, j] * k.toDouble() }
+////    }
+////
+////    override fun divide(a: Matrix<Double>, k: Number): BufferMatrix<Double> {
+////        val aBufferMatrix = a.toBufferMatrix()
+////        return produce(a.rowNum, a.colNum) { i, j -> aBufferMatrix[i, j] / k.toDouble() }
+////    }
 //}
-}
 
 
-/**
+///**
- * Partially optimized real-valued matrix
+// * Partially optimized real-valued matrix
- */
+// */
-public val LinearSpace.Companion.real: RealLinearSpace get() = RealLinearSpace
+//public val LinearSpace.Companion.real: RealLinearSpace get() = RealLinearSpace

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

@@ -52,7 +52,7 @@ public interface NDStructure<T> {
      * optimize operations and performance. If the feature is not present, null is defined.
      */
     @UnstableKMathAPI
-    public fun <T : Any> getFeature(type: KClass<T>): T? = null
+    public fun <F : Any> getFeature(type: KClass<F>): F? = null
 
     public companion object {
         /**

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

@@ -4,6 +4,7 @@ import org.ejml.simple.SimpleMatrix
 import space.kscience.kmath.linear.*
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.getFeature
+import space.kscience.kmath.operations.RealField
 import space.kscience.kmath.operations.ScaleOperations
 
 /**
@@ -11,7 +12,9 @@ import space.kscience.kmath.operations.ScaleOperations
  *
  * @author Iaroslav Postovalov
  */
-public object EjmlLinearSpace : LinearSpace<Double, EjmlMatrix>, ScaleOperations<Matrix<Double>> {
+public object EjmlLinearSpace : LinearSpace<Double, RealField>, ScaleOperations<Matrix<Double>> {
+
+    override val elementAlgebra: RealField get() = RealField
 
     /**
      * Converts this matrix to EJML one.
@@ -25,43 +28,74 @@ public object EjmlLinearSpace : LinearSpace<Double, EjmlMatrix>, ScaleOperations
     /**
      * Converts this vector to EJML one.
      */
-    public fun Point<Double>.toEjml(): EjmlVector =
+    public fun Vector<Double>.toEjml(): EjmlVector = when (this) {
-        if (this is EjmlVector) this else EjmlVector(SimpleMatrix(size, 1).also {
+        is EjmlVector -> this
+        else -> EjmlVector(SimpleMatrix(size, 1).also {
             (0 until it.numRows()).forEach { row -> it[row, 0] = get(row) }
         })
+    }
 
-    override fun buildMatrix(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): EjmlMatrix =
+    override fun buildMatrix(rows: Int, columns: Int, initializer: RealField.(i: Int, j: Int) -> Double): EjmlMatrix =
         EjmlMatrix(SimpleMatrix(rows, columns).also {
             (0 until rows).forEach { row ->
-                (0 until columns).forEach { col -> it[row, col] = initializer(row, col) }
+                (0 until columns).forEach { col -> it[row, col] = RealField.initializer(row, col) }
             }
         })
 
-    override fun buildVector(size: Int, initializer: (Int) -> Double): Point<Double> =
+    override fun buildVector(size: Int, initializer: RealField.(Int) -> Double): Vector<Double> =
         EjmlVector(SimpleMatrix(size, 1).also {
             (0 until it.numRows()).forEach { row -> it[row, 0] = initializer(row) }
         })
 
+    private fun SimpleMatrix.wrapMatrix() = EjmlMatrix(this)
+    private fun SimpleMatrix.wrapVector() = EjmlVector(this)
 
     override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * (-1)
 
     public override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlMatrix =
         EjmlMatrix(toEjml().origin.mult(other.toEjml().origin))
 
-    public override fun Matrix<Double>.dot(vector: Point<Double>): EjmlVector =
+    public override fun Matrix<Double>.dot(vector: Vector<Double>): EjmlVector =
         EjmlVector(toEjml().origin.mult(vector.toEjml().origin))
 
-    public override fun add(a: Matrix<Double>, b: Matrix<Double>): EjmlMatrix =
+    public override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlMatrix =
-        EjmlMatrix(a.toEjml().origin + b.toEjml().origin)
+        EjmlMatrix(toEjml().origin - other.toEjml().origin)
-
-    public override operator fun Matrix<Double>.minus(b: Matrix<Double>): EjmlMatrix =
-        EjmlMatrix(toEjml().origin - b.toEjml().origin)
 
     public override fun scale(a: Matrix<Double>, value: Double): EjmlMatrix =
-        buildMatrix(a.rowNum, a.colNum) { i, j -> a[i, j] * value }
+        a.toEjml().origin.scale(value).wrapMatrix()
 
     public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix =
         EjmlMatrix(toEjml().origin.scale(value))
+
+    override fun Vector<Double>.unaryMinus(): EjmlVector =
+        toEjml().origin.negative().wrapVector()
+
+    override fun Matrix<Double>.plus(other: Matrix<Double>): Matrix<Double> {
+        TODO("Not yet implemented")
+    }
+
+    override fun Vector<Double>.plus(other: Vector<Double>): Vector<Double> {
+        TODO("Not yet implemented")
+    }
+
+    override fun Vector<Double>.minus(other: Vector<Double>): Vector<Double> {
+        TODO("Not yet implemented")
+    }
+
+    override fun Double.times(m: Matrix<Double>): Matrix<Double> {
+        TODO("Not yet implemented")
+    }
+
+    override fun Vector<Double>.times(value: Double): Vector<Double> {
+        TODO("Not yet implemented")
+    }
+
+    override fun Double.times(v: Vector<Double>): Vector<Double> {
+        TODO("Not yet implemented")
+    }
+
+    public override fun add(a: Matrix<Double>, b: Matrix<Double>): EjmlMatrix =
+        EjmlMatrix(a.toEjml().origin + b.toEjml().origin)
 }
 
 /**

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

@@ -85,7 +85,7 @@ public class EjmlMatrix(public val origin: SimpleMatrix) : Matrix<Double> {
 
     override fun equals(other: Any?): Boolean {
         if (this === other) return true
-        if (other !is Matrix<*>) return false
+        if (other !is NDStructure<*>) return false
         return NDStructure.contentEquals(this, other)
     }