kmath-for-real refactoring

CHANGELOG.md

@@ -28,6 +28,8 @@
 - `kmath-prob` renamed to `kmath-stat`
 - Grid generators moved to `kmath-for-real`
 - Use `Point<Double>` instead of specialized type in `kmath-for-real`
+- Optimized dot product for buffer matrices moved to `kmath-for-real`
+- EjmlMatrix context is an object
 
 ### Deprecated
 

examples/src/benchmarks/kotlin/kscience/kmath/structures/MultiplicationBenchmark.kt

@@ -0,0 +1,51 @@
+package kscience.kmath.structures
+
+import kotlinx.benchmark.Benchmark
+import kscience.kmath.commons.linear.CMMatrixContext
+import kscience.kmath.commons.linear.CMMatrixContext.dot
+import kscience.kmath.commons.linear.toCM
+import kscience.kmath.ejml.EjmlMatrixContext
+import kscience.kmath.ejml.toEjml
+import kscience.kmath.linear.real
+import kscience.kmath.operations.RealField
+import kscience.kmath.operations.invoke
+import org.openjdk.jmh.annotations.Scope
+import org.openjdk.jmh.annotations.State
+import kotlin.random.Random
+
+@State(Scope.Benchmark)
+class MultiplicationBenchmark {
+    companion object {
+        val random = Random(12224)
+        val dim = 1000
+
+        //creating invertible matrix
+        val matrix1 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+        val matrix2 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+
+        val cmMatrix1 = matrix1.toCM()
+        val cmMatrix2 = matrix2.toCM()
+
+        val ejmlMatrix1 = matrix1.toEjml()
+        val ejmlMatrix2 = matrix2.toEjml()
+    }
+
+    @Benchmark
+    fun commonsMathMultiplication() {
+        CMMatrixContext.invoke {
+            cmMatrix1 dot cmMatrix2
+        }
+    }
+
+    @Benchmark
+    fun ejmlMultiplication() {
+        EjmlMatrixContext.invoke {
+            ejmlMatrix1 dot ejmlMatrix2
+        }
+    }
+
+    @Benchmark
+    fun bufferedMultiplication() {
+        matrix1 dot matrix2
+    }
+}

examples/src/main/kotlin/kscience/kmath/commons/fit/fitWithAutoDiff.kt

