WIP vector space refactor

This commit is contained in:
Alexander Nozik 2021-03-12 22:52:18 +03:00
parent f449bdd58f
commit 9bc8e8fbf9
25 changed files with 352 additions and 566 deletions

View File

@ -4,11 +4,11 @@ import kotlinx.benchmark.Benchmark
import kotlinx.benchmark.Blackhole import kotlinx.benchmark.Blackhole
import kotlinx.benchmark.Scope import kotlinx.benchmark.Scope
import kotlinx.benchmark.State import kotlinx.benchmark.State
import space.kscience.kmath.commons.linear.CMMatrixContext import space.kscience.kmath.commons.linear.CMLinearSpace
import space.kscience.kmath.ejml.EjmlMatrixContext import space.kscience.kmath.ejml.EjmlLinearSpace
import space.kscience.kmath.linear.BufferMatrixContext import space.kscience.kmath.linear.BufferLinearSpace
import space.kscience.kmath.linear.Matrix import space.kscience.kmath.linear.Matrix
import space.kscience.kmath.linear.RealMatrixContext import space.kscience.kmath.linear.RealLinearSpace
import space.kscience.kmath.operations.RealField import space.kscience.kmath.operations.RealField
import space.kscience.kmath.operations.invoke import space.kscience.kmath.operations.invoke
import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.Buffer
@ -24,44 +24,44 @@ internal class DotBenchmark {
val matrix1 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } val matrix1 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
val matrix2 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 } val matrix2 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
val cmMatrix1 = CMMatrixContext { matrix1.toCM() } val cmMatrix1 = CMLinearSpace { matrix1.toCM() }
val cmMatrix2 = CMMatrixContext { matrix2.toCM() } val cmMatrix2 = CMLinearSpace { matrix2.toCM() }
val ejmlMatrix1 = EjmlMatrixContext { matrix1.toEjml() } val ejmlMatrix1 = EjmlLinearSpace { matrix1.toEjml() }
val ejmlMatrix2 = EjmlMatrixContext { matrix2.toEjml() } val ejmlMatrix2 = EjmlLinearSpace { matrix2.toEjml() }
} }
@Benchmark @Benchmark
fun cmDot(blackhole: Blackhole) { fun cmDot(blackhole: Blackhole) {
CMMatrixContext { CMLinearSpace {
blackhole.consume(cmMatrix1 dot cmMatrix2) blackhole.consume(cmMatrix1 dot cmMatrix2)
} }
} }
@Benchmark @Benchmark
fun ejmlDot(blackhole: Blackhole) { fun ejmlDot(blackhole: Blackhole) {
EjmlMatrixContext { EjmlLinearSpace {
blackhole.consume(ejmlMatrix1 dot ejmlMatrix2) blackhole.consume(ejmlMatrix1 dot ejmlMatrix2)
} }
} }
@Benchmark @Benchmark
fun ejmlDotWithConversion(blackhole: Blackhole) { fun ejmlDotWithConversion(blackhole: Blackhole) {
EjmlMatrixContext { EjmlLinearSpace {
blackhole.consume(matrix1 dot matrix2) blackhole.consume(matrix1 dot matrix2)
} }
} }
@Benchmark @Benchmark
fun bufferedDot(blackhole: Blackhole) { fun bufferedDot(blackhole: Blackhole) {
BufferMatrixContext(RealField, Buffer.Companion::real).invoke { BufferLinearSpace(RealField, Buffer.Companion::real).invoke {
blackhole.consume(matrix1 dot matrix2) blackhole.consume(matrix1 dot matrix2)
} }
} }
@Benchmark @Benchmark
fun realDot(blackhole: Blackhole) { fun realDot(blackhole: Blackhole) {
RealMatrixContext { RealLinearSpace {
blackhole.consume(matrix1 dot matrix2) blackhole.consume(matrix1 dot matrix2)
} }
} }

View File

@ -4,13 +4,12 @@ import kotlinx.benchmark.Benchmark
import kotlinx.benchmark.Blackhole import kotlinx.benchmark.Blackhole
import kotlinx.benchmark.Scope import kotlinx.benchmark.Scope
import kotlinx.benchmark.State import kotlinx.benchmark.State
import space.kscience.kmath.commons.linear.CMMatrixContext import space.kscience.kmath.commons.linear.CMLinearSpace
import space.kscience.kmath.commons.linear.CMMatrixContext.dot
import space.kscience.kmath.commons.linear.inverse import space.kscience.kmath.commons.linear.inverse
import space.kscience.kmath.ejml.EjmlMatrixContext import space.kscience.kmath.ejml.EjmlLinearSpace
import space.kscience.kmath.ejml.inverse import space.kscience.kmath.ejml.inverse
import space.kscience.kmath.linear.LinearSpace
import space.kscience.kmath.linear.Matrix import space.kscience.kmath.linear.Matrix
import space.kscience.kmath.linear.MatrixContext
import space.kscience.kmath.linear.inverseWithLup import space.kscience.kmath.linear.inverseWithLup
import space.kscience.kmath.linear.real import space.kscience.kmath.linear.real
import kotlin.random.Random import kotlin.random.Random
@ -29,19 +28,19 @@ internal class LinearAlgebraBenchmark {
@Benchmark @Benchmark
fun kmathLupInversion(blackhole: Blackhole) { fun kmathLupInversion(blackhole: Blackhole) {
blackhole.consume(MatrixContext.real.inverseWithLup(matrix)) blackhole.consume(LinearSpace.real.inverseWithLup(matrix))
} }
@Benchmark @Benchmark
fun cmLUPInversion(blackhole: Blackhole) { fun cmLUPInversion(blackhole: Blackhole) {
with(CMMatrixContext) { with(CMLinearSpace) {
blackhole.consume(inverse(matrix)) blackhole.consume(inverse(matrix))
} }
} }
@Benchmark @Benchmark
fun ejmlInverse(blackhole: Blackhole) { fun ejmlInverse(blackhole: Blackhole) {
with(EjmlMatrixContext) { with(EjmlLinearSpace) {
blackhole.consume(inverse(matrix)) blackhole.consume(inverse(matrix))
} }
} }

View File

@ -7,7 +7,7 @@ import kotlin.system.measureTimeMillis
@Suppress("UNUSED_VARIABLE") @Suppress("UNUSED_VARIABLE")
fun main() { fun main() {
val n = 6000 val n = 6000
val structure = NDStructure.build(intArrayOf(n, n), Buffer.Companion::auto) { 1.0 } val structure = NDStructure.buffered(intArrayOf(n, n), Buffer.Companion::auto) { 1.0 }
structure.mapToBuffer { it + 1 } // warm-up structure.mapToBuffer { it + 1 } // warm-up
val time1 = measureTimeMillis { val res = structure.mapToBuffer { it + 1 } } val time1 = measureTimeMillis { val res = structure.mapToBuffer { it + 1 } }
println("Structure mapping finished in $time1 millis") println("Structure mapping finished in $time1 millis")

View File

@ -3,6 +3,7 @@ package space.kscience.kmath.commons.linear
import org.apache.commons.math3.linear.* import org.apache.commons.math3.linear.*
import space.kscience.kmath.linear.* import space.kscience.kmath.linear.*
import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.RealField
import space.kscience.kmath.structures.RealBuffer import space.kscience.kmath.structures.RealBuffer
import kotlin.reflect.KClass import kotlin.reflect.KClass
import kotlin.reflect.cast import kotlin.reflect.cast
@ -55,7 +56,7 @@ public inline class CMMatrix(public val origin: RealMatrix) : Matrix<Double> {
public fun RealMatrix.asMatrix(): CMMatrix = CMMatrix(this) public fun RealMatrix.asMatrix(): CMMatrix = CMMatrix(this)
public class CMVector(public val origin: RealVector) : Point<Double> { public class CMVector(public val origin: RealVector) : Vector<Double> {
public override val size: Int get() = origin.dimension public override val size: Int get() = origin.dimension
public override operator fun get(index: Int): Double = origin.getEntry(index) public override operator fun get(index: Int): Double = origin.getEntry(index)
@ -70,8 +71,8 @@ public fun Point<Double>.toCM(): CMVector = if (this is CMVector) this else {
public fun RealVector.toPoint(): CMVector = CMVector(this) public fun RealVector.toPoint(): CMVector = CMVector(this)
public object CMMatrixContext : MatrixContext<Double, CMMatrix> { public object CMLinearSpace : LinearSpace<Double, RealField> {
public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): CMMatrix { public override fun buildMatrix(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): CMMatrix {
val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } } val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } }
return CMMatrix(Array2DRowRealMatrix(array)) return CMMatrix(Array2DRowRealMatrix(array))
} }
@ -86,17 +87,15 @@ public object CMMatrixContext : MatrixContext<Double, CMMatrix> {
} }
} }
override fun scale(a: Matrix<Double>, value: Double): Matrix<Double> = a.toCM().times(value)
public override fun Matrix<Double>.dot(other: Matrix<Double>): CMMatrix = public override fun Matrix<Double>.dot(other: Matrix<Double>): CMMatrix =
CMMatrix(toCM().origin.multiply(other.toCM().origin)) CMMatrix(toCM().origin.multiply(other.toCM().origin))
public override fun Matrix<Double>.dot(vector: Point<Double>): CMVector = public override fun Matrix<Double>.dot(vector: Vector<Double>): CMVector =
CMVector(toCM().origin.preMultiply(vector.toCM().origin)) CMVector(toCM().origin.preMultiply(vector.toCM().origin))
public override operator fun Matrix<Double>.unaryMinus(): CMMatrix = public override operator fun Matrix<Double>.unaryMinus(): CMMatrix =
produce(rowNum, colNum) { i, j -> -get(i, j) } buildMatrix(rowNum, colNum) { i, j -> -get(i, j) }
public override fun add(a: Matrix<Double>, b: Matrix<Double>): CMMatrix = public override fun add(a: Matrix<Double>, b: Matrix<Double>): CMMatrix =
CMMatrix(a.toCM().origin.multiply(b.toCM().origin)) CMMatrix(a.toCM().origin.multiply(b.toCM().origin))
@ -108,7 +107,7 @@ public object CMMatrixContext : MatrixContext<Double, CMMatrix> {
// CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble())) // CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble()))
public override operator fun Matrix<Double>.times(value: Double): CMMatrix = public override operator fun Matrix<Double>.times(value: Double): CMMatrix =
produce(rowNum, colNum) { i, j -> get(i, j) * value } buildMatrix(rowNum, colNum) { i, j -> get(i, j) * value }
} }
public operator fun CMMatrix.plus(other: CMMatrix): CMMatrix = public operator fun CMMatrix.plus(other: CMMatrix): CMMatrix =

