diff --git a/CHANGELOG.md b/CHANGELOG.md index 5d13a360f..89e02d3b1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ - Better trigonometric and hyperbolic functions for `AutoDiffField` (https://github.com/mipt-npm/kmath/pull/140). - Automatic README generation for features (#139) - Native support for `memory`, `core` and `dimensions` +- `kmath-ejml` to supply EJML SimpleMatrix wrapper. ### Changed - Package changed from `scientifik` to `kscience.kmath`. diff --git a/README.md b/README.md index 8bc85bf2b..708bd8eb1 100644 --- a/README.md +++ b/README.md @@ -53,7 +53,9 @@ can be used for a wide variety of purposes from high performance calculations to * **Commons-math wrapper** It is planned to gradually wrap most parts of [Apache commons-math](http://commons.apache.org/proper/commons-math/) library in Kotlin code and maybe rewrite some parts to better suit the Kotlin programming paradigm, however there is no fixed roadmap for that. Feel free to submit a feature request if you want something to be done first. - + +* **EJML wrapper** Provides EJML `SimpleMatrix` wrapper consistent with the core matrix structures. + ## Planned features * **Messaging** A mathematical notation to support multi-language and multi-node communication for mathematical tasks. diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts index a4b14fe96..900da966b 100644 --- a/examples/build.gradle.kts +++ b/examples/build.gradle.kts @@ -26,6 +26,7 @@ dependencies { implementation(project(":kmath-prob")) implementation(project(":kmath-viktor")) implementation(project(":kmath-dimensions")) + implementation(project(":kmath-ejml")) implementation("org.jetbrains.kotlinx:kotlinx-io:0.2.0-npm-dev-11") implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-20") implementation("org.slf4j:slf4j-simple:1.7.30") diff --git a/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt b/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt new file mode 100644 index 000000000..3316f3236 --- /dev/null +++ b/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt @@ -0,0 +1,50 @@ +package kscience.kmath.linear + +import kscience.kmath.commons.linear.CMMatrixContext +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.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(1224) + val dim = 100 + //creating invertible matrix + val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } + val l = Matrix.real(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 } + val matrix = l dot u + + val n = 5000 // iterations + + MatrixContext.real { + repeat(50) { inverse(matrix) } + val inverseTime = measureTimeMillis { repeat(n) { inverse(matrix) } } + println("[kmath] Inversion of $n matrices $dim x $dim finished in $inverseTime millis") + } + + //commons-math + + val commonsTime = measureTimeMillis { + CMMatrixContext { + val cm = matrix.toCM() //avoid overhead on conversion + repeat(n) { inverse(cm) } + } + } + + + println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis") + + val ejmlTime = measureTimeMillis { + (EjmlMatrixContext(RealField)) { + val km = matrix.toEjml() //avoid overhead on conversion + repeat(n) { inverse(km) } + } + } + + println("[ejml] Inversion of $n matrices $dim x $dim finished in $ejmlTime millis") +} diff --git a/examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt b/examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt new file mode 100644 index 000000000..d1011e8f5 --- /dev/null +++ b/examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt @@ -0,0 +1,38 @@ +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") +} diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt index d66530472..f4dbce89a 100644 --- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt +++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt @@ -18,20 +18,52 @@ public interface MatrixContext : SpaceOperations> { */ public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix + public override fun binaryOperation(operation: String, left: Matrix, right: Matrix): Matrix = when (operation) { + "dot" -> left dot right + else -> super.binaryOperation(operation, left, right) + } + + /** + * Computes the dot product of this matrix and another one. + * + * @receiver the multiplicand. + * @param other the multiplier. + * @return the dot product. + */ public infix fun Matrix.dot(other: Matrix): Matrix + /** + * Computes the dot product of this matrix and a vector. + * + * @receiver the multiplicand. + * @param vector the multiplier. + * @return the dot product. + */ public infix fun Matrix.dot(vector: Point): Point + /** + * Multiplies a matrix by its element. + * + * @receiver the multiplicand. + * @param value the multiplier. + * @receiver the product. + */ public operator fun Matrix.times(value: T): Matrix + /** + * Multiplies an element by a matrix of it. + * + * @receiver the multiplicand. + * @param value the multiplier. + * @receiver the product. + */ public operator fun T.times(m: Matrix): Matrix = m * this public companion object { /** * Non-boxing double matrix */ - public val real: RealMatrixContext - get() = RealMatrixContext + public val real: RealMatrixContext = RealMatrixContext /** * A structured matrix with custom buffer @@ -60,7 +92,7 @@ public interface GenericMatrixContext> : MatrixContext { */ public fun point(size: Int, initializer: (Int) -> T): Point - override infix fun Matrix.dot(other: Matrix): Matrix { + public override infix fun Matrix.dot(other: Matrix): Matrix { //TODO add typed error require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" } @@ -71,7 +103,7 @@ public interface GenericMatrixContext> : MatrixContext { } } - override infix fun Matrix.dot(vector: Point): Point { + public override infix fun Matrix.dot(vector: Point): Point { //TODO add typed error require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" } @@ -81,10 +113,10 @@ public interface GenericMatrixContext> : MatrixContext { } } - override operator fun Matrix.unaryMinus(): Matrix = + public override operator fun Matrix.unaryMinus(): Matrix = produce(rowNum, colNum) { i, j -> elementContext { -get(i, j) } } - override fun add(a: Matrix, b: Matrix): Matrix { + public override fun add(a: Matrix, b: Matrix): Matrix { require(a.rowNum == b.rowNum && a.colNum == b.colNum) { "Matrix operation dimension mismatch. [${a.rowNum},${a.colNum}] + [${b.rowNum},${b.colNum}]" } @@ -92,7 +124,7 @@ public interface GenericMatrixContext> : MatrixContext { return produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] + b[i, j] } } } - override operator fun Matrix.minus(b: Matrix): Matrix { + public override operator fun Matrix.minus(b: Matrix): Matrix { require(rowNum == b.rowNum && colNum == b.colNum) { "Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]" } @@ -100,11 +132,11 @@ public interface GenericMatrixContext> : MatrixContext { return produce(rowNum, colNum) { i, j -> elementContext { get(i, j) + b[i, j] } } } - override fun multiply(a: Matrix, k: Number): Matrix = + public override fun multiply(a: Matrix, k: Number): Matrix = produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] * k } } public operator fun Number.times(matrix: FeaturedMatrix): Matrix = matrix * this - override operator fun Matrix.times(value: T): Matrix = + public override operator fun Matrix.times(value: T): Matrix = produce(rowNum, colNum) { i, j -> elementContext { get(i, j) * value } } } diff --git a/kmath-ejml/build.gradle.kts b/kmath-ejml/build.gradle.kts new file mode 100644 index 000000000..fa4aa3e39 --- /dev/null +++ b/kmath-ejml/build.gradle.kts @@ -0,0 +1,8 @@ +plugins { + id("ru.mipt.npm.jvm") +} + +dependencies { + implementation("org.ejml:ejml-simple:0.39") + implementation(project(":kmath-core")) +} diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt new file mode 100644 index 000000000..ed6b1571e --- /dev/null +++ b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt @@ -0,0 +1,71 @@ +package kscience.kmath.ejml + +import org.ejml.dense.row.factory.DecompositionFactory_DDRM +import org.ejml.simple.SimpleMatrix +import kscience.kmath.linear.DeterminantFeature +import kscience.kmath.linear.FeaturedMatrix +import kscience.kmath.linear.LUPDecompositionFeature +import kscience.kmath.linear.MatrixFeature +import kscience.kmath.structures.NDStructure + +/** + * Represents featured matrix over EJML [SimpleMatrix]. + * + * @property origin the underlying [SimpleMatrix]. + * @author Iaroslav Postovalov + */ +public class EjmlMatrix(public val origin: SimpleMatrix, features: Set? = null) : FeaturedMatrix { + public override val rowNum: Int + get() = origin.numRows() + + public override val colNum: Int + get() = origin.numCols() + + public override val shape: IntArray + get() = intArrayOf(origin.numRows(), origin.numCols()) + + public override val features: Set = setOf( + object : LUPDecompositionFeature, DeterminantFeature { + override val determinant: Double + get() = origin.determinant() + + private val lup by lazy { + val ludecompositionF64 = DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols()) + .also { it.decompose(origin.ddrm.copy()) } + + Triple( + EjmlMatrix(SimpleMatrix(ludecompositionF64.getRowPivot(null))), + EjmlMatrix(SimpleMatrix(ludecompositionF64.getLower(null))), + EjmlMatrix(SimpleMatrix(ludecompositionF64.getUpper(null))), + ) + } + + override val l: FeaturedMatrix + get() = lup.second + + override val u: FeaturedMatrix + get() = lup.third + + override val p: FeaturedMatrix + get() = lup.first + } + ) union features.orEmpty() + + public override fun suggestFeature(vararg features: MatrixFeature): EjmlMatrix = + EjmlMatrix(origin, this.features + features) + + public override operator fun get(i: Int, j: Int): Double = origin[i, j] + + public override fun equals(other: Any?): Boolean { + if (other is EjmlMatrix) return origin.isIdentical(other.origin, 0.0) + return NDStructure.equals(this, other as? NDStructure<*> ?: return false) + } + + public override fun hashCode(): Int { + var result = origin.hashCode() + result = 31 * result + features.hashCode() + return result + } + + public override fun toString(): String = "EjmlMatrix(origin=$origin, features=$features)" +} diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt new file mode 100644 index 000000000..52826a7b1 --- /dev/null +++ b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt @@ -0,0 +1,86 @@ +package kscience.kmath.ejml + +import org.ejml.simple.SimpleMatrix +import kscience.kmath.linear.MatrixContext +import kscience.kmath.linear.Point +import kscience.kmath.operations.Space +import kscience.kmath.operations.invoke +import kscience.kmath.structures.Matrix + +/** + * Represents context of basic operations operating with [EjmlMatrix]. + * + * @author Iaroslav Postovalov + */ +public class EjmlMatrixContext(private val space: Space) : MatrixContext { + /** + * Converts this matrix to EJML one. + */ + public fun Matrix.toEjml(): EjmlMatrix = + if (this is EjmlMatrix) this else produce(rowNum, colNum) { i, j -> get(i, j) } + + /** + * Converts this vector to EJML one. + */ + public fun Point.toEjml(): EjmlVector = + if (this is EjmlVector) this else EjmlVector(SimpleMatrix(size, 1).also { + (0 until it.numRows()).forEach { row -> it[row, 0] = get(row) } + }) + + override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): EjmlMatrix = + EjmlMatrix(SimpleMatrix(rows, columns).also { + (0 until it.numRows()).forEach { row -> + (0 until it.numCols()).forEach { col -> it[row, col] = initializer(row, col) } + } + }) + + public override fun Matrix.dot(other: Matrix): EjmlMatrix = + EjmlMatrix(toEjml().origin.mult(other.toEjml().origin)) + + public override fun Matrix.dot(vector: Point): EjmlVector = + EjmlVector(toEjml().origin.mult(vector.toEjml().origin)) + + public override fun add(a: Matrix, b: Matrix): EjmlMatrix = + EjmlMatrix(a.toEjml().origin + b.toEjml().origin) + + public override operator fun Matrix.minus(b: Matrix): EjmlMatrix = + EjmlMatrix(toEjml().origin - b.toEjml().origin) + + public override fun multiply(a: Matrix, k: Number): EjmlMatrix = + produce(a.rowNum, a.colNum) { i, j -> space { a[i, j] * k } } + + public override operator fun Matrix.times(value: Double): EjmlMatrix = EjmlMatrix(toEjml().origin.scale(value)) + + public companion object +} + +/** + * Solves for X in the following equation: x = a^-1*b, where 'a' is base matrix and 'b' is an n by p matrix. + * + * @param a the base matrix. + * @param b n by p matrix. + * @return the solution for 'x' that is n by p. + * @author Iaroslav Postovalov + */ +public fun EjmlMatrixContext.solve(a: Matrix, b: Matrix): EjmlMatrix = + EjmlMatrix(a.toEjml().origin.solve(b.toEjml().origin)) + +/** + * Solves for X in the following equation: x = a^(-1)*b, where 'a' is base matrix and 'b' is an n by p matrix. + * + * @param a the base matrix. + * @param b n by p vector. + * @return the solution for 'x' that is n by p. + * @author Iaroslav Postovalov + */ +public fun EjmlMatrixContext.solve(a: Matrix, b: Point): EjmlVector = + EjmlVector(a.toEjml().origin.solve(b.toEjml().origin)) + +/** + * Returns the inverse of given matrix: b = a^(-1). + * + * @param a the matrix. + * @return the inverse of this matrix. + * @author Iaroslav Postovalov + */ +public fun EjmlMatrixContext.inverse(a: Matrix): EjmlMatrix = EjmlMatrix(a.toEjml().origin.invert()) diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlVector.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlVector.kt new file mode 100644 index 000000000..f7cd1b66d --- /dev/null +++ b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlVector.kt @@ -0,0 +1,40 @@ +package kscience.kmath.ejml + +import org.ejml.simple.SimpleMatrix +import kscience.kmath.linear.Point +import kscience.kmath.structures.Buffer + +/** + * Represents point over EJML [SimpleMatrix]. + * + * @property origin the underlying [SimpleMatrix]. + * @author Iaroslav Postovalov + */ +public class EjmlVector internal constructor(public val origin: SimpleMatrix) : Point { + public override val size: Int + get() = origin.numRows() + + init { + require(origin.numCols() == 1) { "Only single column matrices are allowed" } + } + + public override operator fun get(index: Int): Double = origin[index] + + public override operator fun iterator(): Iterator = object : Iterator { + private var cursor: Int = 0 + + override fun next(): Double { + cursor += 1 + return origin[cursor - 1] + } + + override fun hasNext(): Boolean = cursor < origin.numCols() * origin.numRows() + } + + public override fun contentEquals(other: Buffer<*>): Boolean { + if (other is EjmlVector) return origin.isIdentical(other.origin, 0.0) + return super.contentEquals(other) + } + + public override fun toString(): String = "EjmlVector(origin=$origin)" +} diff --git a/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt b/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt new file mode 100644 index 000000000..e0f15be83 --- /dev/null +++ b/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt @@ -0,0 +1,75 @@ +package kscience.kmath.ejml + +import kscience.kmath.linear.DeterminantFeature +import kscience.kmath.linear.LUPDecompositionFeature +import kscience.kmath.linear.MatrixFeature +import kscience.kmath.linear.getFeature +import org.ejml.dense.row.factory.DecompositionFactory_DDRM +import org.ejml.simple.SimpleMatrix +import kotlin.random.Random +import kotlin.random.asJavaRandom +import kotlin.test.* + +internal class EjmlMatrixTest { + private val random = Random(0) + + private val randomMatrix: SimpleMatrix + get() { + val s = random.nextInt(2, 100) + return SimpleMatrix.random_DDRM(s, s, 0.0, 10.0, random.asJavaRandom()) + } + + @Test + fun rowNum() { + val m = randomMatrix + assertEquals(m.numRows(), EjmlMatrix(m).rowNum) + } + + @Test + fun colNum() { + val m = randomMatrix + assertEquals(m.numCols(), EjmlMatrix(m).rowNum) + } + + @Test + fun shape() { + val m = randomMatrix + val w = EjmlMatrix(m) + assertEquals(listOf(m.numRows(), m.numCols()), w.shape.toList()) + } + + @Test + fun features() { + val m = randomMatrix + val w = EjmlMatrix(m) + val det = w.getFeature>() ?: fail() + assertEquals(m.determinant(), det.determinant) + val lup = w.getFeature>() ?: fail() + + val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows(), m.numCols()) + .also { it.decompose(m.ddrm.copy()) } + + assertEquals(EjmlMatrix(SimpleMatrix(ludecompositionF64.getLower(null))), lup.l) + assertEquals(EjmlMatrix(SimpleMatrix(ludecompositionF64.getUpper(null))), lup.u) + assertEquals(EjmlMatrix(SimpleMatrix(ludecompositionF64.getRowPivot(null))), lup.p) + } + + private object SomeFeature : MatrixFeature {} + + @Test + fun suggestFeature() { + assertNotNull(EjmlMatrix(randomMatrix).suggestFeature(SomeFeature).getFeature()) + } + + @Test + fun get() { + val m = randomMatrix + assertEquals(m[0, 0], EjmlMatrix(m)[0, 0]) + } + + @Test + fun origin() { + val m = randomMatrix + assertSame(m, EjmlMatrix(m).origin) + } +} diff --git a/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlVectorTest.kt b/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlVectorTest.kt new file mode 100644 index 000000000..e27f977d2 --- /dev/null +++ b/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlVectorTest.kt @@ -0,0 +1,47 @@ +package kscience.kmath.ejml + +import org.ejml.simple.SimpleMatrix +import kotlin.random.Random +import kotlin.random.asJavaRandom +import kotlin.test.Test +import kotlin.test.assertEquals +import kotlin.test.assertSame + +internal class EjmlVectorTest { + private val random = Random(0) + + private val randomMatrix: SimpleMatrix + get() = SimpleMatrix.random_DDRM(random.nextInt(2, 100), 1, 0.0, 10.0, random.asJavaRandom()) + + @Test + fun size() { + val m = randomMatrix + val w = EjmlVector(m) + assertEquals(m.numRows(), w.size) + } + + @Test + fun get() { + val m = randomMatrix + val w = EjmlVector(m) + assertEquals(m[0, 0], w[0]) + } + + @Test + fun iterator() { + val m = randomMatrix + val w = EjmlVector(m) + + assertEquals( + m.iterator(true, 0, 0, m.numRows() - 1, 0).asSequence().toList(), + w.iterator().asSequence().toList() + ) + } + + @Test + fun origin() { + val m = randomMatrix + val w = EjmlVector(m) + assertSame(m, w.origin) + } +} diff --git a/settings.gradle.kts b/settings.gradle.kts index 844fa0d6c..7ece3f25c 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -39,5 +39,6 @@ include( ":kmath-for-real", ":kmath-geometry", ":kmath-ast", - ":examples" + ":examples", + ":kmath-ejml" )