kscience/kmath
9
0
Fork 1
Code Issues Pull Requests 6 Packages Projects Releases 16 Wiki Activity

Deprecated Vectors. Working on LUP optimization (not working yet)

Browse Source
This commit is contained in:
Alexander Nozik 2019-04-15 20:50:50 +03:00
parent 98bb72a6a0
commit f1b1010c4d
18 changed files with 347 additions and 247 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
benchmarks/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt
View File

@ -1,10 +1,13 @@
package scientifik.kmath.linear
import koma.matrix.ejml.EJMLMatrixFactory
import scientifik.kmath.operations.RealField
import scientifik.kmath.structures.Matrix
import kotlin.contracts.ExperimentalContracts
import kotlin.random.Random
import kotlin.system.measureTimeMillis
@ExperimentalContracts
fun main() {
val random = Random(12224)
val dim = 100
@ -32,10 +35,8 @@ fun main() {
//commons-math
val cmContext = CMLUPSolver
val commonsTime = measureTimeMillis {
cmContext.run {
CMMatrixContext.run {
val cm = matrix.toCM() //avoid overhead on conversion
repeat(n) {
val res = inverse(cm)
@ -48,10 +49,8 @@ fun main() {
//koma-ejml
val komaContext = KomaMatrixContext(EJMLMatrixFactory())
val komaTime = measureTimeMillis {
komaContext.run {
KomaMatrixContext(EJMLMatrixFactory(), RealField).run {
val km = matrix.toKoma() //avoid overhead on conversion
repeat(n) {
val res = inverse(km)

3
benchmarks/src/main/kotlin/scientifik/kmath/linear/MultiplicationBenchmark.kt
View File

@ -1,6 +1,7 @@
package scientifik.kmath.linear
import koma.matrix.ejml.EJMLMatrixFactory
import scientifik.kmath.operations.RealField
import scientifik.kmath.structures.Matrix
import kotlin.random.Random
import kotlin.system.measureTimeMillis
@ -27,7 +28,7 @@ fun main() {
}
KomaMatrixContext(EJMLMatrixFactory()).run {
KomaMatrixContext(EJMLMatrixFactory(), RealField).run {
val komaMatrix1 = matrix1.toKoma()
val komaMatrix2 = matrix2.toKoma()

2
benchmarks/src/main/kotlin/scientifik/kmath/structures/ComplexNDBenchmark.kt
View File

@ -23,7 +23,7 @@ fun main() {
val complexTime = measureTimeMillis {
complexField.run {
var res: ComplexNDElement = one
var res = one
repeat(n) {
res += 1.0
}

5
build.gradle.kts
View File

@ -1,9 +1,6 @@
import com.moowork.gradle.node.NodeExtension
import com.moowork.gradle.node.npm.NpmTask
import com.moowork.gradle.node.task.NodeTask
import org.jetbrains.kotlin.gradle.dsl.KotlinMultiplatformExtension
import org.jetbrains.kotlin.gradle.tasks.Kotlin2JsCompile
import org.jetbrains.kotlin.gradle.tasks.KotlinCompile
buildscript {
val kotlinVersion: String by rootProject.extra("1.3.21")
@ -194,6 +191,8 @@ subprojects {
sourceSets.all {
languageSettings.progressiveMode = true
languageSettings.enableLanguageFeature("InlineClasses")
languageSettings.useExperimentalAnnotation("ExperimentalContracts")
//languageSettings.enableLanguageFeature("Contracts")
}
}

20
doc/linear.md
View File

@ -1 +1,19 @@
**TODO**
## Basic linear algebra layout
Kmath support for linear algebra organized in a context-oriented way. Meaning that operations are in most cases declared
in context classes, and are not the members of classes that store data. This allows more flexible approach to maintain multiple
back-ends. The new operations added as extensions to contexts instead of being member functions of data structures.
Two major contexts used for linear algebra and hyper-geometry:
* `VectorSpace` forms a mathematical space on top of array-like structure (`Buffer` and its typealias `Point` used for geometry).
* `MatrixContext` forms a space-like context for 2d-structures. It does not store matrix size and therefore does not implement
`Space` interface (it is not possible to create zero element without knowing the matrix size).
## Vector spaces
## Matrix operations
## Back-end overview

29
kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMMatrix.kt
View File

@ -27,7 +27,7 @@ fun Matrix<Double>.toCM(): CMMatrix = if (this is CMMatrix) {
CMMatrix(Array2DRowRealMatrix(array))
}
fun RealMatrix.toMatrix() = CMMatrix(this)
fun RealMatrix.asMatrix() = CMMatrix(this)
class CMVector(val origin: RealVector) : Point<Double> {
override val size: Int get() = origin.dimension
@ -47,7 +47,6 @@ fun Point<Double>.toCM(): CMVector = if (this is CMVector) {
fun RealVector.toPoint() = CMVector(this)
object CMMatrixContext : MatrixContext<Double> {
override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): CMMatrix {
val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } }
return CMMatrix(Array2DRowRealMatrix(array))
@ -59,35 +58,19 @@ object CMMatrixContext : MatrixContext<Double> {
override fun Matrix<Double>.dot(vector: Point<Double>): CMVector =
CMVector(this.toCM().origin.preMultiply(vector.toCM().origin))
override fun Matrix<Double>.unaryMinus(): CMMatrix =
produce(rowNum, colNum) { i, j -> -get(i, j) }
override fun Matrix<Double>.plus(b: Matrix<Double>) =
CMMatrix(this.toCM().origin.multiply(b.toCM().origin))
override fun add(a: Matrix<Double>, b: Matrix<Double>) =
CMMatrix(a.toCM().origin.multiply(b.toCM().origin))
override fun Matrix<Double>.minus(b: Matrix<Double>) =
CMMatrix(this.toCM().origin.subtract(b.toCM().origin))
override fun Matrix<Double>.times(value: Double) =
CMMatrix(this.toCM().origin.scalarMultiply(value.toDouble()))
}
override fun multiply(a: Matrix<Double>, k: Number) =
CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble()))
object CMLUPSolver: LinearSolver<Double>{
override fun solve(a: Matrix<Double>, b: Matrix<Double>): CMMatrix {
val decomposition = LUDecomposition(a.toCM().origin)
return decomposition.solver.solve(b.toCM().origin).toMatrix()
}
override fun solve(a: Matrix<Double>, b: Point<Double>): CMVector {
val decomposition = LUDecomposition(a.toCM().origin)
return decomposition.solver.solve(b.toCM().origin).toPoint()
}
override fun inverse(a: Matrix<Double>): CMMatrix {
val decomposition = LUDecomposition(a.toCM().origin)
return decomposition.solver.inverse.toMatrix()
}
override fun Matrix<Double>.times(value: Double): Matrix<Double> = produce(rowNum,colNum){i,j-> get(i,j)*value}
}
operator fun CMMatrix.plus(other: CMMatrix): CMMatrix = CMMatrix(this.origin.add(other.origin))

39
kmath-commons/src/main/kotlin/scientifik/kmath/linear/CMSolver.kt Normal file
View File

@ -0,0 +1,39 @@
package scientifik.kmath.linear
import org.apache.commons.math3.linear.*
import scientifik.kmath.structures.Matrix
enum class CMDecomposition {
LUP,
QR,
RRQR,
EIGEN,
CHOLESKY
}
fun CMMatrixContext.solver(a: Matrix<Double>, decomposition: CMDecomposition = CMDecomposition.LUP) =
when (decomposition) {
CMDecomposition.LUP -> LUDecomposition(a.toCM().origin).solver
CMDecomposition.RRQR -> RRQRDecomposition(a.toCM().origin).solver
CMDecomposition.QR -> QRDecomposition(a.toCM().origin).solver
CMDecomposition.EIGEN -> EigenDecomposition(a.toCM().origin).solver
CMDecomposition.CHOLESKY -> CholeskyDecomposition(a.toCM().origin).solver
}
fun CMMatrixContext.solve(
a: Matrix<Double>,
b: Matrix<Double>,
decomposition: CMDecomposition = CMDecomposition.LUP
) = solver(a, decomposition).solve(b.toCM().origin).asMatrix()
fun CMMatrixContext.solve(
a: Matrix<Double>,
b: Point<Double>,
decomposition: CMDecomposition = CMDecomposition.LUP
) = solver(a, decomposition).solve(b.toCM().origin).toPoint()
fun CMMatrixContext.inverse(
a: Matrix<Double>,
decomposition: CMDecomposition = CMDecomposition.LUP
) = solver(a, decomposition).inverse.asMatrix()

1
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt
View File

@ -1,6 +1,5 @@
package scientifik.kmath.linear
import scientifik.kmath.operations.RealField
import scientifik.kmath.operations.Ring
import scientifik.kmath.structures.*

103
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/FeaturedMatrix.kt
View File

@ -7,109 +7,6 @@ import scientifik.kmath.structures.*
import scientifik.kmath.structures.Buffer.Companion.boxing
import kotlin.math.sqrt
/**
* Basic operations on matrices. Operates on [Matrix]
*/
interface MatrixContext<T : Any> {
/**
* Produce a matrix with this context and given dimensions
*/
fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T>
infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T>
infix fun Matrix<T>.dot(vector: Point<T>): Point<T>
operator fun Matrix<T>.unaryMinus(): Matrix<T>
operator fun Matrix<T>.plus(b: Matrix<T>): Matrix<T>
operator fun Matrix<T>.minus(b: Matrix<T>): Matrix<T>
operator fun Matrix<T>.times(value: T): Matrix<T>
operator fun T.times(m: Matrix<T>): Matrix<T> = m * this
companion object {
/**
* Non-boxing double matrix
*/
val real = BufferMatrixContext(RealField, Buffer.Companion::auto)
/**
* A structured matrix with custom buffer
*/
fun <T : Any, R : Ring<T>> buffered(
ring: R,
bufferFactory: BufferFactory<T> = ::boxing
): GenericMatrixContext<T, R> =
BufferMatrixContext(ring, bufferFactory)
/**
* Automatic buffered matrix, unboxed if it is possible
*/
inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R> =
buffered(ring, Buffer.Companion::auto)
}
}
interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
/**
* The ring context for matrix elements
*/
val elementContext: R
/**
* Produce a point compatible with matrix space
*/
fun point(size: Int, initializer: (Int) -> T): Point<T>
override infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T> {
//TODO add typed error
if (this.colNum != other.rowNum) error("Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})")
return produce(rowNum, other.colNum) { i, j ->
val row = rows[i]
val column = other.columns[j]
with(elementContext) {
sum(row.asSequence().zip(column.asSequence(), ::multiply))
}
}
}
override infix fun Matrix<T>.dot(vector: Point<T>): Point<T> {
//TODO add typed error
if (this.colNum != vector.size) error("Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})")
return point(rowNum) { i ->
val row = rows[i]
with(elementContext) {
sum(row.asSequence().zip(vector.asSequence(), ::multiply))
}
}
}
override operator fun Matrix<T>.unaryMinus() =
produce(rowNum, colNum) { i, j -> elementContext.run { -get(i, j) } }
override operator fun Matrix<T>.plus(b: Matrix<T>): Matrix<T> {
if (rowNum != b.rowNum || colNum != b.colNum) error("Matrix operation dimension mismatch. [$rowNum,$colNum] + [${b.rowNum},${b.colNum}]")
return produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) + b[i, j] } }
}
override operator fun Matrix<T>.minus(b: Matrix<T>): Matrix<T> {
if (rowNum != b.rowNum || colNum != b.colNum) error("Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]")
return produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) + b[i, j] } }
}
operator fun Matrix<T>.times(number: Number): Matrix<T> =
produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * number } }
operator fun Number.times(matrix: FeaturedMatrix<T>): Matrix<T> = matrix * this
override fun Matrix<T>.times(value: T): Matrix<T> =
produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * value } }
}
/**
* A 2d structure plus optional matrix-specific features
*/

67
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUPDecomposition.kt
View File

@ -4,6 +4,9 @@ import scientifik.kmath.operations.Field
import scientifik.kmath.operations.RealField
import scientifik.kmath.operations.Ring
import scientifik.kmath.structures.*
import kotlin.contracts.ExperimentalContracts
import kotlin.contracts.InvocationKind
import kotlin.contracts.contract
import kotlin.reflect.KClass
/**
@ -63,7 +66,12 @@ class LUPDecomposition<T : Any>(
}
open class BufferAccessor<T : Any>(val type: KClass<T>, val field: Field<T>, val rowNum: Int, val colNum: Int) {
internal open class BufferAccessor<T : Any>(
val type: KClass<T>,
val field: Field<T>,
val rowNum: Int,
val colNum: Int
) {
open operator fun MutableBuffer<T>.get(i: Int, j: Int) = get(i + colNum * j)
open operator fun MutableBuffer<T>.set(i: Int, j: Int, value: T) {
set(i + colNum * j, value)
@ -102,7 +110,8 @@ open class BufferAccessor<T : Any>(val type: KClass<T>, val field: Field<T>, val
/**
* Specialized LU operations for Doubles
*/
class RealBufferAccessor(rowNum: Int, colNum: Int) : BufferAccessor<Double>(Double::class, RealField, rowNum, colNum) {
private class RealBufferAccessor(rowNum: Int, colNum: Int) :
BufferAccessor<Double>(Double::class, RealField, rowNum, colNum) {
override fun MutableBuffer<Double>.get(i: Int, j: Int) = (this as DoubleBuffer).array[i + colNum * j]
override fun MutableBuffer<Double>.set(i: Int, j: Int, value: Double) {
(this as DoubleBuffer).array[i + colNum * j] = value
@ -125,24 +134,33 @@ class RealBufferAccessor(rowNum: Int, colNum: Int) : BufferAccessor<Double>(Doub
}
}
fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.buildAccessor(
type:KClass<T>,
@ExperimentalContracts
private inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.withAccessor(
type: KClass<T>,
rowNum: Int,
colNum: Int
): BufferAccessor<T> {
return if (elementContext == RealField) {
colNum: Int,
block: BufferAccessor<T>.() -> Unit
) {
contract {
callsInPlace(block, InvocationKind.EXACTLY_ONCE)
}
if (elementContext == RealField) {
@Suppress("UNCHECKED_CAST")
RealBufferAccessor(rowNum, colNum) as BufferAccessor<T>
} else {
BufferAccessor(type, elementContext, rowNum, colNum)
}
}.run(block)
}
fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.abs(value: T) =
private fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.abs(value: T) =
if (value > elementContext.zero) value else with(elementContext) { -value }
fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lupDecompose(
/**
* Create a lup decomposition of generic matrix
*/
@ExperimentalContracts
fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
type: KClass<T>,
matrix: Matrix<T>,
checkSingular: (T) -> Boolean
@ -155,7 +173,7 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lupDecompose(
val m = matrix.colNum
val pivot = IntArray(matrix.rowNum)
buildAccessor(type, matrix.rowNum, matrix.colNum).run {
withAccessor(type, matrix.rowNum, matrix.colNum) {
val lu = create(matrix)
@ -170,21 +188,12 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lupDecompose(
// upper
for (row in 0 until col) {
// var sum = lu[row, col]
// for (i in 0 until row) {
// sum -= lu[row, i] * lu[i, col]
// }
val sum = lu.innerProduct(row, col, row)
lu[row, col] = field.run { lu[row, col] - sum }
}
// lower
val max = (col until m).maxBy { row ->
// var sum = lu[row, col]
// for (i in 0 until col) {
// sum -= lu[row, i] * lu[i, col]
// }
// lu[row, col] = sum
val sum = lu.innerProduct(row, col, col)
lu[row, col] = field.run { lu[row, col] - sum }
abs(sum)
@ -214,14 +223,17 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lupDecompose(
//lu[row, col] = lu[row, col] / luDiag
}
}
return scientifik.kmath.linear.LUPDecomposition(elementContext, lu.collect(), pivot, even)
return LUPDecomposition(elementContext, lu.collect(), pivot, even)
}
}
@ExperimentalContracts
fun GenericMatrixContext<Double, RealField>.lup(matrix: Matrix<Double>) = lup(Double::class, matrix) { it < 1e-11 }
/**
* Solve a linear equation **a*x = b**
*/
@ExperimentalContracts
fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.solve(
type: KClass<T>,
a: Matrix<T>,
@ -233,9 +245,9 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.solve(
}
// Use existing decomposition if it is provided by matrix
val decomposition = a.getFeature() ?: lupDecompose(type, a, checkSingular)
val decomposition = a.getFeature() ?: lup(type, a, checkSingular)
buildAccessor(type, a.rowNum, a.colNum).run {
withAccessor(type, a.rowNum, a.colNum) {
val lu = create(decomposition.lu)
@ -271,14 +283,19 @@ fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.solve(
return produce(a.rowNum, a.colNum) { i, j -> bp[i, j] }
}
}
@ExperimentalContracts
fun GenericMatrixContext<Double, RealField>.solve(a: Matrix<Double>, b: Matrix<Double>) =
solve(Double::class, a, b) { it < 1e-11 }
@ExperimentalContracts
inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.inverse(
matrix: Matrix<T>,
noinline checkSingular: (T) -> Boolean
) =
solve(T::class, matrix, one(matrix.rowNum, matrix.colNum), checkSingular)
@ExperimentalContracts
fun GenericMatrixContext<Double, RealField>.inverse(matrix: Matrix<Double>) =
inverse(matrix) { it < 1e-11 }

4
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
View File

@ -13,7 +13,7 @@ import scientifik.kmath.structures.asSequence
*/
interface LinearSolver<T : Any> {
fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T>
fun solve(a: Matrix<T>, b: Point<T>): Point<T> = solve(a, b.toMatrix()).asPoint()
fun solve(a: Matrix<T>, b: Point<T>): Point<T> = solve(a, b.asMatrix()).asPoint()
fun inverse(a: Matrix<T>): Matrix<T>
}
@ -43,4 +43,4 @@ fun <T : Any> Matrix<T>.asPoint(): Point<T> =
error("Can't convert matrix with more than one column to vector")
}
fun <T : Any> Point<T>.toMatrix() = VirtualMatrix(size, 1) { i, _ -> get(i) }
fun <T : Any> Point<T>.asMatrix() = VirtualMatrix(size, 1) { i, _ -> get(i) }

106
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/MatrixContext.kt Normal file
View File

@ -0,0 +1,106 @@
package scientifik.kmath.linear
import scientifik.kmath.operations.RealField
import scientifik.kmath.operations.Ring
import scientifik.kmath.operations.SpaceOperations
import scientifik.kmath.operations.sum
import scientifik.kmath.structures.Buffer
import scientifik.kmath.structures.BufferFactory
import scientifik.kmath.structures.Matrix
import scientifik.kmath.structures.asSequence
/**
* Basic operations on matrices. Operates on [Matrix]
*/
interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
/**
* Produce a matrix with this context and given dimensions
*/
fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T>
infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T>
infix fun Matrix<T>.dot(vector: Point<T>): Point<T>
operator fun Matrix<T>.times(value: T): Matrix<T>
operator fun T.times(m: Matrix<T>): Matrix<T> = m * this
companion object {
/**
* Non-boxing double matrix
*/
val real = BufferMatrixContext(RealField, Buffer.Companion::auto)
/**
* A structured matrix with custom buffer
*/
fun <T : Any, R : Ring<T>> buffered(
ring: R,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
): GenericMatrixContext<T, R> =
BufferMatrixContext(ring, bufferFactory)
/**
* Automatic buffered matrix, unboxed if it is possible
*/
inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R> =
buffered(ring, Buffer.Companion::auto)
}
}
interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
/**
* The ring context for matrix elements
*/
val elementContext: R
/**
* Produce a point compatible with matrix space
*/
fun point(size: Int, initializer: (Int) -> T): Point<T>
override infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T> {
//TODO add typed error
if (this.colNum != other.rowNum) error("Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})")
return produce(rowNum, other.colNum) { i, j ->
val row = rows[i]
val column = other.columns[j]
with(elementContext) {
sum(row.asSequence().zip(column.asSequence(), ::multiply))
}
}
}
override infix fun Matrix<T>.dot(vector: Point<T>): Point<T> {
//TODO add typed error
if (this.colNum != vector.size) error("Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})")
return point(rowNum) { i ->
val row = rows[i]
with(elementContext) {
sum(row.asSequence().zip(vector.asSequence(), ::multiply))
}
}
}
override operator fun Matrix<T>.unaryMinus() =
produce(rowNum, colNum) { i, j -> elementContext.run { -get(i, j) } }
override fun add(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
if (a.rowNum != b.rowNum || a.colNum != b.colNum) error("Matrix operation dimension mismatch. [${a.rowNum},${a.colNum}] + [${b.rowNum},${b.colNum}]")
return produce(a.rowNum, a.colNum) { i, j -> elementContext.run { a.get(i, j) + b[i, j] } }
}
override operator fun Matrix<T>.minus(b: Matrix<T>): Matrix<T> {
if (rowNum != b.rowNum || colNum != b.colNum) error("Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]")
return produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) + b[i, j] } }
}
override fun multiply(a: Matrix<T>, k: Number): Matrix<T> =
produce(a.rowNum, a.colNum) { i, j -> elementContext.run { a.get(i, j) * k } }
operator fun Number.times(matrix: FeaturedMatrix<T>): Matrix<T> = matrix * this
override fun Matrix<T>.times(value: T): Matrix<T> =
produce(rowNum, colNum) { i, j -> elementContext.run { get(i, j) * value } }
}

70
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
View File

@ -4,68 +4,23 @@ import scientifik.kmath.operations.RealField
import scientifik.kmath.operations.Space
import scientifik.kmath.operations.SpaceElement
import scientifik.kmath.structures.Buffer
import scientifik.kmath.structures.BufferFactory
import scientifik.kmath.structures.asSequence
import kotlin.jvm.JvmName
typealias Point<T> = Buffer<T>
/**
* A linear space for vectors.
* Could be used on any point-like structure
*/
interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
val size: Int
val space: S
fun produce(initializer: (Int) -> T): Point<T>
/**
* Produce a space-element of this vector space for expressions
*/
fun produceElement(initializer: (Int) -> T): Vector<T, S>
override val zero: Point<T> get() = produce { space.zero }
override fun add(a: Point<T>, b: Point<T>): Point<T> = produce { with(space) { a[it] + b[it] } }
override fun multiply(a: Point<T>, k: Number): Point<T> = produce { with(space) { a[it] * k } }
//TODO add basis
companion object {
private val realSpaceCache = HashMap<Int, BufferVectorSpace<Double, RealField>>()
/**
* Non-boxing double vector space
*/
fun real(size: Int): BufferVectorSpace<Double, RealField> {
return realSpaceCache.getOrPut(size) { BufferVectorSpace(size, RealField, Buffer.Companion::auto) }
}
/**
* A structured vector space with custom buffer
*/
fun <T : Any, S : Space<T>> buffered(
size: Int,
space: S,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
): VectorSpace<T, S> = BufferVectorSpace(size, space, bufferFactory)
/**
* Automatic buffered vector, unboxed if it is possible
*/
inline fun <reified T : Any, S : Space<T>> smart(size: Int, space: S): VectorSpace<T, S> =
buffered(size, space, Buffer.Companion::auto)
}
}
fun <T : Any, S : Space<T>> BufferVectorSpace<T, S>.produceElement(initializer: (Int) -> T): Vector<T, S> =
BufferVector(this, produce(initializer))
@JvmName("produceRealElement")
fun BufferVectorSpace<Double, RealField>.produceElement(initializer: (Int) -> Double): Vector<Double, RealField> =
BufferVector(this, produce(initializer))
/**
* A point coupled to the linear space
*/
@Deprecated("Use VectorContext instead")
interface Vector<T : Any, S : Space<T>> : SpaceElement<Point<T>, Vector<T, S>, VectorSpace<T, S>>, Point<T> {
override val size: Int get() = context.size
@ -90,16 +45,7 @@ interface Vector<T : Any, S : Space<T>> : SpaceElement<Point<T>, Vector<T, S>, V
}
}
data class BufferVectorSpace<T : Any, S : Space<T>>(
override val size: Int,
override val space: S,
val bufferFactory: BufferFactory<T>
) : VectorSpace<T, S> {
override fun produce(initializer: (Int) -> T) = bufferFactory(size, initializer)
override fun produceElement(initializer: (Int) -> T): Vector<T, S> = BufferVector(this, produce(initializer))
}
@Deprecated("Use VectorContext instead")
data class BufferVector<T : Any, S : Space<T>>(override val context: VectorSpace<T, S>, val buffer: Buffer<T>) :
Vector<T, S> {

75
kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/VectorSpace.kt Normal file
View File

@ -0,0 +1,75 @@
package scientifik.kmath.linear
import scientifik.kmath.operations.RealField
import scientifik.kmath.operations.Space
import scientifik.kmath.structures.Buffer
import scientifik.kmath.structures.BufferFactory
/**
* A linear space for vectors.
* Could be used on any point-like structure
*/
interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
val size: Int
val space: S
fun produce(initializer: (Int) -> T): Point<T>
/**
* Produce a space-element of this vector space for expressions
*/
//fun produceElement(initializer: (Int) -> T): Vector<T, S>
override val zero: Point<T> get() = produce { space.zero }
override fun add(a: Point<T>, b: Point<T>): Point<T> = produce { with(space) { a[it] + b[it] } }
override fun multiply(a: Point<T>, k: Number): Point<T> = produce { with(space) { a[it] * k } }
//TODO add basis
companion object {
private val realSpaceCache = HashMap<Int, BufferVectorSpace<Double, RealField>>()
/**
* Non-boxing double vector space
*/
fun real(size: Int): BufferVectorSpace<Double, RealField> {
return realSpaceCache.getOrPut(size) {
BufferVectorSpace(
size,
RealField,
Buffer.Companion::auto
)
}
}
/**
* A structured vector space with custom buffer
*/
fun <T : Any, S : Space<T>> buffered(
size: Int,
space: S,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
) = BufferVectorSpace(size, space, bufferFactory)
/**
* Automatic buffered vector, unboxed if it is possible
*/
inline fun <reified T : Any, S : Space<T>> auto(size: Int, space: S): VectorSpace<T, S> =
buffered(size, space, Buffer.Companion::auto)
}
}
class BufferVectorSpace<T : Any, S : Space<T>>(
override val size: Int,
override val space: S,
val bufferFactory: BufferFactory<T>
) : VectorSpace<T, S> {
override fun produce(initializer: (Int) -> T) = bufferFactory(size, initializer)
//override fun produceElement(initializer: (Int) -> T): Vector<T, S> = BufferVector(this, produce(initializer))
}

3
kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDAlgebra.kt
View File

@ -73,6 +73,7 @@ interface NDSpace<T, S : Space<T>, N : NDStructure<T>> : Space<N>, NDAlgebra<T,
*/
override fun multiply(a: N, k: Number): N = map(a) { multiply(it, k) }
//TODO move to extensions after KEEP-176
operator fun N.plus(arg: T) = map(this) { value -> add(arg, value) }
operator fun N.minus(arg: T) = map(this) { value -> add(arg, -value) }
@ -90,6 +91,7 @@ interface NDRing<T, R : Ring<T>, N : NDStructure<T>> : Ring<N>, NDSpace<T, R, N>
*/
override fun multiply(a: N, b: N): N = combine(a, b) { aValue, bValue -> multiply(aValue, bValue) }
//TODO move to extensions after KEEP-176
operator fun N.times(arg: T) = map(this) { value -> multiply(arg, value) }
operator fun T.times(arg: N) = map(arg) { value -> multiply(this@times, value) }
}
@ -109,6 +111,7 @@ interface NDField<T, F : Field<T>, N : NDStructure<T>> : Field<N>, NDRing<T, F,
*/
override fun divide(a: N, b: N): N = combine(a, b) { aValue, bValue -> divide(aValue, bValue) }
//TODO move to extensions after KEEP-176
operator fun N.div(arg: T) = map(this) { value -> divide(arg, value) }
operator fun T.div(arg: N) = map(arg) { divide(it, this@div) }

9
kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt
View File

@ -1,5 +1,6 @@
package scientifik.kmath.linear
import scientifik.kmath.structures.Matrix
import kotlin.test.Test
import kotlin.test.assertEquals
@ -16,7 +17,7 @@ class MatrixTest {
@Test
fun testVectorToMatrix() {
val vector = Vector.real(5) { it.toDouble() }
val matrix = vector.toMatrix()
val matrix = vector.asMatrix()
assertEquals(4.0, matrix[4, 0])
}
@ -33,8 +34,8 @@ class MatrixTest {
val vector1 = Vector.real(5) { it.toDouble() }
val vector2 = Vector.real(5) { 5 - it.toDouble() }
val matrix1 = vector1.toMatrix()
val matrix2 = vector2.toMatrix().transpose()
val matrix1 = vector1.asMatrix()
val matrix2 = vector2.asMatrix().transpose()
val product = MatrixContext.real.run { matrix1 dot matrix2 }
@ -44,7 +45,7 @@ class MatrixTest {
@Test
fun testBuilder() {
val matrix = FeaturedMatrix.build<Double>(2, 3)(
val matrix = Matrix.build<Double>(2, 3)(
1.0, 0.0, 0.0,
0.0, 1.0, 2.0
)

15
kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
View File

@ -1,25 +1,28 @@
package scientifik.kmath.linear
import scientifik.kmath.structures.Matrix
import kotlin.contracts.ExperimentalContracts
import kotlin.test.Test
import kotlin.test.assertEquals
@ExperimentalContracts
class RealLUSolverTest {
@Test
fun testInvertOne() {
val matrix = MatrixContext.real.one(2, 2)
val inverted = LUSolver.real.inverse(matrix)
val inverted = MatrixContext.real.inverse(matrix)
assertEquals(matrix, inverted)
}
@Test
fun testInvert() {
val matrix = FeaturedMatrix.square(
val matrix = Matrix.square(
3.0, 1.0,
1.0, 3.0
)
val decomposed = LUSolver.real.decompose(matrix)
val decomposition = decomposed.getFeature<LUPDecomposition<Double>>()!!
val decomposition = MatrixContext.real.lup(matrix)
//Check determinant
assertEquals(8.0, decomposition.determinant)
@ -29,9 +32,9 @@ class RealLUSolverTest {
assertEquals(decomposition.p dot matrix, decomposition.l dot decomposition.u)
}
val inverted = LUSolver.real.inverse(decomposed)
val inverted = MatrixContext.real.inverse(matrix)
val expected = FeaturedMatrix.square(
val expected = Matrix.square(
0.375, -0.125,
-0.125, 0.375
)

32
kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt
View File

@ -2,10 +2,14 @@ package scientifik.kmath.linear
import koma.extensions.fill
import koma.matrix.MatrixFactory
import scientifik.kmath.operations.Space
import scientifik.kmath.structures.Matrix
class KomaMatrixContext<T : Any>(val factory: MatrixFactory<koma.matrix.Matrix<T>>) : MatrixContext<T>,
LinearSolver<T> {
class KomaMatrixContext<T : Any>(
private val factory: MatrixFactory<koma.matrix.Matrix<T>>,
private val space: Space<T>
) :
MatrixContext<T> {
override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T) =
KomaMatrix(factory.zeros(rows, columns).fill(initializer))
@ -32,28 +36,38 @@ class KomaMatrixContext<T : Any>(val factory: MatrixFactory<koma.matrix.Matrix<T
override fun Matrix<T>.unaryMinus() =
KomaMatrix(this.toKoma().origin.unaryMinus())
override fun Matrix<T>.plus(b: Matrix<T>) =
KomaMatrix(this.toKoma().origin + b.toKoma().origin)
override fun add(a: Matrix<T>, b: Matrix<T>) =
KomaMatrix(a.toKoma().origin + b.toKoma().origin)
override fun Matrix<T>.minus(b: Matrix<T>) =
KomaMatrix(this.toKoma().origin - b.toKoma().origin)
override fun multiply(a: Matrix<T>, k: Number): Matrix<T> =
produce(a.rowNum, a.colNum) { i, j -> space.run { a[i, j] * k } }
override fun Matrix<T>.times(value: T) =
KomaMatrix(this.toKoma().origin * value)
companion object {
override fun solve(a: Matrix<T>, b: Matrix<T>) =
KomaMatrix(a.toKoma().origin.solve(b.toKoma().origin))
}
override fun inverse(a: Matrix<T>) =
KomaMatrix(a.toKoma().origin.inv())
}
fun <T : Any> KomaMatrixContext<T>.solve(a: Matrix<T>, b: Matrix<T>) =
KomaMatrix(a.toKoma().origin.solve(b.toKoma().origin))
fun <T : Any> KomaMatrixContext<T>.solve(a: Matrix<T>, b: Point<T>) =
KomaVector(a.toKoma().origin.solve(b.toKoma().origin))
fun <T : Any> KomaMatrixContext<T>.inverse(a: Matrix<T>) =
KomaMatrix(a.toKoma().origin.inv())
class KomaMatrix<T : Any>(val origin: koma.matrix.Matrix<T>, features: Set<MatrixFeature>? = null) : FeaturedMatrix<T> {
override val rowNum: Int get() = origin.numRows()
override val colNum: Int get() = origin.numCols()
override val shape: IntArray get() = intArrayOf(origin.numRows(),origin.numCols())
override val shape: IntArray get() = intArrayOf(origin.numRows(), origin.numCols())
override val features: Set<MatrixFeature> = features ?: setOf(
object : DeterminantFeature<T> {