View File

@ -12,7 +12,7 @@ public enum class CMDecomposition {
CHOLESKY CHOLESKY
} }
public fun CMMatrixContext.solver( public fun CMLinearSpace.solver(
a: Matrix<Double>, a: Matrix<Double>,
decomposition: CMDecomposition = CMDecomposition.LUP decomposition: CMDecomposition = CMDecomposition.LUP
): DecompositionSolver = when (decomposition) { ): DecompositionSolver = when (decomposition) {
@ -23,19 +23,19 @@ public fun CMMatrixContext.solver(
CMDecomposition.CHOLESKY -> CholeskyDecomposition(a.toCM().origin).solver CMDecomposition.CHOLESKY -> CholeskyDecomposition(a.toCM().origin).solver
} }
public fun CMMatrixContext.solve( public fun CMLinearSpace.solve(
a: Matrix<Double>, a: Matrix<Double>,
b: Matrix<Double>, b: Matrix<Double>,
decomposition: CMDecomposition = CMDecomposition.LUP decomposition: CMDecomposition = CMDecomposition.LUP
): CMMatrix = solver(a, decomposition).solve(b.toCM().origin).asMatrix() ): CMMatrix = solver(a, decomposition).solve(b.toCM().origin).asMatrix()
public fun CMMatrixContext.solve( public fun CMLinearSpace.solve(
a: Matrix<Double>, a: Matrix<Double>,
b: Point<Double>, b: Point<Double>,
decomposition: CMDecomposition = CMDecomposition.LUP decomposition: CMDecomposition = CMDecomposition.LUP
): CMVector = solver(a, decomposition).solve(b.toCM().origin).toPoint() ): CMVector = solver(a, decomposition).solve(b.toCM().origin).toPoint()
public fun CMMatrixContext.inverse( public fun CMLinearSpace.inverse(
a: Matrix<Double>, a: Matrix<Double>,
decomposition: CMDecomposition = CMDecomposition.LUP decomposition: CMDecomposition = CMDecomposition.LUP
): CMMatrix = solver(a, decomposition).inverse.asMatrix() ): CMMatrix = solver(a, decomposition).inverse.asMatrix()

View File

@ -1,144 +0,0 @@
package space.kscience.kmath.linear
import space.kscience.kmath.nd.NDStructure
import space.kscience.kmath.nd.Structure2D
import space.kscience.kmath.operations.Ring
import space.kscience.kmath.operations.ScaleOperations
import space.kscience.kmath.operations.invoke
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.BufferFactory
import space.kscience.kmath.structures.asSequence
/**
* Alias for [Structure2D] with more familiar name.
*
* @param T the type of items.
*/
public typealias Matrix<T> = Structure2D<T>
/**
* Basic implementation of Matrix space based on [NDStructure]
*/
public class BufferMatrixContext<T : Any, A>(
public override val elementContext: A,
private val bufferFactory: BufferFactory<T>,
) : GenericMatrixContext<T, A, BufferMatrix<T>> where A : Ring<T>, A : ScaleOperations<T> {
public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): BufferMatrix<T> {
val buffer = bufferFactory(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
return BufferMatrix(rows, columns, buffer)
}
override fun scale(a: Matrix<T>, value: Double): Matrix<T> = elementContext {
produce(a.rowNum, a.colNum) { i, j ->
a[i, j] * value
}
}
public override fun point(size: Int, initializer: (Int) -> T): Point<T> = bufferFactory(size, initializer)
private fun Matrix<T>.toBufferMatrix(): BufferMatrix<T> = if (this is BufferMatrix) this else {
produce(rowNum, colNum) { i, j -> get(i, j) }
}
public fun one(rows: Int, columns: Int): Matrix<Double> = VirtualMatrix(rows, columns) { i, j ->
if (i == j) 1.0 else 0.0
} + DiagonalFeature
public override infix fun Matrix<T>.dot(other: Matrix<T>): BufferMatrix<T> {
require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
val bufferMatrix = toBufferMatrix()
val otherBufferMatrix = other.toBufferMatrix()
return elementContext {
produce(rowNum, other.colNum) { i, j ->
var res = one
for (l in 0 until colNum) {
res += bufferMatrix[i, l] * otherBufferMatrix[l, j]
}
res
}
}
}
public override infix fun Matrix<T>.dot(vector: Point<T>): Point<T> {
require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
val bufferMatrix = toBufferMatrix()
return elementContext {
bufferFactory(rowNum) { i ->
var res = one
for (j in 0 until colNum) {
res += bufferMatrix[i, j] * vector[j]
}
res
}
}
}
override fun add(a: Matrix<T>, b: Matrix<T>): BufferMatrix<T> {
require(a.rowNum == b.rowNum) { "Row number mismatch in matrix addition. Left side: ${a.rowNum}, right side: ${b.rowNum}" }
require(a.colNum == b.colNum) { "Column number mismatch in matrix addition. Left side: ${a.colNum}, right side: ${b.colNum}" }
val aBufferMatrix = a.toBufferMatrix()
val bBufferMatrix = b.toBufferMatrix()
return elementContext {
produce(a.rowNum, a.colNum) { i, j ->
aBufferMatrix[i, j] + bBufferMatrix[i, j]
}
}
}
// override fun multiply(a: Matrix<T>, k: Number): BufferMatrix<T> {
// val aBufferMatrix = a.toBufferMatrix()
// return elementContext {
// produce(a.rowNum, a.colNum) { i, j -> aBufferMatrix[i, j] * k.toDouble() }
// }
// }
public companion object
}
public class BufferMatrix<T : Any>(
public override val rowNum: Int,
public override val colNum: Int,
public val buffer: Buffer<T>,
) : Matrix<T> {
init {
require(buffer.size == rowNum * colNum) { "Dimension mismatch for matrix structure" }
}
override val shape: IntArray get() = intArrayOf(rowNum, colNum)
public override operator fun get(index: IntArray): T = get(index[0], index[1])
public override operator fun get(i: Int, j: Int): T = buffer[i * colNum + j]
public override fun elements(): Sequence<Pair<IntArray, T>> = sequence {
for (i in 0 until rowNum) for (j in 0 until colNum) yield(intArrayOf(i, j) to get(i, j))
}
public override fun equals(other: Any?): Boolean {
if (this === other) return true
return when (other) {
is NDStructure<*> -> NDStructure.contentEquals(this, other)
else -> false
}
}
override fun hashCode(): Int {
var result = rowNum
result = 31 * result + colNum
result = 31 * result + buffer.hashCode()
return result
}
public override fun toString(): String {
return if (rowNum <= 5 && colNum <= 5)
"Matrix(rowsNum = $rowNum, colNum = $colNum)\n" +
rows.asSequence().joinToString(prefix = "(", postfix = ")", separator = "\n ") { buffer ->
buffer.asSequence().joinToString(separator = "\t") { it.toString() }
}
else "Matrix(rowsNum = $rowNum, colNum = $colNum)"
}
}

View File

@ -1,7 +1,7 @@
package space.kscience.kmath.linear package space.kscience.kmath.linear
import space.kscience.kmath.nd.as1D
import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.VirtualBuffer
public typealias Point<T> = Buffer<T> public typealias Point<T> = Buffer<T>
@ -10,16 +10,16 @@ public typealias Point<T> = Buffer<T>
*/ */
public interface LinearSolver<T : Any> { public interface LinearSolver<T : Any> {
public fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T> public fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T>
public fun solve(a: Matrix<T>, b: Point<T>): Point<T> = solve(a, b.asMatrix()).asPoint() public fun solve(a: Matrix<T>, b: Point<T>): Point<T> = solve(a, b.asMatrix()).asVector()
public fun inverse(a: Matrix<T>): Matrix<T> public fun inverse(a: Matrix<T>): Matrix<T>
} }
/** /**
* Convert matrix to vector if it is possible * Convert matrix to vector if it is possible
*/ */
public fun <T : Any> Matrix<T>.asPoint(): Point<T> = public fun <T : Any> Matrix<T>.asVector(): Vector<T> =
if (this.colNum == 1) if (this.colNum == 1)
VirtualBuffer(rowNum) { get(it, 0) } as1D()
else else
error("Can't convert matrix with more than one column to vector") error("Can't convert matrix with more than one column to vector")

View File

@ -0,0 +1,213 @@
package space.kscience.kmath.linear
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.nd.*
import space.kscience.kmath.operations.*
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.BufferFactory
import kotlin.reflect.KClass
/**
* Alias for [Structure2D] with more familiar name.
*
* @param T the type of items.
*/
public typealias Matrix<T> = Structure2D<T>
/**
* Alias for [Structure1D] with more familiar name.
*
* @param T the type of items.
*/
public typealias Vector<T> = Structure1D<T>
/**
* Basic operations on matrices and vectors. Operates on [Matrix].
*
* @param T the type of items in the matrices.
* @param M the type of operated matrices.
*/
public interface LinearSpace<T : Any, A : Ring<T>> {
public val elementAlgebra: A
/**
* Produces a matrix with this context and given dimensions.
*/
public fun buildMatrix(rows: Int, columns: Int, initializer: A.(i: Int, j: Int) -> T): Matrix<T>
/**
* Produces a point compatible with matrix space (and possibly optimized for it).
*/
public fun buildVector(size: Int, initializer: A.(Int) -> T): Vector<T> =
buildMatrix(1, size) { _, j -> initializer(j) }.as1D()
public operator fun Matrix<T>.unaryMinus(): Matrix<T> = buildMatrix(rowNum, colNum) { i, j ->
-get(i, j)
}
public operator fun Vector<T>.unaryMinus(): Vector<T> = buildVector(size) {
-get(it)
}
/**
* Matrix sum
*/
public operator fun Matrix<T>.plus(other: Matrix<T>): Matrix<T> = buildMatrix(rowNum, colNum) { i, j ->
get(i, j) + other[i, j]
}
/**
* Vector sum
*/
public operator fun Vector<T>.plus(other: Vector<T>): Vector<T> = buildVector(size) {
get(it) + other[it]
}
/**
* Matrix subtraction
*/
public operator fun Matrix<T>.minus(other: Matrix<T>): Matrix<T> = buildMatrix(rowNum, colNum) { i, j ->
get(i, j) - other[i, j]
}
/**
* Vector subtraction
*/
public operator fun Vector<T>.minus(other: Vector<T>): Vector<T> = buildVector(size) {
get(it) - other[it]
}
/**
* 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<T>.dot(other: Matrix<T>): Matrix<T> {
require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
return elementAlgebra {
buildMatrix(rowNum, other.colNum) { i, j ->
var res = zero
for (l in 0 until colNum) {
res += this@dot[i, l] * other[l, j]
}
res
}
}
}
/**
* 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<T>.dot(vector: Vector<T>): Vector<T> {
require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
return elementAlgebra {
buildVector(rowNum) { i ->
var res = one
for (j in 0 until colNum) {
res += this@dot[i, j] * vector[j]
}
res
}
}
}
/**
* Multiplies a matrix by its element.
*
* @receiver the multiplicand.
* @param value the multiplier.
* @receiver the product.
*/
public operator fun Matrix<T>.times(value: T): Matrix<T> =
buildMatrix(rowNum, colNum) { i, j -> get(i, j) * value }
/**
* Multiplies an element by a matrix of it.
*
* @receiver the multiplicand.
* @param m the multiplier.
* @receiver the product.
*/
public operator fun T.times(m: Matrix<T>): Matrix<T> = m * this
/**
* Multiplies a vector by its element.
*
* @receiver the multiplicand.
* @param value the multiplier.
* @receiver the product.
*/
public operator fun Vector<T>.times(value: T): Vector<T> =
buildVector(size) { i -> get(i) * value }
/**
* Multiplies an element by a vector of it.
*
* @receiver the multiplicand.
* @param v the multiplier.
* @receiver the product.
*/
public operator fun T.times(v: Vector<T>): Vector<T> = v * this
/**
* Gets a feature from the matrix. This function may return some additional features to
* [space.kscience.kmath.nd.NDStructure.getFeature].
*
* @param F the type of feature.
* @param m the matrix.
* @param type the [KClass] instance of [F].
* @return a feature object or `null` if it isn't present.
*/
@UnstableKMathAPI
public fun <F : Any> getFeature(m: Matrix<T>, type: KClass<F>): F? = m.getFeature(type)
public companion object {
/**
* A structured matrix with custom buffer
*/
public fun <T : Any, A : Ring<T>> buffered(
algebra: A,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
): LinearSpace<T, A> = object : LinearSpace<T, A> {
override val elementAlgebra: A = algebra
override fun buildMatrix(
rows: Int, columns: Int,
initializer: A.(i: Int, j: Int) -> T,
): Matrix<T> = NDStructure.buffered(intArrayOf(rows, columns)) { (i, j) ->
algebra.initializer(i, j)
}.as2D()
}
/**
* Automatic buffered matrix, unboxed if it is possible
*/
public inline fun <reified T : Any, A : Ring<T>> auto(ring: A): LinearSpace<T, A> =
buffered(ring, Buffer.Companion::auto)
}
}
/**
* Gets a feature from the matrix. This function may return some additional features to
* [space.kscience.kmath.nd.NDStructure.getFeature].
*
* @param T the type of items in the matrices.
* @param M the type of operated matrices.
* @param F the type of feature.
* @receiver the [LinearSpace] of [T].
* @param m the matrix.
* @return a feature object or `null` if it isn't present.
*/
@UnstableKMathAPI
public inline fun <T : Any, reified F : Any> LinearSpace<T, *>.getFeature(m: Matrix<T>): F? = getFeature(m, F::class)

View File

@ -12,7 +12,7 @@ import space.kscience.kmath.structures.MutableBufferFactory
* Common implementation of [LupDecompositionFeature]. * Common implementation of [LupDecompositionFeature].
*/ */
public class LupDecomposition<T : Any>( public class LupDecomposition<T : Any>(
public val context: MatrixContext<T, Matrix<T>>, public val context: LinearSpace<T, *>,
public val elementContext: Field<T>, public val elementContext: Field<T>,
public val lu: Matrix<T>, public val lu: Matrix<T>,
public val pivot: IntArray, public val pivot: IntArray,
@ -62,15 +62,14 @@ public class LupDecomposition<T : Any>(
} }
@PublishedApi @PublishedApi
internal fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, *>.abs(value: T): T = internal fun <T : Comparable<T>> LinearSpace<T, Ring<T>>.abs(value: T): T =
if (value > elementContext.zero) value else elementContext { -value } if (value > elementAlgebra.zero) value else elementAlgebra { -value }
/** /**
* Create a lup decomposition of generic matrix. * Create a lup decomposition of generic matrix.
*/ */
public fun <T : Comparable<T>> MatrixContext<T, Matrix<T>>.lup( public fun <T : Comparable<T>> LinearSpace<T, Field<T>>.lup(
factory: MutableBufferFactory<T>, factory: MutableBufferFactory<T>,
elementContext: Field<T>,
matrix: Matrix<T>, matrix: Matrix<T>,
checkSingular: (T) -> Boolean, checkSingular: (T) -> Boolean,
): LupDecomposition<T> { ): LupDecomposition<T> {
@ -80,7 +79,7 @@ public fun <T : Comparable<T>> MatrixContext<T, Matrix<T>>.lup(
//TODO just waits for KEEP-176 //TODO just waits for KEEP-176
BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run { BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {
elementContext { elementAlgebra {
val lu = create(matrix) val lu = create(matrix)
// Initialize permutation array and parity // Initialize permutation array and parity
@ -142,18 +141,18 @@ public fun <T : Comparable<T>> MatrixContext<T, Matrix<T>>.lup(
for (row in col + 1 until m) lu[row, col] /= luDiag for (row in col + 1 until m) lu[row, col] /= luDiag
} }
return LupDecomposition(this@lup, elementContext, lu.collect(), pivot, even) return LupDecomposition(this@lup, elementAlgebra, lu.collect(), pivot, even)
} }
} }
} }
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, Matrix<T>>.lup( public inline fun <reified T : Comparable<T>> LinearSpace<T, Field<T>>.lup(
matrix: Matrix<T>, matrix: Matrix<T>,
noinline checkSingular: (T) -> Boolean, noinline checkSingular: (T) -> Boolean,
): LupDecomposition<T> = lup(MutableBuffer.Companion::auto, elementContext, matrix, checkSingular) ): LupDecomposition<T> = lup(MutableBuffer.Companion::auto, matrix, checkSingular)
public fun MatrixContext<Double, Matrix<Double>>.lup(matrix: Matrix<Double>): LupDecomposition<Double> = public fun LinearSpace<Double, RealField>.lup(matrix: Matrix<Double>): LupDecomposition<Double> =
lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 } lup(Buffer.Companion::real, matrix) { it < 1e-11 }
public fun <T : Any> LupDecomposition<T>.solveWithLup( public fun <T : Any> LupDecomposition<T>.solveWithLup(
factory: MutableBufferFactory<T>, factory: MutableBufferFactory<T>,
@ -198,7 +197,7 @@ public fun <T : Any> LupDecomposition<T>.solveWithLup(
} }
} }
return context.produce(pivot.size, matrix.colNum) { i, j -> bp[i, j] } return context.buildMatrix(pivot.size, matrix.colNum) { i, j -> bp[i, j] }
} }
} }
} }
@ -210,18 +209,18 @@ public inline fun <reified T : Any> LupDecomposition<T>.solveWithLup(matrix: Mat
* Solves a system of linear equations *ax = b** using LUP decomposition. * Solves a system of linear equations *ax = b** using LUP decomposition.
*/ */
@OptIn(UnstableKMathAPI::class) @OptIn(UnstableKMathAPI::class)
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, Matrix<T>>.solveWithLup( public inline fun <reified T : Comparable<T>> LinearSpace<T, Field<T>>.solveWithLup(
a: Matrix<T>, a: Matrix<T>,
b: Matrix<T>, b: Matrix<T>,
noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto, noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
noinline checkSingular: (T) -> Boolean, noinline checkSingular: (T) -> Boolean,
): Matrix<T> { ): Matrix<T> {
// Use existing decomposition if it is provided by matrix // Use existing decomposition if it is provided by matrix
val decomposition = a.getFeature() ?: lup(bufferFactory, elementContext, a, checkSingular) val decomposition = a.getFeature() ?: lup(bufferFactory, a, checkSingular)
return decomposition.solveWithLup(bufferFactory, b) return decomposition.solveWithLup(bufferFactory, b)
} }
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, Matrix<T>>.inverseWithLup( public inline fun <reified T : Comparable<T>> LinearSpace<T, Field<T>>.inverseWithLup(
matrix: Matrix<T>, matrix: Matrix<T>,
noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto, noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
noinline checkSingular: (T) -> Boolean, noinline checkSingular: (T) -> Boolean,
@ -229,15 +228,15 @@ public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext
@OptIn(UnstableKMathAPI::class) @OptIn(UnstableKMathAPI::class)
public fun RealMatrixContext.solveWithLup(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> { public fun RealLinearSpace.solveWithLup(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> {
// Use existing decomposition if it is provided by matrix // Use existing decomposition if it is provided by matrix
val bufferFactory: MutableBufferFactory<Double> = MutableBuffer.Companion::real val bufferFactory: MutableBufferFactory<Double> = MutableBuffer.Companion::real
val decomposition: LupDecomposition<Double> = a.getFeature() ?: lup(bufferFactory, RealField, a) { it < 1e-11 } val decomposition: LupDecomposition<Double> = a.getFeature() ?: lup(bufferFactory, a) { it < 1e-11 }
return decomposition.solveWithLup(bufferFactory, b) return decomposition.solveWithLup(bufferFactory, b)
} }
/** /**
* Inverses a square matrix using LUP decomposition. Non square matrix will throw a error. * Inverses a square matrix using LUP decomposition. Non square matrix will throw a error.
*/ */
public fun RealMatrixContext.inverseWithLup(matrix: Matrix<Double>): Matrix<Double> = public fun RealLinearSpace.inverseWithLup(matrix: Matrix<Double>): Matrix<Double> =
solveWithLup(matrix, one(matrix.rowNum, matrix.colNum)) solveWithLup(matrix, one(matrix.rowNum, matrix.colNum))

View File

@ -1,46 +1,30 @@
package space.kscience.kmath.linear package space.kscience.kmath.linear
import space.kscience.kmath.nd.Structure2D import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.structures.Buffer import space.kscience.kmath.operations.Ring
import space.kscience.kmath.structures.BufferFactory
import space.kscience.kmath.structures.asBuffer
public class MatrixBuilder(public val rows: Int, public val columns: Int) {
public operator fun <T : Any> invoke(vararg elements: T): Matrix<T> {
require(rows * columns == elements.size) { "The number of elements ${elements.size} is not equal $rows * $columns" }
val buffer = elements.asBuffer()
return BufferMatrix(rows, columns, buffer)
}
//TODO add specific matrix builder functions like diagonal, etc @UnstableKMathAPI
public fun <T : Any> LinearSpace<T, Ring<T>>.matrix(rows: Int, columns: Int, vararg elements: T): Matrix<T> {
require(rows * columns == elements.size) { "The number of elements ${elements.size} is not equal $rows * $columns" }
return buildMatrix(rows, columns) { i, j -> elements[i * columns + j] }
} }
public fun Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder = MatrixBuilder(rows, columns) @UnstableKMathAPI
public fun <T : Any> LinearSpace<T, Ring<T>>.vector(vararg elements: T): Vector<T> {
public fun <T : Any> Structure2D.Companion.row(vararg values: T): Matrix<T> { return buildVector(elements.size, elements::get)
val buffer = values.asBuffer()
return BufferMatrix(1, values.size, buffer)
} }
public inline fun <reified T : Any> Structure2D.Companion.row( public inline fun <T : Any> LinearSpace<T, Ring<T>>.row(
size: Int, size: Int,
factory: BufferFactory<T> = Buffer.Companion::auto, crossinline builder: (Int) -> T,
noinline builder: (Int) -> T, ): Matrix<T> = buildMatrix(1, size) { _, j -> builder(j) }
): Matrix<T> {
val buffer = factory(size, builder)
return BufferMatrix(1, size, buffer)
}
public fun <T : Any> Structure2D.Companion.column(vararg values: T): Matrix<T> { public fun <T : Any> LinearSpace<T, Ring<T>>.row(vararg values: T): Matrix<T> = row(values.size, values::get)
val buffer = values.asBuffer()
return BufferMatrix(values.size, 1, buffer)
}
public inline fun <reified T : Any> Structure2D.Companion.column( public inline fun <T : Any> LinearSpace<T, Ring<T>>.column(
size: Int, size: Int,
factory: BufferFactory<T> = Buffer.Companion::auto, crossinline builder: (Int) -> T,
noinline builder: (Int) -> T, ): Matrix<T> = buildMatrix(size, 1) { i, _ -> builder(i) }
): Matrix<T> {
val buffer = factory(size, builder) public fun <T : Any> LinearSpace<T, Ring<T>>.column(vararg values: T): Matrix<T> = column(values.size, values::get)
return BufferMatrix(size, 1, buffer)
}

View File

@ -1,173 +0,0 @@
package space.kscience.kmath.linear
import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.*
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.BufferFactory
import space.kscience.kmath.structures.asSequence
import kotlin.reflect.KClass
/**
* Basic operations on matrices. Operates on [Matrix].
*
* @param T the type of items in the matrices.
* @param M the type of operated matrices.
*/
public interface MatrixContext<T : Any, out M : Matrix<T>> : GroupOperations<Matrix<T>>, ScaleOperations<Matrix<T>> {
/**
* Produces a matrix with this context and given dimensions.
*/
public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): M
/**
* Produces a point compatible with matrix space (and possibly optimized for it).
*/
public fun point(size: Int, initializer: (Int) -> T): Point<T> = Buffer.boxing(size, initializer)
@Suppress("UNCHECKED_CAST")
public override fun binaryOperationFunction(operation: String): (left: Matrix<T>, right: Matrix<T>) -> M =
when (operation) {
"dot" -> { left, right -> left dot right }
else -> super<GroupOperations>.binaryOperationFunction(operation) as (Matrix<T>, Matrix<T>) -> M
}
/**
* 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<T>.dot(other: Matrix<T>): M
/**
* 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<T>.dot(vector: Point<T>): Point<T>
/**
* Multiplies a matrix by its element.
*
* @receiver the multiplicand.
* @param value the multiplier.
* @receiver the product.
*/
public operator fun Matrix<T>.times(value: T): M
/**
* Multiplies an element by a matrix of it.
*
* @receiver the multiplicand.
* @param m the multiplier.
* @receiver the product.
*/
public operator fun T.times(m: Matrix<T>): M = m * this
/**
* Gets a feature from the matrix. This function may return some additional features to
* [kscience.kmath.nd.NDStructure.getFeature].
*
* @param F the type of feature.
* @param m the matrix.
* @param type the [KClass] instance of [F].
* @return a feature object or `null` if it isn't present.
*/
@UnstableKMathAPI
public fun <F : Any> getFeature(m: Matrix<T>, type: KClass<F>): F? = m.getFeature(type)
public companion object {
/**
* A structured matrix with custom buffer
*/
public fun <T : Any, A> buffered(
ring: A,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
): GenericMatrixContext<T, A, BufferMatrix<T>> where A : Ring<T>, A: ScaleOperations<T> = BufferMatrixContext(ring, bufferFactory)
/**
* Automatic buffered matrix, unboxed if it is possible
*/
public inline fun <reified T : Any, A> auto(ring: A): GenericMatrixContext<T, A, BufferMatrix<T>> where A : Ring<T>, A: ScaleOperations<T> =
buffered(ring, Buffer.Companion::auto)
}
}
/**
* Gets a feature from the matrix. This function may return some additional features to
* [kscience.kmath.nd.NDStructure.getFeature].
*
* @param T the type of items in the matrices.
* @param M the type of operated matrices.
* @param F the type of feature.
* @receiver the [MatrixContext] of [T].
* @param m the matrix.
* @return a feature object or `null` if it isn't present.
*/
@UnstableKMathAPI
public inline fun <T : Any, reified F : Any> MatrixContext<T, *>.getFeature(m: Matrix<T>): F? =
getFeature(m, F::class)
/**
* Partial implementation of [MatrixContext] for matrices of [Ring].
*
* @param T the type of items in the matrices.
* @param A the type of ring of matrix elements.
* @param M the type of operated matrices.
*/
public interface GenericMatrixContext<T : Any, A, out M : Matrix<T>> : MatrixContext<T, M> where A : Ring<T>, A : ScaleOperations<T>{
/**
* The ring over matrix elements.
*/
public val elementContext: A
public override infix fun Matrix<T>.dot(other: Matrix<T>): M {
//TODO add typed error
require(colNum == other.rowNum) { "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]
elementContext { sum(row.asSequence().zip(column.asSequence(), ::multiply)) }
}
}
public override infix fun Matrix<T>.dot(vector: Point<T>): Point<T> {
//TODO add typed error
require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
return point(rowNum) { i ->
val row = rows[i]
elementContext { sum(row.asSequence().zip(vector.asSequence(), ::multiply)) }
}
}
public override operator fun Matrix<T>.unaryMinus(): M =
produce(rowNum, colNum) { i, j -> elementContext { -get(i, j) } }
public override fun add(a: Matrix<T>, b: Matrix<T>): M {
require(a.rowNum == b.rowNum && a.colNum == b.colNum) {
"Matrix operation dimension mismatch. [${a.rowNum},${a.colNum}] + [${b.rowNum},${b.colNum}]"
}
return produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] + b[i, j] } }
}
public override operator fun Matrix<T>.minus(b: Matrix<T>): M {
require(rowNum == b.rowNum && colNum == b.colNum) {
"Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]"
}
return produce(rowNum, colNum) { i, j -> elementContext { get(i, j) + b[i, j] } }
}
//
// public override fun multiply(a: Matrix<T>, k: Number): M =
// produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] * k } }
public override operator fun Matrix<T>.times(value: T): M =
produce(rowNum, colNum) { i, j -> elementContext { get(i, j) * value } }
}

View File

@ -1,14 +1,9 @@
package space.kscience.kmath.linear package space.kscience.kmath.linear
import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.nd.Structure2D
import space.kscience.kmath.nd.getFeature import space.kscience.kmath.nd.getFeature
import space.kscience.kmath.operations.Ring import space.kscience.kmath.operations.Ring
import space.kscience.kmath.operations.ScaleOperations
import space.kscience.kmath.structures.asBuffer
import kotlin.math.sqrt
import kotlin.reflect.KClass import kotlin.reflect.KClass
import kotlin.reflect.safeCast
/** /**
* A [Matrix] that holds [MatrixFeature] objects. * A [Matrix] that holds [MatrixFeature] objects.
@ -24,7 +19,8 @@ public class MatrixWrapper<T : Any> internal constructor(
* Get the first feature matching given class. Does not guarantee that matrix has only one feature matching the criteria * Get the first feature matching given class. Does not guarantee that matrix has only one feature matching the criteria
*/ */
@UnstableKMathAPI @UnstableKMathAPI
override fun <T : Any> getFeature(type: KClass<T>): T? = type.safeCast(features.find { type.isInstance(it) }) @Suppress("UNCHECKED_CAST")
override fun <T : Any> getFeature(type: KClass<T>): T? = features.singleOrNull { type.isInstance(it) } as? T
?: origin.getFeature(type) ?: origin.getFeature(type)
override fun equals(other: Any?): Boolean = origin == other override fun equals(other: Any?): Boolean = origin == other
@ -61,35 +57,25 @@ public operator fun <T : Any> Matrix<T>.plus(newFeatures: Collection<MatrixFeatu
MatrixWrapper(this, newFeatures.toSet()) MatrixWrapper(this, newFeatures.toSet())
} }
/**
* Build a square matrix from given elements.
*/
public fun <T : Any> Structure2D.Companion.square(vararg elements: T): Matrix<T> {
val size: Int = sqrt(elements.size.toDouble()).toInt()
require(size * size == elements.size) { "The number of elements ${elements.size} is not a full square" }
val buffer = elements.asBuffer()
return BufferMatrix(size, size, buffer)
}
/** /**
* Diagonal matrix of ones. The matrix is virtual no actual matrix is created * Diagonal matrix of ones. The matrix is virtual no actual matrix is created
*/ */
public fun <T : Any, A> GenericMatrixContext<T, A, *>.one( public fun <T : Any> LinearSpace<T, Ring<T>>.one(
rows: Int, rows: Int,
columns: Int, columns: Int,
): Matrix<T> where A : Ring<T>, A : ScaleOperations<T> = VirtualMatrix(rows, columns) { i, j -> ): Matrix<T> = VirtualMatrix(rows, columns) { i, j ->
if (i == j) elementContext.one else elementContext.zero if (i == j) elementAlgebra.one else elementAlgebra.zero
} + UnitFeature } + UnitFeature
/** /**
* A virtual matrix of zeroes * A virtual matrix of zeroes
*/ */
public fun <T : Any, A> GenericMatrixContext<T, A, *>.zero( public fun <T : Any> LinearSpace<T, Ring<T>>.zero(
rows: Int, rows: Int,
columns: Int, columns: Int,
): Matrix<T> where A : Ring<T>, A : ScaleOperations<T> = VirtualMatrix(rows, columns) { _, _ -> ): Matrix<T> = VirtualMatrix(rows, columns) { _, _ ->
elementContext.zero elementAlgebra.zero
} + ZeroFeature } + ZeroFeature
public class TransposedFeature<T : Any>(public val original: Matrix<T>) : MatrixFeature public class TransposedFeature<T : Any>(public val original: Matrix<T>) : MatrixFeature