@@ -70,6 +70,7 @@ fun main() {
     //minimize the chi^2 in given starting point. Derivatives are not required, they are already included.
     val result: OptimizationResult<Double> = chi2.minimize(a to 1.5, b to 0.9, c to 1.0)
 
+    //display a page with plot and numerical results
     val page = Plotly.page {
         plot {
             scatter {

examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt

@@ -1,10 +1,12 @@
 package kscience.kmath.linear
 
 import kscience.kmath.commons.linear.CMMatrixContext
+import kscience.kmath.commons.linear.CMMatrixContext.dot
 import kscience.kmath.commons.linear.inverse
 import kscience.kmath.commons.linear.toCM
 import kscience.kmath.ejml.EjmlMatrixContext
 import kscience.kmath.ejml.inverse
+import kscience.kmath.ejml.toEjml
 import kscience.kmath.operations.RealField
 import kscience.kmath.operations.invoke
 import kscience.kmath.structures.Matrix
@@ -40,7 +42,7 @@ fun main() {
     println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis")
 
     val ejmlTime = measureTimeMillis {
-        (EjmlMatrixContext(RealField)) {
+        EjmlMatrixContext {
             val km = matrix.toEjml()      //avoid overhead on conversion
             repeat(n) { inverse(km) }
         }

examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt

@@ -1,38 +0,0 @@
-package kscience.kmath.linear
-
-import kscience.kmath.commons.linear.CMMatrixContext
-import kscience.kmath.commons.linear.toCM
-import kscience.kmath.ejml.EjmlMatrixContext
-import kscience.kmath.operations.RealField
-import kscience.kmath.operations.invoke
-import kscience.kmath.structures.Matrix
-import kotlin.random.Random
-import kotlin.system.measureTimeMillis
-
-fun main() {
-    val random = Random(12224)
-    val dim = 1000
-    //creating invertible matrix
-    val matrix1 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
-    val matrix2 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
-
-//    //warmup
-//    matrix1 dot matrix2
-
-    CMMatrixContext {
-        val cmMatrix1 = matrix1.toCM()
-        val cmMatrix2 = matrix2.toCM()
-        val cmTime = measureTimeMillis { cmMatrix1 dot cmMatrix2 }
-        println("CM implementation time: $cmTime")
-    }
-
-    (EjmlMatrixContext(RealField)) {
-        val ejmlMatrix1 = matrix1.toEjml()
-        val ejmlMatrix2 = matrix2.toEjml()
-        val ejmlTime = measureTimeMillis { ejmlMatrix1 dot ejmlMatrix2 }
-        println("EJML implementation time: $ejmlTime")
-    }
-
-    val genericTime = measureTimeMillis { val res = matrix1 dot matrix2 }
-    println("Generic implementation time: $genericTime")
-}

kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt

@@ -9,7 +9,7 @@ import kscience.kmath.structures.*
  */
 public class BufferMatrixContext<T : Any, R : Ring<T>>(
     public override val elementContext: R,
-    private val bufferFactory: BufferFactory<T>
+    private val bufferFactory: BufferFactory<T>,
 ) : GenericMatrixContext<T, R> {
     public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): BufferMatrix<T> {
         val buffer = bufferFactory(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
@@ -29,8 +29,8 @@ public object RealMatrixContext : GenericMatrixContext<Double, RealField> {
     public override inline fun produce(
         rows: Int,
         columns: Int,
-        initializer: (i: Int, j: Int) -> Double
-    ): Matrix<Double> {
+        initializer: (i: Int, j: Int) -> Double,
+    ): BufferMatrix<Double> {
         val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
         return BufferMatrix(rows, columns, buffer)
     }
@@ -43,15 +43,15 @@ public class BufferMatrix<T : Any>(
     public override val rowNum: Int,
     public override val colNum: Int,
     public val buffer: Buffer<out T>,
-    public override val features: Set<MatrixFeature> = emptySet()
+    public override val features: Set<MatrixFeature> = emptySet(),
 ) : FeaturedMatrix<T> {
-    override val shape: IntArray
-        get() = intArrayOf(rowNum, colNum)
 
     init {
         require(buffer.size == rowNum * colNum) { "Dimension mismatch for matrix structure" }
     }
 
+    override val shape: IntArray get() = intArrayOf(rowNum, colNum)
+
     public override fun suggestFeature(vararg features: MatrixFeature): BufferMatrix<T> =
         BufferMatrix(rowNum, colNum, buffer, this.features + features)
 
@@ -86,28 +86,3 @@ public class BufferMatrix<T : Any>(
         else "Matrix(rowsNum = $rowNum, colNum = $colNum, features=$features)"
     }
 }
-
-/**
- * Optimized dot product for real matrices
- */
-public infix fun BufferMatrix<Double>.dot(other: BufferMatrix<Double>): BufferMatrix<Double> {
-    require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
-    val array = DoubleArray(this.rowNum * other.colNum)
-
-    //convert to array to insure there is not memory indirection
-    fun Buffer<out Double>.unsafeArray() = if (this is RealBuffer)
-        array
-    else
-        DoubleArray(size) { get(it) }
-
-    val a = this.buffer.unsafeArray()
-    val b = other.buffer.unsafeArray()
-
-    for (i in (0 until rowNum))
-        for (j in (0 until other.colNum))
-            for (k in (0 until colNum))
-                array[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j]
-
-    val buffer = RealBuffer(array)
-    return BufferMatrix(rowNum, other.colNum, buffer)
-}

kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt

@@ -27,9 +27,8 @@ public interface FeaturedMatrix<T : Any> : Matrix<T> {
 public inline fun Structure2D.Companion.real(
     rows: Int,
     columns: Int,
-    initializer: (Int, Int) -> Double
-): Matrix<Double> =
-    MatrixContext.real.produce(rows, columns, initializer)
+    initializer: (Int, Int) -> Double,
+): BufferMatrix<Double> = MatrixContext.real.produce(rows, columns, initializer)
 
 /**
  * Build a square matrix from given elements.
@@ -82,5 +81,3 @@ public fun <T : Any> Matrix<T>.transpose(): Matrix<T> {
         setOf(TransposedFeature(this))
     ) { i, j -> get(j, i) }
 }
-
-public infix fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = with(MatrixContext.real) { dot(other) }

kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt

@@ -1,5 +1,6 @@
 package kscience.kmath.linear
 
+import kscience.kmath.operations.invoke
 import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.NDStructure
 import kscience.kmath.structures.as2D
@@ -38,7 +39,7 @@ class MatrixTest {
         infix fun Matrix<Double>.pow(power: Int): Matrix<Double> {
             var res = this
             repeat(power - 1) {
-                res = res dot this
+                res = RealMatrixContext.invoke { res dot this@pow }
             }
             return res
         }

kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt

@@ -7,17 +7,18 @@ import kscience.kmath.operations.Space
 import kscience.kmath.operations.invoke
 import kscience.kmath.structures.Matrix
 
+/**
+ * Converts this matrix to EJML one.
+ */
+public fun Matrix<Double>.toEjml(): EjmlMatrix =
+    if (this is EjmlMatrix) this else EjmlMatrixContext.produce(rowNum, colNum) { i, j -> get(i, j) }
+
 /**
  * Represents context of basic operations operating with [EjmlMatrix].
  *
  * @author Iaroslav Postovalov
  */
-public class EjmlMatrixContext(private val space: Space<Double>) : MatrixContext<Double> {
-    /**
-     * Converts this matrix to EJML one.
-     */
-    public fun Matrix<Double>.toEjml(): EjmlMatrix =
-        if (this is EjmlMatrix) this else produce(rowNum, colNum) { i, j -> get(i, j) }
+public object EjmlMatrixContext : MatrixContext<Double> {
 
     /**
      * Converts this vector to EJML one.
@@ -47,11 +48,10 @@ public class EjmlMatrixContext(private val space: Space<Double>) : MatrixContext
         EjmlMatrix(toEjml().origin - b.toEjml().origin)
 
     public override fun multiply(a: Matrix<Double>, k: Number): EjmlMatrix =
-        produce(a.rowNum, a.colNum) { i, j -> space { a[i, j] * k } }
+        produce(a.rowNum, a.colNum) { i, j -> a[i, j] * k.toDouble() }
 
-    public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix = EjmlMatrix(toEjml().origin.scale(value))
-
-    public companion object
+    public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix =
+        EjmlMatrix(toEjml().origin.scale(value))
 }
 
 /**

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

@@ -6,10 +6,7 @@ import kscience.kmath.linear.VirtualMatrix
 import kscience.kmath.misc.UnstableKMathAPI
 import kscience.kmath.operations.invoke
 import kscience.kmath.operations.sum
-import kscience.kmath.structures.Buffer
-import kscience.kmath.structures.Matrix
-import kscience.kmath.structures.RealBuffer
-import kscience.kmath.structures.asIterable
+import kscience.kmath.structures.*
 import kotlin.math.pow
 
 /*
@@ -86,18 +83,6 @@ public operator fun Double.minus(matrix: RealMatrix): RealMatrix =
 //    row, col -> matrix[row, col] / this
 //}
 
-/*
- *  Per-element (!) square and power operations
- */
-
-public fun RealMatrix.square(): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { row, col ->
-    this[row, col].pow(2)
-}
-
-public fun RealMatrix.pow(n: Int): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { i, j ->
-    this[i, j].pow(n)
-}
-
 /*
  * Operations on two matrices (per-element!)
  */
@@ -157,3 +142,30 @@ public fun RealMatrix.sum(): Double = elements().map { (_, value) -> value }.sum
 public fun RealMatrix.min(): Double? = elements().map { (_, value) -> value }.minOrNull()
 public fun RealMatrix.max(): Double? = elements().map { (_, value) -> value }.maxOrNull()
 public fun RealMatrix.average(): Double = elements().map { (_, value) -> value }.average()
+
+public inline fun RealMatrix.map(transform: (Double) -> Double): RealMatrix =
+    MatrixContext.real.produce(rowNum, colNum) { i, j ->
+        transform(get(i, j))
+    }
+
+//extended operations
+
+public fun RealMatrix.pow(p: Double): RealMatrix = map { it.pow(p) }
+
+public fun RealMatrix.pow(p: Int): RealMatrix = map { it.pow(p) }
+
+public fun exp(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.exp(it) }
+
+public fun sqrt(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.sqrt(it) }
+
+public fun RealMatrix.square(): RealMatrix = map { it.pow(2) }
+
+public fun sin(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.sin(it) }
+
+public fun cos(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.cos(it) }
+
+public fun tan(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.tan(it) }
+
+public fun ln(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.ln(it) }
+
+public fun log10(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.log10(it) }

kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt

@@ -3,7 +3,6 @@ package kscience.kmath.real
 import kscience.kmath.linear.Point
 import kscience.kmath.operations.Norm
 import kscience.kmath.structures.Buffer
-import kscience.kmath.structures.RealBuffer
 import kscience.kmath.structures.asBuffer
 import kscience.kmath.structures.asIterable
 import kotlin.math.pow
@@ -11,17 +10,11 @@ import kotlin.math.sqrt
 
 public typealias RealVector = Point<Double>
 
-public inline fun RealVector(size: Int, init: (Int) -> Double): RealVector = RealBuffer(size, init)
-public fun RealVector(vararg doubles: Double): RealVector = RealBuffer(doubles)
-
-public fun DoubleArray.asVector(): RealVector = asBuffer()
-public fun List<Double>.asVector(): RealVector = asBuffer()
-
-
 public object VectorL2Norm : Norm<Point<out Number>, Double> {
     override fun norm(arg: Point<out Number>): Double = sqrt(arg.asIterable().sumByDouble(Number::toDouble))
 }
 
+public operator fun Buffer.Companion.invoke(vararg doubles: Double): RealVector = doubles.asBuffer()
 
 /**
  * Fill the vector of given [size] with given [value]
@@ -36,12 +29,6 @@ public inline fun RealVector.map(transform: (Double) -> Double): RealVector =
 public inline fun RealVector.mapIndexed(transform: (index: Int, value: Double) -> Double): RealVector =
     Buffer.real(size) { transform(it, get(it)) }
 
-public fun RealVector.pow(p: Double): RealVector = map { it.pow(p) }
-
-public fun RealVector.pow(p: Int): RealVector = map { it.pow(p) }
-
-public fun exp(vector: RealVector): RealVector = vector.map { kotlin.math.exp(it) }
-
 public operator fun RealVector.plus(other: RealVector): RealVector =
     mapIndexed { index, value -> value + other[index] }
 
@@ -71,3 +58,25 @@ public operator fun RealVector.div(other: RealVector): RealVector =
 public operator fun RealVector.div(number: Number): RealVector = map { it / number.toDouble() }
 
 public operator fun Number.div(vector: RealVector): RealVector = vector.map { toDouble() / it }
+
+//extended operations
+
+public fun RealVector.pow(p: Double): RealVector = map { it.pow(p) }
+
+public fun RealVector.pow(p: Int): RealVector = map { it.pow(p) }
+
+public fun exp(vector: RealVector): RealVector = vector.map { kotlin.math.exp(it) }
+
+public fun sqrt(vector: RealVector): RealVector = vector.map { kotlin.math.sqrt(it) }
+
+public fun RealVector.square(): RealVector = map { it.pow(2) }
+
+public fun sin(vector: RealVector): RealVector = vector.map { kotlin.math.sin(it) }
+
+public fun cos(vector: RealVector): RealVector = vector.map { kotlin.math.cos(it) }
+
+public fun tan(vector: RealVector): RealVector = vector.map { kotlin.math.tan(it) }
+
+public fun ln(vector: RealVector): RealVector = vector.map { kotlin.math.ln(it) }
+
+public fun log10(vector: RealVector): RealVector = vector.map { kotlin.math.log10(it) }

kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/dot.kt

@@ -0,0 +1,31 @@
+package kscience.kmath.real
+
+import kscience.kmath.linear.BufferMatrix
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.RealBuffer
+
+
+/**
+ * Optimized dot product for real matrices
+ */
+public infix fun BufferMatrix<Double>.dot(other: BufferMatrix<Double>): BufferMatrix<Double> {
+    require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
+    val resultArray = DoubleArray(this.rowNum * other.colNum)
+
+    //convert to array to insure there is no memory indirection
+    fun Buffer<out Double>.unsafeArray() = if (this is RealBuffer)
+        this.array
+    else
+        DoubleArray(size) { get(it) }
+
+    val a = this.buffer.unsafeArray()
+    val b = other.buffer.unsafeArray()
+
+    for (i in (0 until rowNum))
+        for (j in (0 until other.colNum))
+            for (k in (0 until colNum))
+                resultArray[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j]
+
+    val buffer = RealBuffer(resultArray)
+    return BufferMatrix(rowNum, other.colNum, buffer)
+}

kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt

@@ -1,34 +1,33 @@
 package kaceince.kmath.real
 
-import kscience.kmath.linear.MatrixContext
-import kscience.kmath.linear.asMatrix
-import kscience.kmath.linear.transpose
+import kscience.kmath.linear.*
 import kscience.kmath.operations.invoke
 import kscience.kmath.real.RealVector
 import kscience.kmath.real.plus
+import kscience.kmath.structures.Buffer
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
 internal class RealVectorTest {
     @Test
     fun testSum() {
-        val vector1 = RealVector(5) { it.toDouble() }
-        val vector2 = RealVector(5) { 5 - it.toDouble() }
+        val vector1 = Buffer.real(5) { it.toDouble() }
+        val vector2 =  Buffer.real(5) { 5 - it.toDouble() }
         val sum = vector1 + vector2
         assertEquals(5.0, sum[2])
     }
 
     @Test
     fun testVectorToMatrix() {
-        val vector = RealVector(5) { it.toDouble() }
+        val vector =  Buffer.real(5) { it.toDouble() }
         val matrix = vector.asMatrix()
         assertEquals(4.0, matrix[4, 0])
     }
 
     @Test
     fun testDot() {
-        val vector1 = RealVector(5) { it.toDouble() }
-        val vector2 = RealVector(5) { 5 - it.toDouble() }
+        val vector1 =  Buffer.real(5) { it.toDouble() }
+        val vector2 =  Buffer.real(5) { 5 - it.toDouble() }
         val matrix1 = vector1.asMatrix()
         val matrix2 = vector2.asMatrix().transpose()
         val product = MatrixContext.real { matrix1 dot matrix2 }

kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt

@@ -3,7 +3,6 @@ package kscience.kmath.histogram
 import kscience.kmath.linear.Point
 import kscience.kmath.operations.SpaceOperations
 import kscience.kmath.operations.invoke
-import kscience.kmath.real.asVector
 import kscience.kmath.structures.*
 import kotlin.math.floor
 
@@ -123,8 +122,8 @@ public class RealHistogram(
          *```
          */
         public fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): RealHistogram = RealHistogram(
-            ranges.map(ClosedFloatingPointRange<Double>::start).asVector(),
-            ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asVector()
+            ranges.map(ClosedFloatingPointRange<Double>::start).asBuffer(),
+            ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asBuffer()
         )
 
         /**

kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt

@@ -4,6 +4,8 @@ import kscience.kmath.histogram.RealHistogram
 import kscience.kmath.histogram.fill
 import kscience.kmath.histogram.put
 import kscience.kmath.real.RealVector
+import kscience.kmath.real.invoke
+import kscience.kmath.structures.Buffer
 import kotlin.random.Random
 import kotlin.test.*
 

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

@@ -1,8 +1,8 @@
 package kscience.kmath.histogram
 
 import kscience.kmath.real.RealVector
-import kscience.kmath.real.asVector
 import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.asBuffer
 import java.util.*
 import kotlin.math.floor
 
@@ -16,7 +16,7 @@ public class UnivariateBin(
     //TODO add weighting
     public override val value: Number get() = counter.sum()
 
-    public override val center: RealVector get() = doubleArrayOf(position).asVector()
+    public override val center: RealVector get() = doubleArrayOf(position).asBuffer()
     public override val dimension: Int get() = 1
 
     public operator fun contains(value: Double): Boolean = value in (position - size / 2)..(position + size / 2)