Code generation of EJML linear spaces

This commit is contained in:
Iaroslav Postovalov 2021-05-15 02:23:28 +07:00
parent a0621bba5b
commit eb3a8655fb
No known key found for this signature in database
GPG Key ID: 46E15E4A31B3BCD7
9 changed files with 1457 additions and 250 deletions

View File

@ -0,0 +1,5 @@
plugins {
`kotlin-dsl`
}
repositories.mavenCentral()

View File

@ -0,0 +1,415 @@
/*
* 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/LICENSE.txt file.
*/
@file:Suppress("KDocUnresolvedReference")
package space.kscience.kmath.ejml.codegen
import org.intellij.lang.annotations.Language
import java.io.File
private fun Appendable.appendEjmlVector(type: String, ejmlMatrixType: String) {
@Language("kotlin") val text = """/**
* [EjmlVector] specialization for [$type].
*/
public class Ejml${type}Vector<out M : $ejmlMatrixType>(public override val origin: M) : EjmlVector<$type, M>(origin) {
init {
require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" }
}
public override operator fun get(index: Int): $type = origin[0, index]
}"""
appendLine(text)
appendLine()
}
private fun Appendable.appendEjmlMatrix(type: String, ejmlMatrixType: String) {
val text = """/**
* [EjmlMatrix] specialization for [$type].
*/
public class Ejml${type}Matrix<out M : $ejmlMatrixType>(public override val origin: M) : EjmlMatrix<$type, M>(origin) {
public override operator fun get(i: Int, j: Int): $type = origin[i, j]
}"""
appendLine(text)
appendLine()
}
private fun Appendable.appendEjmlLinearSpace(
type: String,
kmathAlgebra: String,
ejmlMatrixParentTypeMatrix: String,
ejmlMatrixType: String,
ejmlMatrixDenseType: String,
ops: String,
denseOps: String,
isDense: Boolean,
) {
@Language("kotlin") val text = """/**
* [EjmlLinearSpace] implementation based on [CommonOps_$ops], [DecompositionFactory_${ops}] operations and
* [${ejmlMatrixType}] matrices.
*/
public object EjmlLinearSpace${ops} : EjmlLinearSpace<${type}, ${kmathAlgebra}, $ejmlMatrixType>() {
/**
* The [${kmathAlgebra}] reference.
*/
public override val elementAlgebra: $kmathAlgebra get() = $kmathAlgebra
@Suppress("UNCHECKED_CAST")
public override fun Matrix<${type}>.toEjml(): Ejml${type}Matrix<${ejmlMatrixType}> = when {
this is Ejml${type}Matrix<*> && origin is $ejmlMatrixType -> this as Ejml${type}Matrix<${ejmlMatrixType}>
else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
}
@Suppress("UNCHECKED_CAST")
public override fun Point<${type}>.toEjml(): Ejml${type}Vector<${ejmlMatrixType}> = when {
this is Ejml${type}Vector<*> && origin is $ejmlMatrixType -> this as Ejml${type}Vector<${ejmlMatrixType}>
else -> Ejml${type}Vector(${ejmlMatrixType}(size, 1).also {
(0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
})
}
public override fun buildMatrix(
rows: Int,
columns: Int,
initializer: ${kmathAlgebra}.(i: Int, j: Int) -> ${type},
): Ejml${type}Matrix<${ejmlMatrixType}> = ${ejmlMatrixType}(rows, columns).also {
(0 until rows).forEach { row ->
(0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
}
}.wrapMatrix()
public override fun buildVector(
size: Int,
initializer: ${kmathAlgebra}.(Int) -> ${type},
): Ejml${type}Vector<${ejmlMatrixType}> = Ejml${type}Vector(${ejmlMatrixType}(size, 1).also {
(0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
})
private fun <T : ${ejmlMatrixParentTypeMatrix}> T.wrapMatrix() = Ejml${type}Matrix(this)
private fun <T : ${ejmlMatrixParentTypeMatrix}> T.wrapVector() = Ejml${type}Vector(this)
public override fun Matrix<${type}>.unaryMinus(): Matrix<${type}> = this * elementAlgebra { -one }
public override fun Matrix<${type}>.dot(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
val out = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.mult(toEjml().origin, other.toEjml().origin, out)
return out.wrapMatrix()
}
public override fun Matrix<${type}>.dot(vector: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
val out = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.mult(toEjml().origin, vector.toEjml().origin, out)
return out.wrapVector()
}
public override operator fun Matrix<${type}>.minus(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
val out = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.add(
elementAlgebra.one,
toEjml().origin,
elementAlgebra { -one },
other.toEjml().origin,
out,${
if (isDense) "" else
"""
null,
null,"""
}
)
return out.wrapMatrix()
}
public override operator fun Matrix<${type}>.times(value: ${type}): Ejml${type}Matrix<${ejmlMatrixType}> {
val res = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.scale(value, toEjml().origin, res)
return res.wrapMatrix()
}
public override fun Point<${type}>.unaryMinus(): Ejml${type}Vector<${ejmlMatrixType}> {
val res = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.changeSign(toEjml().origin, res)
return res.wrapVector()
}
public override fun Matrix<${type}>.plus(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
val out = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.add(
elementAlgebra.one,
toEjml().origin,
elementAlgebra.one,
other.toEjml().origin,
out,${
if (isDense) "" else
"""
null,
null,"""
}
)
return out.wrapMatrix()
}
public override fun Point<${type}>.plus(other: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
val out = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.add(
elementAlgebra.one,
toEjml().origin,
elementAlgebra.one,
other.toEjml().origin,
out,${
if (isDense) "" else
"""
null,
null,"""
}
)
return out.wrapVector()
}
public override fun Point<${type}>.minus(other: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
val out = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.add(
elementAlgebra.one,
toEjml().origin,
elementAlgebra { -one },
other.toEjml().origin,
out,${
if (isDense) "" else
"""
null,
null,"""
}
)
return out.wrapVector()
}
public override fun ${type}.times(m: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> = m * this
public override fun Point<${type}>.times(value: ${type}): Ejml${type}Vector<${ejmlMatrixType}> {
val res = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.scale(value, toEjml().origin, res)
return res.wrapVector()
}
public override fun ${type}.times(v: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> = v * this
@UnstableKMathAPI
public override fun <F : StructureFeature> getFeature(structure: Matrix<${type}>, type: KClass<out F>): F? {
structure.getFeature(type)?.let { return it }
val origin = structure.toEjml().origin
return when (type) {
${
if (isDense)
""" InverseMatrixFeature::class -> object : InverseMatrixFeature<${type}> {
override val inverse: Matrix<${type}> by lazy {
val res = origin.copy()
CommonOps_${ops}.invert(res)
res.wrapMatrix()
}
}
DeterminantFeature::class -> object : DeterminantFeature<${type}> {
override val determinant: $type by lazy { CommonOps_${ops}.det(origin) }
}
SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<${type}> {
private val svd by lazy {
DecompositionFactory_${ops}.svd(origin.numRows, origin.numCols, true, true, false)
.apply { decompose(origin.copy()) }
}
override val u: Matrix<${type}> by lazy { svd.getU(null, false).wrapMatrix() }
override val s: Matrix<${type}> by lazy { svd.getW(null).wrapMatrix() }
override val v: Matrix<${type}> by lazy { svd.getV(null, false).wrapMatrix() }
override val singularValues: Point<${type}> by lazy { ${type}Buffer(svd.singularValues) }
}
QRDecompositionFeature::class -> object : QRDecompositionFeature<${type}> {
private val qr by lazy {
DecompositionFactory_${ops}.qr().apply { decompose(origin.copy()) }
}
override val q: Matrix<${type}> by lazy {
qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
}
override val r: Matrix<${type}> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
}
CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<${type}> {
override val l: Matrix<${type}> by lazy {
val cholesky =
DecompositionFactory_${ops}.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
cholesky.getT(null).wrapMatrix() + LFeature
}
}
LupDecompositionFeature::class -> object : LupDecompositionFeature<${type}> {
private val lup by lazy {
DecompositionFactory_${ops}.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
}
override val l: Matrix<${type}> by lazy {
lup.getLower(null).wrapMatrix() + LFeature
}
override val u: Matrix<${type}> by lazy {
lup.getUpper(null).wrapMatrix() + UFeature
}
override val p: Matrix<${type}> by lazy { lup.getRowPivot(null).wrapMatrix() }
}""" else """ QRDecompositionFeature::class -> object : QRDecompositionFeature<$type> {
private val qr by lazy {
DecompositionFactory_${ops}.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
}
override val q: Matrix<${type}> by lazy {
qr.getQ(null, false).wrapMatrix() + OrthogonalFeature
}
override val r: Matrix<${type}> by lazy { qr.getR(null, false).wrapMatrix() + UFeature }
}
CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<${type}> {
override val l: Matrix<${type}> by lazy {
val cholesky =
DecompositionFactory_${ops}.cholesky().apply { decompose(origin.copy()) }
(cholesky.getT(null) as ${ejmlMatrixParentTypeMatrix}).wrapMatrix() + LFeature
}
}
LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
LUDecompositionFeature<${type}>, DeterminantFeature<${type}>, InverseMatrixFeature<${type}> {
private val lu by lazy {
DecompositionFactory_${ops}.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
}
override val l: Matrix<${type}> by lazy {
lu.getLower(null).wrapMatrix() + LFeature
}
override val u: Matrix<${type}> by lazy {
lu.getUpper(null).wrapMatrix() + UFeature
}
override val inverse: Matrix<${type}> by lazy {
var a = origin
val inverse = ${ejmlMatrixDenseType}(1, 1)
val solver = LinearSolverFactory_${ops}.lu(FillReducing.NONE)
if (solver.modifiesA()) a = a.copy()
val i = CommonOps_${denseOps}.identity(a.numRows)
solver.solve(i, inverse)
inverse.wrapMatrix()
}
override val determinant: $type 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<${type}>, b: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
val res = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.solve(${ejmlMatrixType}(a.toEjml().origin), ${ejmlMatrixType}(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<${type}>, b: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
val res = ${ejmlMatrixType}(1, 1)
CommonOps_${ops}.solve(${ejmlMatrixType}(a.toEjml().origin), ${ejmlMatrixType}(b.toEjml().origin), res)
return Ejml${type}Vector(res)
}
}"""
appendLine(text)
appendLine()
}
/**
* Generates routine EJML classes.
*/
fun ejmlCodegen(outputFile: String): Unit = File(outputFile).run {
parentFile.mkdirs()
writer().use {
it.appendLine("/*")
it.appendLine(" * Copyright 2018-2021 KMath contributors.")
it.appendLine(" * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.")
it.appendLine(" */")
it.appendLine()
it.appendLine("/* This file is generated with buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt */")
it.appendLine()
it.appendLine("package space.kscience.kmath.ejml")
it.appendLine()
it.appendLine("import org.ejml.data.*")
it.appendLine("import space.kscience.kmath.linear.*")
it.appendLine("import space.kscience.kmath.operations.*")
it.appendLine("import space.kscience.kmath.structures.*")
it.appendLine("import space.kscience.kmath.misc.*")
it.appendLine("import kotlin.reflect.*")
it.appendLine("import org.ejml.dense.row.*")
it.appendLine("import org.ejml.dense.row.factory.*")
it.appendLine("import org.ejml.sparse.*")
it.appendLine("import org.ejml.sparse.csc.*")
it.appendLine("import org.ejml.sparse.csc.factory.*")
it.appendLine("import space.kscience.kmath.nd.*")
it.appendLine("import space.kscience.kmath.linear.Matrix")
it.appendLine()
it.appendEjmlVector("Double", "DMatrix")
it.appendEjmlVector("Float", "FMatrix")
it.appendEjmlMatrix("Double", "DMatrix")
it.appendEjmlMatrix("Float", "FMatrix")
it.appendEjmlLinearSpace("Double", "DoubleField", "DMatrix", "DMatrixRMaj", "DMatrixRMaj", "DDRM", "DDRM", true)
it.appendEjmlLinearSpace("Float", "FloatField", "FMatrix", "FMatrixRMaj", "FMatrixRMaj", "FDRM", "FDRM", true)
it.appendEjmlLinearSpace(
type = "Double",
kmathAlgebra = "DoubleField",
ejmlMatrixParentTypeMatrix = "DMatrix",
ejmlMatrixType = "DMatrixSparseCSC",
ejmlMatrixDenseType = "DMatrixRMaj",
ops = "DSCC",
denseOps = "DDRM",
isDense = false,
)
it.appendEjmlLinearSpace(
type = "Float",
kmathAlgebra = "FloatField",
ejmlMatrixParentTypeMatrix = "FMatrix",
ejmlMatrixType = "FMatrixSparseCSC",
ejmlMatrixDenseType = "FMatrixRMaj",
ops = "FSCC",
denseOps = "FDRM",
isDense = false,
)
}
}

View File

@ -75,6 +75,23 @@ public object LFeature : MatrixFeature
*/
public object UFeature : MatrixFeature
/**
* Matrices with this feature support LU factorization: *a = [l] &middot; [u]* where *a* is the owning matrix.
*
* @param T the type of matrices' items.
*/
public interface LUDecompositionFeature<T : Any> : MatrixFeature {
/**
* The lower triangular matrix in this decomposition. It may have [LFeature].
*/
public val l: Matrix<T>
/**
* The upper triangular matrix in this decomposition. It may have [UFeature].
*/
public val u: Matrix<T>
}
/**
* Matrices with this feature support LU factorization with partial pivoting: *[p] &middot; a = [l] &middot; [u]* where
* *a* is the owning matrix.

View File

@ -1,3 +1,5 @@
import space.kscience.kmath.ejml.codegen.ejmlCodegen
plugins {
kotlin("jvm")
id("ru.mipt.npm.gradle.common")
@ -5,6 +7,9 @@ plugins {
dependencies {
api("org.ejml:ejml-ddense:0.40")
api("org.ejml:ejml-fdense:0.40")
api("org.ejml:ejml-dsparse:0.40")
api("org.ejml:ejml-fsparse:0.40")
api(project(":kmath-core"))
}
@ -14,19 +19,24 @@ readme {
feature(
id = "ejml-vector",
description = "Point implementations.",
ref = "src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt"
)
) { "Point implementations." }
feature(
id = "ejml-matrix",
description = "Matrix implementation.",
ref = "src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt"
)
) { "Matrix implementation." }
feature(
id = "ejml-linear-space",
description = "LinearSpace implementations.",
ref = "src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt"
)
) { "LinearSpace implementations." }
}
kotlin.sourceSets.main {
val codegen by tasks.creating {
ejmlCodegen(kotlin.srcDirs.first().absolutePath + "/space/kscience/kmath/ejml/_generated.kt")
}
kotlin.srcDirs(files().builtBy(codegen))
}

View File

@ -5,19 +5,10 @@
package space.kscience.kmath.ejml
import org.ejml.data.DMatrix
import org.ejml.data.DMatrixD1
import org.ejml.data.DMatrixRMaj
import org.ejml.dense.row.CommonOps_DDRM
import org.ejml.dense.row.factory.DecompositionFactory_DDRM
import space.kscience.kmath.linear.*
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.nd.StructureFeature
import space.kscience.kmath.operations.DoubleField
import space.kscience.kmath.linear.LinearSpace
import space.kscience.kmath.linear.Matrix
import space.kscience.kmath.linear.Point
import space.kscience.kmath.operations.Ring
import space.kscience.kmath.structures.DoubleBuffer
import kotlin.reflect.KClass
import kotlin.reflect.cast
/**
* [LinearSpace] implementation specialized for a certain EJML type.
@ -27,7 +18,7 @@ import kotlin.reflect.cast
* @param M the EJML matrix type.
* @author Iaroslav Postovalov
*/
public abstract class EjmlLinearSpace<T : Any, out A : Ring<T>, M : org.ejml.data.Matrix> : LinearSpace<T, A> {
public abstract class EjmlLinearSpace<T : Any, out A : Ring<T>, out M : org.ejml.data.Matrix> : LinearSpace<T, A> {
/**
* Converts this matrix to EJML one.
*/
@ -46,209 +37,3 @@ public abstract class EjmlLinearSpace<T : Any, out A : Ring<T>, M : org.ejml.dat
public abstract override fun buildVector(size: Int, initializer: A.(Int) -> T): EjmlVector<T, M>
}
/**
* [EjmlLinearSpace] implementation based on [CommonOps_DDRM], [DecompositionFactory_DDRM] operations and
* [DMatrixRMaj] matrices.
*
* @author Iaroslav Postovalov
*/
public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, DoubleField, DMatrixRMaj>() {
/**
* The [DoubleField] reference.
*/
public override val elementAlgebra: DoubleField get() = DoubleField
@Suppress("UNCHECKED_CAST")
public 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")
public 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) }
})
}
public override fun buildMatrix(
rows: Int,
columns: Int,
initializer: DoubleField.(i: Int, j: Int) -> Double,
): EjmlDoubleMatrix<DMatrixRMaj> = EjmlDoubleMatrix(DMatrixRMaj(rows, columns).also {
(0 until rows).forEach { row ->
(0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) }
}
})
public 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 : DMatrixD1> T.wrapVector() = EjmlDoubleVector(this)
public override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * (-1.0)
public 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()
}
public 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()
}
public override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
val out = DMatrixRMaj(1, 1)
CommonOps_DDRM.subtract(toEjml().origin, other.toEjml().origin, out)
return out.wrapMatrix()
}
public override operator fun Matrix<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixRMaj> {
val res = this.toEjml().origin.copy()
CommonOps_DDRM.scale(value, res)
return res.wrapMatrix()
}
public override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixRMaj> {
val out = toEjml().origin.copy()
CommonOps_DDRM.changeSign(out)
return out.wrapVector()
}
public override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
val out = DMatrixRMaj(1, 1)
CommonOps_DDRM.add(toEjml().origin, other.toEjml().origin, out)
return out.wrapMatrix()
}
public override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
val out = DMatrixRMaj(1, 1)
CommonOps_DDRM.add(toEjml().origin, other.toEjml().origin, out)
return out.wrapVector()
}
public override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
val out = DMatrixRMaj(1, 1)
CommonOps_DDRM.subtract(toEjml().origin, other.toEjml().origin, out)
return out.wrapVector()
}
public override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> = m * this
public override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixRMaj> {
val res = this.toEjml().origin.copy()
CommonOps_DDRM.scale(value, res)
return res.wrapVector()
}
public override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixRMaj> = v * this
@UnstableKMathAPI
public override fun <F : StructureFeature> getFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
// Return the feature if it is intrinsic to the structure
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)
EjmlDoubleMatrix(res)
}
}
DeterminantFeature::class -> object : DeterminantFeature<Double> {
override val determinant: Double by lazy { CommonOps_DDRM.det(DMatrixRMaj(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 { EjmlDoubleMatrix(svd.getU(null, false)) }
override val s: Matrix<Double> by lazy { EjmlDoubleMatrix(svd.getW(null)) }
override val v: Matrix<Double> by lazy { EjmlDoubleMatrix(svd.getV(null, false)) }
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 {
EjmlDoubleMatrix(qr.getQ(null, false)) + OrthogonalFeature
}
override val r: Matrix<Double> by lazy { EjmlDoubleMatrix(qr.getR(null, false)) + 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()) }
EjmlDoubleMatrix(cholesky.getT(null)) + 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 {
EjmlDoubleMatrix(lup.getLower(null)) + LFeature
}
override val u: Matrix<Double> by lazy {
EjmlDoubleMatrix(lup.getUpper(null)) + UFeature
}
override val p: Matrix<Double> by lazy { EjmlDoubleMatrix(lup.getRowPivot(null)) }
}
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 EjmlDoubleMatrix(res)
}
/**
* 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)
}
}

View File

@ -5,7 +5,6 @@
package space.kscience.kmath.ejml
import org.ejml.data.DMatrix
import org.ejml.data.Matrix
import space.kscience.kmath.nd.Structure2D
@ -21,12 +20,3 @@ public abstract class EjmlMatrix<T, out M : Matrix>(public open val origin: M) :
public override val rowNum: Int get() = origin.numRows
public override val colNum: Int get() = origin.numCols
}
/**
* [EjmlMatrix] specialization for [Double].
*
* @author Iaroslav Postovalov
*/
public class EjmlDoubleMatrix<out M : DMatrix>(public override val origin: M) : EjmlMatrix<Double, M>(origin) {
public override operator fun get(i: Int, j: Int): Double = origin[i, j]
}

View File

@ -5,7 +5,6 @@
package space.kscience.kmath.ejml
import org.ejml.data.DMatrixD1
import org.ejml.data.Matrix
import space.kscience.kmath.linear.Point
@ -14,12 +13,12 @@ import space.kscience.kmath.linear.Point
*
* @param T the type of elements contained in the buffer.
* @param M the type of EJML matrix.
* @property origin The underlying matrix.
* @property origin The underlying matrix, must have only one row.
* @author Iaroslav Postovalov
*/
public abstract class EjmlVector<out T, out M : Matrix>(public open val origin: M) : Point<T> {
public override val size: Int
get() = origin.numRows
get() = origin.numCols
public override operator fun iterator(): Iterator<T> = object : Iterator<T> {
private var cursor: Int = 0
@ -34,12 +33,3 @@ public abstract class EjmlVector<out T, out M : Matrix>(public open val origin:
public override fun toString(): String = "EjmlVector(origin=$origin)"
}
/**
* [EjmlVector] specialization for [Double].
*
* @author Iaroslav Postovalov
*/
public class EjmlDoubleVector<out M : DMatrixD1>(public override val origin: M) : EjmlVector<Double, M>(origin) {
public override operator fun get(index: Int): Double = origin[index]
}

View File

@ -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/LICENSE.txt 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>(public 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" }
}
public override operator fun get(index: Int): Double = origin[0, index]
}
/**
* [EjmlVector] specialization for [Float].
*/
public class EjmlFloatVector<out M : FMatrix>(public 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" }
}
public override operator fun get(index: Int): Float = origin[0, index]
}
/**
* [EjmlMatrix] specialization for [Double].
*/
public class EjmlDoubleMatrix<out M : DMatrix>(public override val origin: M) : EjmlMatrix<Double, M>(origin) {
public override operator fun get(i: Int, j: Int): Double = origin[i, j]
}
/**
* [EjmlMatrix] specialization for [Float].
*/
public class EjmlFloatMatrix<out M : FMatrix>(public override val origin: M) : EjmlMatrix<Float, M>(origin) {
public 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.
*/
public override val elementAlgebra: DoubleField get() = DoubleField
@Suppress("UNCHECKED_CAST")
public 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")
public 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) }
})
}
public 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()
public 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)
public override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
public 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()
}
public 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()
}
public 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()
}
public 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()
}
public override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixRMaj> {
val res = DMatrixRMaj(1, 1)
CommonOps_DDRM.changeSign(toEjml().origin, res)
return res.wrapVector()
}
public 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()
}
public 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()
}
public 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()
}
public override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> = m * this
public override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixRMaj> {
val res = DMatrixRMaj(1, 1)
CommonOps_DDRM.scale(value, toEjml().origin, res)
return res.wrapVector()
}
public override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixRMaj> = v * this
@UnstableKMathAPI
public override fun <F : StructureFeature> getFeature(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() + OrthogonalFeature
}
override val r: Matrix<Double> by lazy { qr.getR(null, false).wrapMatrix() + 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() + 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() + LFeature
}
override val u: Matrix<Double> by lazy {
lup.getUpper(null).wrapMatrix() + 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.
*/
public override val elementAlgebra: FloatField get() = FloatField
@Suppress("UNCHECKED_CAST")
public 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")
public 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) }
})
}
public 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()
public 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)
public override fun Matrix<Float>.unaryMinus(): Matrix<Float> = this * elementAlgebra { -one }
public 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()
}
public 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()
}
public 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()
}
public 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()
}
public override fun Point<Float>.unaryMinus(): EjmlFloatVector<FMatrixRMaj> {
val res = FMatrixRMaj(1, 1)
CommonOps_FDRM.changeSign(toEjml().origin, res)
return res.wrapVector()
}
public 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()
}
public 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()
}
public 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()
}
public override fun Float.times(m: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> = m * this
public override fun Point<Float>.times(value: Float): EjmlFloatVector<FMatrixRMaj> {
val res = FMatrixRMaj(1, 1)
CommonOps_FDRM.scale(value, toEjml().origin, res)
return res.wrapVector()
}
public override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixRMaj> = v * this
@UnstableKMathAPI
public override fun <F : StructureFeature> getFeature(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() + OrthogonalFeature
}
override val r: Matrix<Float> by lazy { qr.getR(null, false).wrapMatrix() + 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() + 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() + LFeature
}
override val u: Matrix<Float> by lazy {
lup.getUpper(null).wrapMatrix() + 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.
*/
public override val elementAlgebra: DoubleField get() = DoubleField
@Suppress("UNCHECKED_CAST")
public 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")
public 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) }
})
}
public 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()
public 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)
public override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
public 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()
}
public 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()
}
public 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()
}
public 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()
}
public override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixSparseCSC> {
val res = DMatrixSparseCSC(1, 1)
CommonOps_DSCC.changeSign(toEjml().origin, res)
return res.wrapVector()
}
public 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()
}
public 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()
}
public 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()
}
public override fun Double.times(m: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> = m * this
public override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixSparseCSC> {
val res = DMatrixSparseCSC(1, 1)
CommonOps_DSCC.scale(value, toEjml().origin, res)
return res.wrapVector()
}
public override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> = v * this
@UnstableKMathAPI
public override fun <F : StructureFeature> getFeature(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() + OrthogonalFeature
}
override val r: Matrix<Double> by lazy { qr.getR(null, false).wrapMatrix() + 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() + 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() + LFeature
}
override val u: Matrix<Double> by lazy {
lu.getUpper(null).wrapMatrix() + 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.
*/
public override val elementAlgebra: FloatField get() = FloatField
@Suppress("UNCHECKED_CAST")
public 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")
public 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) }
})
}
public 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()
public 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)
public override fun Matrix<Float>.unaryMinus(): Matrix<Float> = this * elementAlgebra { -one }
public 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()
}
public 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()
}
public 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()
}
public 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()
}
public override fun Point<Float>.unaryMinus(): EjmlFloatVector<FMatrixSparseCSC> {
val res = FMatrixSparseCSC(1, 1)
CommonOps_FSCC.changeSign(toEjml().origin, res)
return res.wrapVector()
}
public 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()
}
public 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()
}
public 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()
}
public override fun Float.times(m: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> = m * this
public override fun Point<Float>.times(value: Float): EjmlFloatVector<FMatrixSparseCSC> {
val res = FMatrixSparseCSC(1, 1)
CommonOps_FSCC.scale(value, toEjml().origin, res)
return res.wrapVector()
}
public override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> = v * this
@UnstableKMathAPI
public override fun <F : StructureFeature> getFeature(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() + OrthogonalFeature
}
override val r: Matrix<Float> by lazy { qr.getR(null, false).wrapMatrix() + 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() + 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() + LFeature
}
override val u: Matrix<Float> by lazy {
lu.getUpper(null).wrapMatrix() + 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)
}
}

View File

@ -18,7 +18,7 @@ internal class EjmlVectorTest {
private val randomMatrix: DMatrixRMaj
get() {
val d = DMatrixRMaj(random.nextInt(2, 100), 1)
val d = DMatrixRMaj(1, random.nextInt(2, 100))
RandomMatrices_DDRM.fillUniform(d, random.asJavaRandom())
return d
}
@ -27,7 +27,7 @@ internal class EjmlVectorTest {
fun size() {
val m = randomMatrix
val w = EjmlDoubleVector(m)
assertEquals(m.numRows, w.size)
assertEquals(m.numCols, w.size)
}
@Test
@ -43,7 +43,7 @@ internal class EjmlVectorTest {
val w = EjmlDoubleVector(m)
assertEquals(
m.iterator(true, 0, 0, m.numRows - 1, 0).asSequence().toList(),
m.iterator(true, 0, 0, 0, m.numCols - 1).asSequence().toList(),
w.iterator().asSequence().toList()
)
}