View File

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

View File

@ -1,72 +0,0 @@
package space.kscience.kmath.linear
import space.kscience.kmath.operations.Group
import space.kscience.kmath.operations.RealField
import space.kscience.kmath.operations.ScaleOperations
import space.kscience.kmath.operations.invoke
import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.BufferFactory
/**
* A linear space for vectors.
* Could be used on any point-like structure
*/
public interface VectorSpace<T : Any, A> : Group<Point<T>>, ScaleOperations<Point<T>>
where A : Group<T>, A : ScaleOperations<T> {
public val size: Int
public val algebra: A
override val zero: Point<T> get() = produce { algebra.zero }
public fun produce(initializer: A.(Int) -> T): Point<T>
override fun add(a: Point<T>, b: Point<T>): Point<T> = produce { algebra { a[it] + b[it] } }
override fun scale(a: Point<T>, value: Double): Point<T> = produce { algebra.scale(a[it], value) }
override fun Point<T>.unaryMinus(): Point<T> = produce { -get(it) }
//TODO add basis
public companion object {
private val realSpaceCache: MutableMap<Int, BufferVectorSpace<Double, RealField>> = hashMapOf()
/**
* Non-boxing double vector space
*/
public fun real(size: Int): BufferVectorSpace<Double, RealField> = realSpaceCache.getOrPut(size) {
BufferVectorSpace(
size,
RealField,
Buffer.Companion::auto
)
}
/**
* A structured vector space with custom buffer
*/
public fun <T : Any, A> buffered(
size: Int,
space: A,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
): BufferVectorSpace<T, A> where A : Group<T>, A : ScaleOperations<T> =
BufferVectorSpace(size, space, bufferFactory)
/**
* Automatic buffered vector, unboxed if it is possible
*/
public inline fun <reified T : Any, A> auto(
size: Int,
space: A,
): VectorSpace<T, A> where A : Group<T>, A : ScaleOperations<T> =
buffered(size, space, Buffer.Companion::auto)
}
}
public class BufferVectorSpace<T : Any, A>(
override val size: Int,
override val algebra: A,
public val bufferFactory: BufferFactory<T>,
) : VectorSpace<T, A> where A : Group<T>, A : ScaleOperations<T> {
override fun produce(initializer: A.(Int) -> T): Buffer<T> = bufferFactory(size) { algebra.initializer(it) }
}

View File

@ -31,6 +31,4 @@ public class VirtualMatrix<T : Any>(
result = 31 * result + generator.hashCode() result = 31 * result + generator.hashCode()
return result return result
} }
} }

View File

@ -74,7 +74,7 @@ public interface NDStructure<T> {
* *
* Strides should be reused if possible. * Strides should be reused if possible.
*/ */
public fun <T> build( public fun <T> buffered(
strides: Strides, strides: Strides,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing, bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
initializer: (IntArray) -> T, initializer: (IntArray) -> T,
@ -94,11 +94,11 @@ public interface NDStructure<T> {
crossinline initializer: (IntArray) -> T, crossinline initializer: (IntArray) -> T,
): NDBuffer<T> = NDBuffer(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) }) ): NDBuffer<T> = NDBuffer(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) })
public fun <T> build( public fun <T> buffered(
shape: IntArray, shape: IntArray,
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing, bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
initializer: (IntArray) -> T, initializer: (IntArray) -> T,
): NDBuffer<T> = build(DefaultStrides(shape), bufferFactory, initializer) ): NDBuffer<T> = buffered(DefaultStrides(shape), bufferFactory, initializer)
public inline fun <reified T : Any> auto( public inline fun <reified T : Any> auto(
shape: IntArray, shape: IntArray,

View File

@ -46,7 +46,11 @@ private inline class Buffer1DWrapper<T>(val buffer: Buffer<T>) : Structure1D<T>
* Represent a [NDStructure] as [Structure1D]. Throw error in case of dimension mismatch * Represent a [NDStructure] as [Structure1D]. Throw error in case of dimension mismatch
*/ */
public fun <T> NDStructure<T>.as1D(): Structure1D<T> = if (shape.size == 1) { public fun <T> NDStructure<T>.as1D(): Structure1D<T> = if (shape.size == 1) {
if (this is NDBuffer) Buffer1DWrapper(this.buffer) else Structure1DWrapper(this) when (this) {
is Structure1DWrapper -> this
is NDBuffer -> Buffer1DWrapper(this.buffer)
else -> Structure1DWrapper(this)
}
} else } else
error("Can't create 1d-structure from ${shape.size}d-structure") error("Can't create 1d-structure from ${shape.size}d-structure")

View File

@ -1,7 +1,5 @@
package space.kscience.kmath.nd package space.kscience.kmath.nd
import space.kscience.kmath.linear.BufferMatrix
import space.kscience.kmath.linear.RealMatrixContext
import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.Buffer
import space.kscience.kmath.structures.VirtualBuffer import space.kscience.kmath.structures.VirtualBuffer
@ -54,15 +52,7 @@ public interface Structure2D<T> : NDStructure<T> {
for (j in 0 until colNum) yield(intArrayOf(i, j) to get(i, j)) for (j in 0 until colNum) yield(intArrayOf(i, j) to get(i, j))
} }
public companion object { public companion object
public inline fun real(
rows: Int,
columns: Int,
crossinline init: (i: Int, j: Int) -> Double,
): BufferMatrix<Double> = RealMatrixContext.produce(rows,columns) { i, j ->
init(i, j)
}
}
} }
/** /**

View File

@ -25,7 +25,7 @@ internal class BufferAccessor2D<T : Any>(
public fun create(mat: Structure2D<T>): MutableBuffer<T> = create { i, j -> mat[i, j] } public fun create(mat: Structure2D<T>): MutableBuffer<T> = create { i, j -> mat[i, j] }
//TODO optimize wrapper //TODO optimize wrapper
public fun MutableBuffer<T>.collect(): Structure2D<T> = NDStructure.build( public fun MutableBuffer<T>.collect(): Structure2D<T> = NDStructure.buffered(
DefaultStrides(intArrayOf(rowNum, colNum)), DefaultStrides(intArrayOf(rowNum, colNum)),
factory factory
) { (i, j) -> ) { (i, j) ->

View File

@ -10,7 +10,7 @@ import kotlin.test.assertEquals
class MatrixTest { class MatrixTest {
@Test @Test
fun testTranspose() { fun testTranspose() {
val matrix = MatrixContext.real.one(3, 3) val matrix = LinearSpace.real.one(3, 3)
val transposed = matrix.transpose() val transposed = matrix.transpose()
assertEquals(matrix, transposed) assertEquals(matrix, transposed)
} }
@ -39,7 +39,7 @@ class MatrixTest {
infix fun Matrix<Double>.pow(power: Int): Matrix<Double> { infix fun Matrix<Double>.pow(power: Int): Matrix<Double> {
var res = this var res = this
repeat(power - 1) { repeat(power - 1) {
res = RealMatrixContext.invoke { res dot this@pow } res = RealLinearSpace.invoke { res dot this@pow }
} }
return res return res
} }
@ -52,7 +52,7 @@ class MatrixTest {
val firstMatrix = NDStructure.auto(2, 3) { (i, j) -> (i + j).toDouble() }.as2D() val firstMatrix = NDStructure.auto(2, 3) { (i, j) -> (i + j).toDouble() }.as2D()
val secondMatrix = NDStructure.auto(3, 2) { (i, j) -> (i + j).toDouble() }.as2D() val secondMatrix = NDStructure.auto(3, 2) { (i, j) -> (i + j).toDouble() }.as2D()
MatrixContext.real.run { LinearSpace.real.run {
// val firstMatrix = produce(2, 3) { i, j -> (i + j).toDouble() } // val firstMatrix = produce(2, 3) { i, j -> (i + j).toDouble() }
// val secondMatrix = produce(3, 2) { i, j -> (i + j).toDouble() } // val secondMatrix = produce(3, 2) { i, j -> (i + j).toDouble() }
val result = firstMatrix dot secondMatrix val result = firstMatrix dot secondMatrix

View File

@ -7,8 +7,8 @@ class RealLUSolverTest {
@Test @Test
fun testInvertOne() { fun testInvertOne() {
val matrix = MatrixContext.real.one(2, 2) val matrix = LinearSpace.real.one(2, 2)
val inverted = MatrixContext.real.inverseWithLup(matrix) val inverted = LinearSpace.real.inverseWithLup(matrix)
assertEquals(matrix, inverted) assertEquals(matrix, inverted)
} }
@ -19,7 +19,7 @@ class RealLUSolverTest {
1.0, 3.0 1.0, 3.0
) )
MatrixContext.real.run { LinearSpace.real.run {
val lup = lup(matrix) val lup = lup(matrix)
//Check determinant //Check determinant
@ -36,7 +36,7 @@ class RealLUSolverTest {
1.0, 3.0 1.0, 3.0
) )
val inverted = MatrixContext.real.inverseWithLup(matrix) val inverted = LinearSpace.real.inverseWithLup(matrix)
val expected = Matrix.square( val expected = Matrix.square(
0.375, -0.125, 0.375, -0.125,

View File

@ -77,7 +77,7 @@ public inline class DPointWrapper<T, D : Dimension>(public val point: Point<T>)
/** /**
* Basic operations on dimension-safe matrices. Operates on [Matrix] * Basic operations on dimension-safe matrices. Operates on [Matrix]
*/ */
public inline class DMatrixContext<T : Any>(public val context: MatrixContext<T, Matrix<T>>) { public inline class DMatrixContext<T : Any>(public val context: LinearSpace<T, Matrix<T>>) {
public inline fun <reified R : Dimension, reified C : Dimension> Matrix<T>.coerce(): DMatrix<T, R, C> { public inline fun <reified R : Dimension, reified C : Dimension> Matrix<T>.coerce(): DMatrix<T, R, C> {
require(rowNum == Dimension.dim<R>().toInt()) { require(rowNum == Dimension.dim<R>().toInt()) {
"Row number mismatch: expected ${Dimension.dim<R>()} but found $rowNum" "Row number mismatch: expected ${Dimension.dim<R>()} but found $rowNum"
@ -96,14 +96,14 @@ public inline class DMatrixContext<T : Any>(public val context: MatrixContext<T,
public inline fun <reified R : Dimension, reified C : Dimension> produce(noinline initializer: (i: Int, j: Int) -> T): DMatrix<T, R, C> { public inline fun <reified R : Dimension, reified C : Dimension> produce(noinline initializer: (i: Int, j: Int) -> T): DMatrix<T, R, C> {
val rows = Dimension.dim<R>() val rows = Dimension.dim<R>()
val cols = Dimension.dim<C>() val cols = Dimension.dim<C>()
return context.produce(rows.toInt(), cols.toInt(), initializer).coerce<R, C>() return context.buildMatrix(rows.toInt(), cols.toInt(), initializer).coerce<R, C>()
} }
public inline fun <reified D : Dimension> point(noinline initializer: (Int) -> T): DPoint<T, D> { public inline fun <reified D : Dimension> point(noinline initializer: (Int) -> T): DPoint<T, D> {
val size = Dimension.dim<D>() val size = Dimension.dim<D>()
return DPoint.coerceUnsafe( return DPoint.coerceUnsafe(
context.point( context.buildVector(
size.toInt(), size.toInt(),
initializer initializer
) )
@ -136,7 +136,7 @@ public inline class DMatrixContext<T : Any>(public val context: MatrixContext<T,
context { (this@transpose as Matrix<T>).transpose() }.coerce() context { (this@transpose as Matrix<T>).transpose() }.coerce()
public companion object { public companion object {
public val real: DMatrixContext<Double> = DMatrixContext(MatrixContext.real) public val real: DMatrixContext<Double> = DMatrixContext(LinearSpace.real)
} }
} }

View File

@ -11,7 +11,7 @@ import space.kscience.kmath.operations.ScaleOperations
* *
* @author Iaroslav Postovalov * @author Iaroslav Postovalov
*/ */
public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix>, ScaleOperations<Matrix<Double>> { public object EjmlLinearSpace : LinearSpace<Double, EjmlMatrix>, ScaleOperations<Matrix<Double>> {
/** /**
* Converts this matrix to EJML one. * Converts this matrix to EJML one.
@ -19,7 +19,7 @@ public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix>, ScaleOperat
@OptIn(UnstableKMathAPI::class) @OptIn(UnstableKMathAPI::class)
public fun Matrix<Double>.toEjml(): EjmlMatrix = when (val matrix = origin) { public fun Matrix<Double>.toEjml(): EjmlMatrix = when (val matrix = origin) {
is EjmlMatrix -> matrix is EjmlMatrix -> matrix
else -> produce(rowNum, colNum) { i, j -> get(i, j) } else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
} }
/** /**
@ -30,14 +30,14 @@ public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix>, ScaleOperat
(0 until it.numRows()).forEach { row -> it[row, 0] = get(row) } (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 = override fun buildMatrix(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): EjmlMatrix =
EjmlMatrix(SimpleMatrix(rows, columns).also { EjmlMatrix(SimpleMatrix(rows, columns).also {
(0 until rows).forEach { row -> (0 until rows).forEach { row ->
(0 until columns).forEach { col -> it[row, col] = initializer(row, col) } (0 until columns).forEach { col -> it[row, col] = initializer(row, col) }
} }
}) })
override fun point(size: Int, initializer: (Int) -> Double): Point<Double> = override fun buildVector(size: Int, initializer: (Int) -> Double): Point<Double> =
EjmlVector(SimpleMatrix(size, 1).also { EjmlVector(SimpleMatrix(size, 1).also {
(0 until it.numRows()).forEach { row -> it[row, 0] = initializer(row) } (0 until it.numRows()).forEach { row -> it[row, 0] = initializer(row) }
}) })
@ -58,7 +58,7 @@ public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix>, ScaleOperat
EjmlMatrix(toEjml().origin - b.toEjml().origin) EjmlMatrix(toEjml().origin - b.toEjml().origin)
public override fun scale(a: Matrix<Double>, value: Double): EjmlMatrix = public override fun scale(a: Matrix<Double>, value: Double): EjmlMatrix =
produce(a.rowNum, a.colNum) { i, j -> a[i, j] * value } buildMatrix(a.rowNum, a.colNum) { i, j -> a[i, j] * value }
public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix = public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix =
EjmlMatrix(toEjml().origin.scale(value)) EjmlMatrix(toEjml().origin.scale(value))
@ -72,7 +72,7 @@ public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix>, ScaleOperat
* @return the solution for 'x' that is n by p. * @return the solution for 'x' that is n by p.
* @author Iaroslav Postovalov * @author Iaroslav Postovalov
*/ */
public fun EjmlMatrixContext.solve(a: Matrix<Double>, b: Matrix<Double>): EjmlMatrix = public fun EjmlLinearSpace.solve(a: Matrix<Double>, b: Matrix<Double>): EjmlMatrix =
EjmlMatrix(a.toEjml().origin.solve(b.toEjml().origin)) EjmlMatrix(a.toEjml().origin.solve(b.toEjml().origin))
/** /**
@ -83,10 +83,10 @@ public fun EjmlMatrixContext.solve(a: Matrix<Double>, b: Matrix<Double>): EjmlMa
* @return the solution for 'x' that is n by p. * @return the solution for 'x' that is n by p.
* @author Iaroslav Postovalov * @author Iaroslav Postovalov
*/ */
public fun EjmlMatrixContext.solve(a: Matrix<Double>, b: Point<Double>): EjmlVector = public fun EjmlLinearSpace.solve(a: Matrix<Double>, b: Point<Double>): EjmlVector =
EjmlVector(a.toEjml().origin.solve(b.toEjml().origin)) EjmlVector(a.toEjml().origin.solve(b.toEjml().origin))
@OptIn(UnstableKMathAPI::class) @OptIn(UnstableKMathAPI::class)
public fun EjmlMatrix.inverted(): EjmlMatrix = getFeature<InverseMatrixFeature<Double>>()!!.inverse as EjmlMatrix public fun EjmlMatrix.inverted(): EjmlMatrix = getFeature<InverseMatrixFeature<Double>>()!!.inverse as EjmlMatrix
public fun EjmlMatrixContext.inverse(matrix: Matrix<Double>): Matrix<Double> = matrix.toEjml().inverted() public fun EjmlLinearSpace.inverse(matrix: Matrix<Double>): Matrix<Double> = matrix.toEjml().inverted()

View File

@ -22,14 +22,14 @@ import kotlin.math.pow
public typealias RealMatrix = Matrix<Double> public typealias RealMatrix = Matrix<Double>
public fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double): RealMatrix = public fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum, initializer) LinearSpace.real.buildMatrix(rowNum, colNum, initializer)
public fun Array<DoubleArray>.toMatrix(): RealMatrix { public fun Array<DoubleArray>.toMatrix(): RealMatrix {
return MatrixContext.real.produce(size, this[0].size) { row, col -> this[row][col] } return LinearSpace.real.buildMatrix(size, this[0].size) { row, col -> this[row][col] }
} }
public fun Sequence<DoubleArray>.toMatrix(): RealMatrix = toList().let { public fun Sequence<DoubleArray>.toMatrix(): RealMatrix = toList().let {
MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] } LinearSpace.real.buildMatrix(it.size, it[0].size) { row, col -> it[row][col] }
} }
public fun RealMatrix.repeatStackVertical(n: Int): RealMatrix = public fun RealMatrix.repeatStackVertical(n: Int): RealMatrix =
@ -42,37 +42,37 @@ public fun RealMatrix.repeatStackVertical(n: Int): RealMatrix =
*/ */
public operator fun RealMatrix.times(double: Double): RealMatrix = public operator fun RealMatrix.times(double: Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col -> LinearSpace.real.buildMatrix(rowNum, colNum) { row, col ->
this[row, col] * double this[row, col] * double
} }
public operator fun RealMatrix.plus(double: Double): RealMatrix = public operator fun RealMatrix.plus(double: Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col -> LinearSpace.real.buildMatrix(rowNum, colNum) { row, col ->
this[row, col] + double this[row, col] + double
} }
public operator fun RealMatrix.minus(double: Double): RealMatrix = public operator fun RealMatrix.minus(double: Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col -> LinearSpace.real.buildMatrix(rowNum, colNum) { row, col ->
this[row, col] - double this[row, col] - double
} }
public operator fun RealMatrix.div(double: Double): RealMatrix = public operator fun RealMatrix.div(double: Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col -> LinearSpace.real.buildMatrix(rowNum, colNum) { row, col ->
this[row, col] / double this[row, col] / double
} }
public operator fun Double.times(matrix: RealMatrix): RealMatrix = public operator fun Double.times(matrix: RealMatrix): RealMatrix =
MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> LinearSpace.real.buildMatrix(matrix.rowNum, matrix.colNum) { row, col ->
this * matrix[row, col] this * matrix[row, col]
} }
public operator fun Double.plus(matrix: RealMatrix): RealMatrix = public operator fun Double.plus(matrix: RealMatrix): RealMatrix =
MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> LinearSpace.real.buildMatrix(matrix.rowNum, matrix.colNum) { row, col ->
this + matrix[row, col] this + matrix[row, col]
} }
public operator fun Double.minus(matrix: RealMatrix): RealMatrix = public operator fun Double.minus(matrix: RealMatrix): RealMatrix =
MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col -> LinearSpace.real.buildMatrix(matrix.rowNum, matrix.colNum) { row, col ->
this - matrix[row, col] this - matrix[row, col]
} }
@ -87,20 +87,20 @@ public operator fun Double.minus(matrix: RealMatrix): RealMatrix =
@UnstableKMathAPI @UnstableKMathAPI
public operator fun RealMatrix.times(other: RealMatrix): RealMatrix = public operator fun RealMatrix.times(other: RealMatrix): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col -> this[row, col] * other[row, col] } LinearSpace.real.buildMatrix(rowNum, colNum) { row, col -> this[row, col] * other[row, col] }
public operator fun RealMatrix.plus(other: RealMatrix): RealMatrix = public operator fun RealMatrix.plus(other: RealMatrix): RealMatrix =
MatrixContext.real.add(this, other) LinearSpace.real.add(this, other)
public operator fun RealMatrix.minus(other: RealMatrix): RealMatrix = public operator fun RealMatrix.minus(other: RealMatrix): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { row, col -> this[row, col] - other[row, col] } LinearSpace.real.buildMatrix(rowNum, colNum) { row, col -> this[row, col] - other[row, col] }
/* /*
* Operations on columns * Operations on columns
*/ */
public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer<Double>) -> Double): RealMatrix = public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer<Double>) -> Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum + 1) { row, col -> LinearSpace.real.buildMatrix(rowNum, colNum + 1) { row, col ->
if (col < colNum) if (col < colNum)
this[row, col] this[row, col]
else else
@ -108,7 +108,7 @@ public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer<Double>) -
} }
public fun RealMatrix.extractColumns(columnRange: IntRange): RealMatrix = public fun RealMatrix.extractColumns(columnRange: IntRange): RealMatrix =
MatrixContext.real.produce(rowNum, columnRange.count()) { row, col -> LinearSpace.real.buildMatrix(rowNum, columnRange.count()) { row, col ->
this[row, columnRange.first + col] this[row, columnRange.first + col]
} }
@ -141,14 +141,14 @@ public fun RealMatrix.max(): Double? = elements().map { (_, value) -> value }.ma
public fun RealMatrix.average(): Double = elements().map { (_, value) -> value }.average() public fun RealMatrix.average(): Double = elements().map { (_, value) -> value }.average()
public inline fun RealMatrix.map(crossinline transform: (Double) -> Double): RealMatrix = public inline fun RealMatrix.map(crossinline transform: (Double) -> Double): RealMatrix =
MatrixContext.real.produce(rowNum, colNum) { i, j -> LinearSpace.real.buildMatrix(rowNum, colNum) { i, j ->
transform(get(i, j)) transform(get(i, j))
} }
/** /**
* Inverse a square real matrix using LUP decomposition * Inverse a square real matrix using LUP decomposition
*/ */
public fun RealMatrix.inverseWithLup(): RealMatrix = MatrixContext.real.inverseWithLup(this) public fun RealMatrix.inverseWithLup(): RealMatrix = LinearSpace.real.inverseWithLup(this)
//extended operations //extended operations

View File

@ -1,6 +1,6 @@
package kaceince.kmath.real package kaceince.kmath.real
import space.kscience.kmath.linear.MatrixContext import space.kscience.kmath.linear.LinearSpace
import space.kscience.kmath.linear.asMatrix import space.kscience.kmath.linear.asMatrix
import space.kscience.kmath.linear.real import space.kscience.kmath.linear.real
import space.kscience.kmath.linear.transpose import space.kscience.kmath.linear.transpose
@ -32,7 +32,7 @@ internal class RealVectorTest {
val vector2 = Buffer.real(5) { 5 - it.toDouble() } val vector2 = Buffer.real(5) { 5 - it.toDouble() }
val matrix1 = vector1.asMatrix() val matrix1 = vector1.asMatrix()
val matrix2 = vector2.asMatrix().transpose() val matrix2 = vector2.asMatrix().transpose()
val product = MatrixContext.real { matrix1 dot matrix2 } val product = LinearSpace.real { matrix1 dot matrix2 }
assertEquals(5.0, product[1, 0]) assertEquals(5.0, product[1, 0])
assertEquals(6.0, product[2, 2]) assertEquals(6.0, product[2, 2])
} }