Matrix revision
This commit is contained in:
parent
037735c210
commit
7e83b080ad
@ -2,31 +2,39 @@ package scientifik.kmath.linear
|
|||||||
|
|
||||||
import scientifik.kmath.operations.Field
|
import scientifik.kmath.operations.Field
|
||||||
import scientifik.kmath.operations.RealField
|
import scientifik.kmath.operations.RealField
|
||||||
|
import scientifik.kmath.operations.Ring
|
||||||
|
import scientifik.kmath.structures.*
|
||||||
import scientifik.kmath.structures.MutableBuffer.Companion.boxing
|
import scientifik.kmath.structures.MutableBuffer.Companion.boxing
|
||||||
import scientifik.kmath.structures.MutableNDStructure
|
|
||||||
import scientifik.kmath.structures.NDStructure
|
|
||||||
import scientifik.kmath.structures.get
|
|
||||||
import scientifik.kmath.structures.mutableNdStructure
|
|
||||||
import kotlin.math.absoluteValue
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Implementation based on Apache common-maths LU-decomposition
|
* Matrix LUP decomposition
|
||||||
*/
|
*/
|
||||||
abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matrix<T, F>) {
|
interface LUPDecompositionFeature<T : Any> : DeterminantFeature<T> {
|
||||||
|
/**
|
||||||
|
* A reference to L-matrix
|
||||||
|
*/
|
||||||
|
val l: Matrix<T>
|
||||||
|
/**
|
||||||
|
* A reference to u-matrix
|
||||||
|
*/
|
||||||
|
val u: Matrix<T>
|
||||||
|
/**
|
||||||
|
* Pivoting points for each row
|
||||||
|
*/
|
||||||
|
val pivot: IntArray
|
||||||
|
/**
|
||||||
|
* Permutation matrix based on [pivot]
|
||||||
|
*/
|
||||||
|
val p: Matrix<T>
|
||||||
|
}
|
||||||
|
|
||||||
private val field get() = matrix.context.ring
|
|
||||||
/** Entries of LU decomposition. */
|
|
||||||
internal val lu: NDStructure<T>
|
|
||||||
/** Pivot permutation associated with LU decomposition. */
|
|
||||||
internal val pivot: IntArray
|
|
||||||
/** Parity of the permutation associated with the LU decomposition. */
|
|
||||||
private var even: Boolean = false
|
|
||||||
|
|
||||||
init {
|
private class LUPDecomposition<T : Comparable<T>, R : Ring<T>>(
|
||||||
val pair = calculateLU()
|
val context: R,
|
||||||
lu = pair.first
|
val lu: NDStructure<T>,
|
||||||
pivot = pair.second
|
override val pivot: IntArray,
|
||||||
}
|
private val even: Boolean
|
||||||
|
) : LUPDecompositionFeature<T> {
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Returns the matrix L of the decomposition.
|
* Returns the matrix L of the decomposition.
|
||||||
@ -34,13 +42,11 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
|
|||||||
* L is a lower-triangular matrix
|
* L is a lower-triangular matrix
|
||||||
* @return the L matrix (or null if decomposed matrix is singular)
|
* @return the L matrix (or null if decomposed matrix is singular)
|
||||||
*/
|
*/
|
||||||
val l: Matrix<out T, F> by lazy {
|
override val l: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j ->
|
||||||
matrix.context.produce(matrix.numRows, matrix.numCols) { i, j ->
|
when {
|
||||||
when {
|
j < i -> lu[i, j]
|
||||||
j < i -> lu[i, j]
|
j == i -> context.one
|
||||||
j == i -> matrix.context.ring.one
|
else -> context.zero
|
||||||
else -> matrix.context.ring.zero
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -51,12 +57,11 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
|
|||||||
* U is an upper-triangular matrix
|
* U is an upper-triangular matrix
|
||||||
* @return the U matrix (or null if decomposed matrix is singular)
|
* @return the U matrix (or null if decomposed matrix is singular)
|
||||||
*/
|
*/
|
||||||
val u: Matrix<out T, F> by lazy {
|
override val u: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j ->
|
||||||
matrix.context.produce(matrix.numRows, matrix.numCols) { i, j ->
|
if (j >= i) lu[i, j] else context.zero
|
||||||
if (j >= i) lu[i, j] else field.zero
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Returns the P rows permutation matrix.
|
* Returns the P rows permutation matrix.
|
||||||
*
|
*
|
||||||
@ -67,59 +72,64 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
|
|||||||
* @return the P rows permutation matrix (or null if decomposed matrix is singular)
|
* @return the P rows permutation matrix (or null if decomposed matrix is singular)
|
||||||
* @see .getPivot
|
* @see .getPivot
|
||||||
*/
|
*/
|
||||||
val p: Matrix<out T, F> by lazy {
|
override val p: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j ->
|
||||||
matrix.context.produce(matrix.numRows, matrix.numCols) { i, j ->
|
if (j == pivot[i]) context.one else context.zero
|
||||||
//TODO ineffective. Need sparse matrix for that
|
|
||||||
if (j == pivot[i]) field.one else field.zero
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Return the determinant of the matrix
|
* Return the determinant of the matrix
|
||||||
* @return determinant of the matrix
|
* @return determinant of the matrix
|
||||||
*/
|
*/
|
||||||
val determinant: T
|
override val determinant: T by lazy {
|
||||||
get() {
|
with(context) {
|
||||||
with(matrix.context.ring) {
|
(0 until lu.shape[0]).fold(if (even) one else -one) { value, i -> value * lu[i, i] }
|
||||||
var determinant = if (even) one else -one
|
|
||||||
for (i in 0 until matrix.numRows) {
|
|
||||||
determinant *= lu[i, i]
|
|
||||||
}
|
|
||||||
return determinant
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Implementation based on Apache common-maths LU-decomposition
|
||||||
|
*/
|
||||||
|
class LUPDecompositionBuilder<T : Comparable<T>, F : Field<T>>(val context: F, val bufferFactory: MutableBufferFactory<T> = ::boxing, val singularityCheck: (T) -> Boolean) {
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* In-place transformation for [MutableNDStructure], using given transformation for each element
|
* In-place transformation for [MutableNDStructure], using given transformation for each element
|
||||||
*/
|
*/
|
||||||
operator fun <T> MutableNDStructure<T>.set(i: Int, j: Int, value: T) {
|
private operator fun <T> MutableNDStructure<T>.set(i: Int, j: Int, value: T) {
|
||||||
this[intArrayOf(i, j)] = value
|
this[intArrayOf(i, j)] = value
|
||||||
}
|
}
|
||||||
|
|
||||||
abstract fun isSingular(value: T): Boolean
|
private fun abs(value: T) = if (value > context.zero) value else with(context) { -value }
|
||||||
|
|
||||||
private fun abs(value: T) = if (value > matrix.context.ring.zero) value else with(matrix.context.ring) { -value }
|
fun decompose(matrix: Matrix<T>): LUPDecompositionFeature<T> {
|
||||||
|
// Use existing decomposition if it is provided by matrix
|
||||||
|
matrix.features.find { it is LUPDecompositionFeature<*> }?.let {
|
||||||
|
@Suppress("UNCHECKED_CAST")
|
||||||
|
return it as LUPDecompositionFeature<T>
|
||||||
|
}
|
||||||
|
|
||||||
private fun calculateLU(): Pair<NDStructure<T>, IntArray> {
|
if (matrix.rowNum != matrix.colNum) {
|
||||||
if (matrix.numRows != matrix.numCols) {
|
|
||||||
error("LU decomposition supports only square matrices")
|
error("LU decomposition supports only square matrices")
|
||||||
}
|
}
|
||||||
|
|
||||||
val m = matrix.numCols
|
val m = matrix.colNum
|
||||||
val pivot = IntArray(matrix.numRows)
|
val pivot = IntArray(matrix.rowNum)
|
||||||
//TODO fix performance
|
//TODO replace by custom optimized 2d structure
|
||||||
val lu: MutableNDStructure<T> = mutableNdStructure(
|
val lu: MutableNDStructure<T> = mutableNdStructure(
|
||||||
intArrayOf(matrix.numRows, matrix.numCols),
|
intArrayOf(matrix.rowNum, matrix.colNum),
|
||||||
::boxing
|
bufferFactory
|
||||||
) { index: IntArray -> matrix[index[0], index[1]] }
|
) { index: IntArray -> matrix[index[0], index[1]] }
|
||||||
|
|
||||||
|
|
||||||
with(matrix.context.ring) {
|
with(context) {
|
||||||
// Initialize permutation array and parity
|
// Initialize permutation array and parity
|
||||||
for (row in 0 until m) {
|
for (row in 0 until m) {
|
||||||
pivot[row] = row
|
pivot[row] = row
|
||||||
}
|
}
|
||||||
even = true
|
var even = true
|
||||||
|
|
||||||
// Loop over columns
|
// Loop over columns
|
||||||
for (col in 0 until m) {
|
for (col in 0 until m) {
|
||||||
@ -139,22 +149,18 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
|
|||||||
for (i in 0 until col) {
|
for (i in 0 until col) {
|
||||||
sum -= lu[row, i] * lu[i, col]
|
sum -= lu[row, i] * lu[i, col]
|
||||||
}
|
}
|
||||||
//luRow[col] = sum
|
|
||||||
lu[row, col] = sum
|
lu[row, col] = sum
|
||||||
|
|
||||||
abs(sum)
|
abs(sum)
|
||||||
} ?: col
|
} ?: col
|
||||||
|
|
||||||
// Singularity check
|
// Singularity check
|
||||||
if (isSingular(lu[max, col])) {
|
if (singularityCheck(lu[max, col])) {
|
||||||
error("Singular matrix")
|
error("Singular matrix")
|
||||||
}
|
}
|
||||||
|
|
||||||
// Pivot if necessary
|
// Pivot if necessary
|
||||||
if (max != col) {
|
if (max != col) {
|
||||||
//var tmp = zero
|
|
||||||
//val luMax = lu[max]
|
|
||||||
//val luCol = lu[col]
|
|
||||||
for (i in 0 until m) {
|
for (i in 0 until m) {
|
||||||
lu[max, i] = lu[col, i]
|
lu[max, i] = lu[col, i]
|
||||||
lu[col, i] = lu[max, i]
|
lu[col, i] = lu[max, i]
|
||||||
@ -169,86 +175,70 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
|
|||||||
val luDiag = lu[col, col]
|
val luDiag = lu[col, col]
|
||||||
for (row in col + 1 until m) {
|
for (row in col + 1 until m) {
|
||||||
lu[row, col] = lu[row, col] / luDiag
|
lu[row, col] = lu[row, col] / luDiag
|
||||||
// lu[row, col] /= luDiag
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
return LUPDecomposition(context, lu, pivot, even)
|
||||||
}
|
}
|
||||||
return Pair(lu, pivot)
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Returns the pivot permutation vector.
|
|
||||||
* @return the pivot permutation vector
|
|
||||||
* @see .getP
|
|
||||||
*/
|
|
||||||
fun getPivot(): IntArray = pivot.copyOf()
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
class RealLUDecomposition(matrix: RealMatrix, private val singularityThreshold: Double = DEFAULT_TOO_SMALL) :
|
|
||||||
LUDecomposition<Double, RealField>(matrix) {
|
|
||||||
override fun isSingular(value: Double): Boolean {
|
|
||||||
return value.absoluteValue < singularityThreshold
|
|
||||||
}
|
}
|
||||||
|
|
||||||
companion object {
|
companion object {
|
||||||
/** Default bound to determine effective singularity in LU decomposition. */
|
val real: LUPDecompositionBuilder<Double, RealField> = LUPDecompositionBuilder(RealField) { it < 1e-11 }
|
||||||
private const val DEFAULT_TOO_SMALL = 1e-11
|
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
/** Specialized solver. */
|
//class LUSolver<T : Comparable<T>, F : Field<T>>(val singularityCheck: (T) -> Boolean) : LinearSolver<T, F> {
|
||||||
object RealLUSolver : LinearSolver<Double, RealField> {
|
//
|
||||||
|
//
|
||||||
fun decompose(mat: Matrix<Double, RealField>, threshold: Double = 1e-11): RealLUDecomposition =
|
// override fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
|
||||||
RealLUDecomposition(mat, threshold)
|
// val decomposition = LUPDecompositionBuilder(ring, singularityCheck).decompose(a)
|
||||||
|
//
|
||||||
override fun solve(a: RealMatrix, b: RealMatrix): RealMatrix {
|
// if (b.rowNum != a.colNum) {
|
||||||
val decomposition = decompose(a, a.context.ring.zero)
|
// error("Matrix dimension mismatch expected ${a.rowNum}, but got ${b.colNum}")
|
||||||
|
// }
|
||||||
if (b.numRows != a.numCols) {
|
//
|
||||||
error("Matrix dimension mismatch expected ${a.numRows}, but got ${b.numCols}")
|
//
|
||||||
}
|
//// val bp = Array(a.rowNum) { Array<T>(b.colNum){ring.zero} }
|
||||||
|
//// for (row in 0 until a.rowNum) {
|
||||||
// Apply permutations to b
|
//// val bpRow = bp[row]
|
||||||
val bp = Array(a.numRows) { DoubleArray(b.numCols) }
|
//// val pRow = decomposition.pivot[row]
|
||||||
for (row in 0 until a.numRows) {
|
//// for (col in 0 until b.colNum) {
|
||||||
val bpRow = bp[row]
|
//// bpRow[col] = b[pRow, col]
|
||||||
val pRow = decomposition.pivot[row]
|
//// }
|
||||||
for (col in 0 until b.numCols) {
|
//// }
|
||||||
bpRow[col] = b[pRow, col]
|
//
|
||||||
}
|
// // Apply permutations to b
|
||||||
}
|
// val bp = produce(a.rowNum, a.colNum) { i, j -> b[decomposition.pivot[i], j] }
|
||||||
|
//
|
||||||
// Solve LY = b
|
// // Solve LY = b
|
||||||
for (col in 0 until a.numRows) {
|
// for (col in 0 until a.rowNum) {
|
||||||
val bpCol = bp[col]
|
// val bpCol = bp[col]
|
||||||
for (i in col + 1 until a.numRows) {
|
// for (i in col + 1 until a.rowNum) {
|
||||||
val bpI = bp[i]
|
// val bpI = bp[i]
|
||||||
val luICol = decomposition.lu[i, col]
|
// val luICol = decomposition.lu[i, col]
|
||||||
for (j in 0 until b.numCols) {
|
// for (j in 0 until b.colNum) {
|
||||||
bpI[j] -= bpCol[j] * luICol
|
// bpI[j] -= bpCol[j] * luICol
|
||||||
}
|
// }
|
||||||
}
|
// }
|
||||||
}
|
// }
|
||||||
|
//
|
||||||
// Solve UX = Y
|
// // Solve UX = Y
|
||||||
for (col in a.numRows - 1 downTo 0) {
|
// for (col in a.rowNum - 1 downTo 0) {
|
||||||
val bpCol = bp[col]
|
// val bpCol = bp[col]
|
||||||
val luDiag = decomposition.lu[col, col]
|
// val luDiag = decomposition.lu[col, col]
|
||||||
for (j in 0 until b.numCols) {
|
// for (j in 0 until b.colNum) {
|
||||||
bpCol[j] /= luDiag
|
// bpCol[j] /= luDiag
|
||||||
}
|
// }
|
||||||
for (i in 0 until col) {
|
// for (i in 0 until col) {
|
||||||
val bpI = bp[i]
|
// val bpI = bp[i]
|
||||||
val luICol = decomposition.lu[i, col]
|
// val luICol = decomposition.lu[i, col]
|
||||||
for (j in 0 until b.numCols) {
|
// for (j in 0 until b.colNum) {
|
||||||
bpI[j] -= bpCol[j] * luICol
|
// bpI[j] -= bpCol[j] * luICol
|
||||||
}
|
// }
|
||||||
}
|
// }
|
||||||
}
|
// }
|
||||||
|
//
|
||||||
return a.context.produce(a.numRows, a.numCols) { i, j -> bp[i][j] }
|
// return produce(a.rowNum, a.colNum) { i, j -> bp[i][j] }
|
||||||
}
|
// }
|
||||||
}
|
//}
|
||||||
|
@ -4,18 +4,26 @@ import scientifik.kmath.operations.Field
|
|||||||
import scientifik.kmath.operations.Norm
|
import scientifik.kmath.operations.Norm
|
||||||
import scientifik.kmath.operations.RealField
|
import scientifik.kmath.operations.RealField
|
||||||
import scientifik.kmath.operations.Ring
|
import scientifik.kmath.operations.Ring
|
||||||
import scientifik.kmath.structures.Buffer.Companion.boxing
|
|
||||||
import scientifik.kmath.structures.asSequence
|
import scientifik.kmath.structures.asSequence
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors
|
* A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors
|
||||||
*/
|
*/
|
||||||
interface LinearSolver<T : Any, F : Field<T>> {
|
interface LinearSolver<T : Any, R : Ring<T>> : MatrixContext<T, R> {
|
||||||
fun solve(a: Matrix<T, F>, b: Matrix<T, F>): Matrix<T, F>
|
/**
|
||||||
fun solve(a: Matrix<T, F>, b: Vector<T, F>): Vector<T, F> = solve(a, b.toMatrix()).toVector()
|
* Convert matrix to vector if it is possible
|
||||||
fun inverse(a: Matrix<T, F>): Matrix<T, F> = solve(a, a.context.one)
|
*/
|
||||||
|
fun Matrix<T>.toVector(): Point<T> =
|
||||||
|
if (this.colNum == 1) {
|
||||||
|
point(rowNum){ get(it, 0) }
|
||||||
|
} else error("Can't convert matrix with more than one column to vector")
|
||||||
|
|
||||||
|
fun Point<T>.toMatrix(): Matrix<T> = produce(size, 1) { i, _ -> get(i) }
|
||||||
|
|
||||||
|
fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T>
|
||||||
|
fun solve(a: Matrix<T>, b: Point<T>): Point<T> = solve(a, b.toMatrix()).toVector()
|
||||||
|
fun inverse(a: Matrix<T>): Matrix<T> = solve(a, one(a.rowNum, a.colNum))
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
@ -26,39 +34,10 @@ fun <T : Any> Array<T>.toVector(field: Field<T>) = Vector.generic(size, field) {
|
|||||||
fun DoubleArray.toVector() = Vector.real(this.size) { this[it] }
|
fun DoubleArray.toVector() = Vector.real(this.size) { this[it] }
|
||||||
fun List<Double>.toVector() = Vector.real(this.size) { this[it] }
|
fun List<Double>.toVector() = Vector.real(this.size) { this[it] }
|
||||||
|
|
||||||
/**
|
object VectorL2Norm : Norm<Point<out Number>, Double> {
|
||||||
* Convert matrix to vector if it is possible
|
override fun norm(arg: Point<out Number>): Double =
|
||||||
*/
|
kotlin.math.sqrt(arg.asSequence().sumByDouble { it.toDouble() })
|
||||||
fun <T : Any, F : Ring<T>> Matrix<T, F>.toVector(): Vector<T, F> =
|
|
||||||
if (this.numCols == 1) {
|
|
||||||
// if (this is ArrayMatrix) {
|
|
||||||
// //Reuse existing underlying array
|
|
||||||
// ArrayVector(ArrayVectorSpace(rows, context.field, context.ndFactory), array)
|
|
||||||
// } else {
|
|
||||||
// //Generic vector
|
|
||||||
// vector(rows, context.field) { get(it, 0) }
|
|
||||||
// }
|
|
||||||
Vector.generic(numRows, context.ring) { get(it, 0) }
|
|
||||||
} else error("Can't convert matrix with more than one column to vector")
|
|
||||||
|
|
||||||
fun <T : Any, R : Ring<T>> Vector<T, R>.toMatrix(): Matrix<T, R> {
|
|
||||||
// val context = StructureMatrixContext(size, 1, context.space)
|
|
||||||
//
|
|
||||||
// return if (this is ArrayVector) {
|
|
||||||
// //Reuse existing underlying array
|
|
||||||
// StructureMatrix(context,this.buffer)
|
|
||||||
// } else {
|
|
||||||
// //Generic vector
|
|
||||||
// matrix(size, 1, context.field) { i, j -> get(i) }
|
|
||||||
// }
|
|
||||||
//return Matrix.of(size, 1, context.space) { i, _ -> get(i) }
|
|
||||||
return StructureMatrixSpace(size, 1, context.space, ::boxing).produce { i, _ -> get(i) }
|
|
||||||
}
|
|
||||||
|
|
||||||
object VectorL2Norm : Norm<Vector<out Number, *>, Double> {
|
|
||||||
override fun norm(arg: Vector<out Number, *>): Double =
|
|
||||||
kotlin.math.sqrt(arg.asSequence().sumByDouble { it.toDouble() })
|
|
||||||
}
|
}
|
||||||
|
|
||||||
typealias RealVector = Vector<Double, RealField>
|
typealias RealVector = Vector<Double, RealField>
|
||||||
typealias RealMatrix = Matrix<Double, RealField>
|
typealias RealMatrix = Matrix<Double>
|
||||||
|
@ -2,64 +2,91 @@ package scientifik.kmath.linear
|
|||||||
|
|
||||||
import scientifik.kmath.operations.RealField
|
import scientifik.kmath.operations.RealField
|
||||||
import scientifik.kmath.operations.Ring
|
import scientifik.kmath.operations.Ring
|
||||||
import scientifik.kmath.operations.Space
|
|
||||||
import scientifik.kmath.operations.SpaceElement
|
|
||||||
import scientifik.kmath.structures.*
|
import scientifik.kmath.structures.*
|
||||||
import scientifik.kmath.structures.Buffer.Companion.DoubleBufferFactory
|
import scientifik.kmath.structures.Buffer.Companion.DoubleBufferFactory
|
||||||
import scientifik.kmath.structures.Buffer.Companion.boxing
|
import scientifik.kmath.structures.Buffer.Companion.boxing
|
||||||
|
|
||||||
|
|
||||||
interface MatrixSpace<T : Any, R : Ring<T>> : Space<Matrix<T, R>> {
|
interface MatrixContext<T : Any, R : Ring<T>> {
|
||||||
/**
|
/**
|
||||||
* The ring context for matrix elements
|
* The ring context for matrix elements
|
||||||
*/
|
*/
|
||||||
val ring: R
|
val ring: R
|
||||||
|
|
||||||
val rowNum: Int
|
|
||||||
val colNum: Int
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Produce a matrix with this context and given dimensions
|
* Produce a matrix with this context and given dimensions
|
||||||
*/
|
*/
|
||||||
fun produce(rows: Int = rowNum, columns: Int = colNum, initializer: (i: Int, j: Int) -> T): Matrix<T, R>
|
fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T>
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Produce a point compatible with matrix space
|
* Produce a point compatible with matrix space
|
||||||
*/
|
*/
|
||||||
fun point(size: Int, initializer: (Int) -> T): Point<T>
|
fun point(size: Int, initializer: (Int) -> T): Point<T>
|
||||||
|
|
||||||
override val zero: Matrix<T, R> get() = produce { _, _ -> ring.zero }
|
fun scale(a: Matrix<T>, k: Number): Matrix<T> {
|
||||||
|
//TODO create a special wrapper class for scaled matrices
|
||||||
|
return produce(a.rowNum, a.colNum) { i, j -> ring.run { a[i, j] * k } }
|
||||||
|
}
|
||||||
|
|
||||||
val one get() = produce { i, j -> if (i == j) ring.one else ring.zero }
|
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(ring) {
|
||||||
|
row.asSequence().zip(column.asSequence(), ::multiply).sum()
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
override fun add(a: Matrix<T, R>, b: Matrix<T, R>): Matrix<T, R> =
|
infix fun Matrix<T>.dot(vector: Point<T>): Point<T> {
|
||||||
produce(rowNum, colNum) { i, j -> ring.run { a[i, j] + b[i, j] } }
|
//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(ring) {
|
||||||
|
row.asSequence().zip(vector.asSequence(), ::multiply).sum()
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
operator fun Matrix<T>.unaryMinus() =
|
||||||
|
produce(rowNum, colNum) { i, j -> ring.run { -get(i, j) } }
|
||||||
|
|
||||||
|
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 -> ring.run { get(i, j) + b[i, j] } }
|
||||||
|
}
|
||||||
|
|
||||||
|
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 -> ring.run { get(i, j) + b[i, j] } }
|
||||||
|
}
|
||||||
|
|
||||||
|
operator fun Matrix<T>.times(number: Number): Matrix<T> =
|
||||||
|
produce(rowNum, colNum) { i, j -> ring.run { get(i, j) * number } }
|
||||||
|
|
||||||
|
operator fun Number.times(m: Matrix<T>): Matrix<T> = m * this
|
||||||
|
|
||||||
override fun multiply(a: Matrix<T, R>, k: Number): Matrix<T, R> =
|
|
||||||
produce(rowNum, colNum) { i, j -> ring.run { a[i, j] * k } }
|
|
||||||
|
|
||||||
companion object {
|
companion object {
|
||||||
/**
|
/**
|
||||||
* Non-boxing double matrix
|
* Non-boxing double matrix
|
||||||
*/
|
*/
|
||||||
fun real(rows: Int, columns: Int): MatrixSpace<Double, RealField> =
|
val real: MatrixContext<Double, RealField> = StructureMatrixContext(RealField, DoubleBufferFactory)
|
||||||
StructureMatrixSpace(rows, columns, RealField, DoubleBufferFactory)
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* A structured matrix with custom buffer
|
* A structured matrix with custom buffer
|
||||||
*/
|
*/
|
||||||
fun <T : Any, R : Ring<T>> buffered(
|
fun <T : Any, R : Ring<T>> buffered(ring: R, bufferFactory: BufferFactory<T> = ::boxing): MatrixContext<T, R> =
|
||||||
rows: Int,
|
StructureMatrixContext(ring, bufferFactory)
|
||||||
columns: Int,
|
|
||||||
ring: R,
|
|
||||||
bufferFactory: BufferFactory<T> = ::boxing
|
|
||||||
): MatrixSpace<T, R> = StructureMatrixSpace(rows, columns, ring, bufferFactory)
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Automatic buffered matrix, unboxed if it is possible
|
* Automatic buffered matrix, unboxed if it is possible
|
||||||
*/
|
*/
|
||||||
inline fun <reified T : Any, R : Ring<T>> auto(rows: Int, columns: Int, ring: R): MatrixSpace<T, R> =
|
inline fun <reified T : Any, R : Ring<T>> auto(ring: R): MatrixContext<T, R> =
|
||||||
buffered(rows, columns, ring, Buffer.Companion::auto)
|
buffered(ring, Buffer.Companion::auto)
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -69,108 +96,102 @@ interface MatrixSpace<T : Any, R : Ring<T>> : Space<Matrix<T, R>> {
|
|||||||
*/
|
*/
|
||||||
interface MatrixFeature
|
interface MatrixFeature
|
||||||
|
|
||||||
|
object DiagonalFeature : MatrixFeature
|
||||||
|
|
||||||
|
object ZeroFeature : MatrixFeature
|
||||||
|
|
||||||
|
object UnitFeature : MatrixFeature
|
||||||
|
|
||||||
|
interface InverseMatrixFeature<T : Any> : MatrixFeature {
|
||||||
|
val inverse: Matrix<T>
|
||||||
|
}
|
||||||
|
|
||||||
|
interface DeterminantFeature<T : Any> : MatrixFeature {
|
||||||
|
val determinant: T
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Specialized 2-d structure
|
* Specialized 2-d structure
|
||||||
*/
|
*/
|
||||||
interface Matrix<T : Any, R : Ring<T>> : NDStructure<T>, SpaceElement<Matrix<T, R>, Matrix<T, R>, MatrixSpace<T, R>> {
|
interface Matrix<T : Any> : NDStructure<T> {
|
||||||
|
val rowNum: Int
|
||||||
|
val colNum: Int
|
||||||
|
|
||||||
|
val features: Set<MatrixFeature>
|
||||||
|
|
||||||
operator fun get(i: Int, j: Int): T
|
operator fun get(i: Int, j: Int): T
|
||||||
|
|
||||||
override fun get(index: IntArray): T = get(index[0], index[1])
|
override fun get(index: IntArray): T = get(index[0], index[1])
|
||||||
|
|
||||||
val numRows get() = context.rowNum
|
override val shape: IntArray get() = intArrayOf(rowNum, colNum)
|
||||||
val numCols get() = context.colNum
|
|
||||||
|
|
||||||
//TODO replace by lazy buffers
|
|
||||||
val rows: Point<Point<T>>
|
val rows: Point<Point<T>>
|
||||||
get() = ListBuffer((0 until numRows).map { i ->
|
get() = VirtualBuffer(rowNum) { i ->
|
||||||
context.point(numCols) { j -> get(i, j) }
|
VirtualBuffer(colNum) { j -> get(i, j) }
|
||||||
})
|
}
|
||||||
|
|
||||||
val columns: Point<Point<T>>
|
val columns: Point<Point<T>>
|
||||||
get() = ListBuffer((0 until numCols).map { j ->
|
get() = VirtualBuffer(colNum) { j ->
|
||||||
context.point(numRows) { i -> get(i, j) }
|
VirtualBuffer(rowNum) { i -> get(i, j) }
|
||||||
})
|
}
|
||||||
|
|
||||||
val features: Set<MatrixFeature>
|
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))
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
companion object {
|
companion object {
|
||||||
fun real(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
|
fun real(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
|
||||||
MatrixSpace.real(rows, columns).produce(rows, columns, initializer)
|
MatrixContext.real.produce(rows, columns, initializer)
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
infix fun <T : Any, R : Ring<T>> Matrix<T, R>.dot(other: Matrix<T, R>): Matrix<T, R> {
|
/**
|
||||||
//TODO add typed error
|
* Diagonal matrix of ones. The matrix is virtual no actual matrix is created
|
||||||
if (this.numCols != other.numRows) error("Matrix dot operation dimension mismatch: ($numRows, $numCols) x (${other.numRows}, ${other.numCols})")
|
*/
|
||||||
return context.produce(numRows, other.numCols) { i, j ->
|
fun <T : Any, R : Ring<T>> MatrixContext<T, R>.one(rows: Int, columns: Int): Matrix<T> {
|
||||||
val row = rows[i]
|
return object : Matrix<T> {
|
||||||
val column = other.columns[j]
|
override val rowNum: Int get() = rows
|
||||||
with(context.ring) {
|
override val colNum: Int get() = columns
|
||||||
row.asSequence().zip(column.asSequence(), ::multiply).sum()
|
override val features: Set<MatrixFeature> get() = setOf(DiagonalFeature, UnitFeature)
|
||||||
}
|
override fun get(i: Int, j: Int): T = if (i == j) ring.one else ring.zero
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
infix fun <T : Any, R : Ring<T>> Matrix<T, R>.dot(vector: Point<T>): Point<T> {
|
/**
|
||||||
//TODO add typed error
|
* A virtual matrix of zeroes
|
||||||
if (this.numCols != vector.size) error("Matrix dot vector operation dimension mismatch: ($numRows, $numCols) x (${vector.size})")
|
*/
|
||||||
return context.point(numRows) { i ->
|
fun <T : Any, R : Ring<T>> MatrixContext<T, R>.zero(rows: Int, columns: Int): Matrix<T> {
|
||||||
val row = rows[i]
|
return object : Matrix<T> {
|
||||||
with(context.ring) {
|
override val rowNum: Int get() = rows
|
||||||
row.asSequence().zip(vector.asSequence(), ::multiply).sum()
|
override val colNum: Int get() = columns
|
||||||
}
|
override val features: Set<MatrixFeature> get() = setOf(ZeroFeature)
|
||||||
|
override fun get(i: Int, j: Int): T = ring.zero
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
data class StructureMatrixSpace<T : Any, R : Ring<T>>(
|
|
||||||
override val rowNum: Int,
|
|
||||||
override val colNum: Int,
|
|
||||||
override val ring: R,
|
|
||||||
private val bufferFactory: BufferFactory<T>
|
|
||||||
) : MatrixSpace<T, R> {
|
|
||||||
|
|
||||||
val shape: IntArray = intArrayOf(rowNum, colNum)
|
inline class TransposedMatrix<T : Any>(val original: Matrix<T>) : Matrix<T> {
|
||||||
|
override val rowNum: Int get() = original.colNum
|
||||||
|
override val colNum: Int get() = original.rowNum
|
||||||
|
override val features: Set<MatrixFeature> get() = emptySet() //TODO retain some features
|
||||||
|
|
||||||
private val strides = DefaultStrides(shape)
|
override fun get(i: Int, j: Int): T = original[j, i]
|
||||||
|
|
||||||
override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T, R> {
|
override fun elements(): Sequence<Pair<IntArray, T>> =
|
||||||
return if (rows == rowNum && columns == colNum) {
|
original.elements().map { (key, value) -> intArrayOf(key[1], key[0]) to value }
|
||||||
val structure = ndStructure(strides, bufferFactory) { initializer(it[0], it[1]) }
|
|
||||||
StructureMatrix(this, structure)
|
|
||||||
} else {
|
|
||||||
val context = StructureMatrixSpace(rows, columns, ring, bufferFactory)
|
|
||||||
val structure = ndStructure(context.strides, bufferFactory) { initializer(it[0], it[1]) }
|
|
||||||
StructureMatrix(context, structure)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
override fun point(size: Int, initializer: (Int) -> T): Point<T> = bufferFactory(size, initializer)
|
|
||||||
}
|
}
|
||||||
|
|
||||||
data class StructureMatrix<T : Any, R : Ring<T>>(
|
/**
|
||||||
override val context: StructureMatrixSpace<T, R>,
|
* Create a virtual transposed matrix without copying anything. `A.transpose().transpose() === A`
|
||||||
val structure: NDStructure<T>,
|
*/
|
||||||
override val features: Set<MatrixFeature> = emptySet()
|
fun <T : Any, R : Ring<T>> Matrix<T>.transpose(): Matrix<T> {
|
||||||
) : Matrix<T, R> {
|
return if (this is TransposedMatrix) {
|
||||||
init {
|
original
|
||||||
if (structure.shape.size != 2 || structure.shape[0] != context.rowNum || structure.shape[1] != context.colNum) {
|
} else {
|
||||||
error("Dimension mismatch for structure, (${context.rowNum}, ${context.colNum}) expected, but ${structure.shape} found")
|
TransposedMatrix(this)
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
override fun unwrap(): Matrix<T, R> = this
|
|
||||||
|
|
||||||
override fun Matrix<T, R>.wrap(): Matrix<T, R> = this
|
|
||||||
|
|
||||||
override val shape: IntArray get() = structure.shape
|
|
||||||
|
|
||||||
override fun get(index: IntArray): T = structure[index]
|
|
||||||
|
|
||||||
override fun get(i: Int, j: Int): T = structure[i, j]
|
|
||||||
|
|
||||||
override fun elements(): Sequence<Pair<IntArray, T>> = structure.elements()
|
|
||||||
}
|
}
|
||||||
|
|
||||||
//TODO produce transposed matrix via reference without creating new space and structure
|
|
||||||
fun <T : Any, R : Ring<T>> Matrix<T, R>.transpose(): Matrix<T, R> =
|
|
||||||
context.produce(numCols, numRows) { i, j -> get(j, i) }
|
|
@ -0,0 +1,50 @@
|
|||||||
|
package scientifik.kmath.linear
|
||||||
|
|
||||||
|
import scientifik.kmath.operations.Ring
|
||||||
|
import scientifik.kmath.structures.BufferFactory
|
||||||
|
import scientifik.kmath.structures.NDStructure
|
||||||
|
import scientifik.kmath.structures.get
|
||||||
|
import scientifik.kmath.structures.ndStructure
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Basic implementation of Matrix space based on [NDStructure]
|
||||||
|
*/
|
||||||
|
class StructureMatrixContext<T : Any, R : Ring<T>>(
|
||||||
|
override val ring: R,
|
||||||
|
private val bufferFactory: BufferFactory<T>
|
||||||
|
) : MatrixContext<T, R> {
|
||||||
|
|
||||||
|
override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T> {
|
||||||
|
val structure =
|
||||||
|
ndStructure(intArrayOf(rows, columns), bufferFactory) { index -> initializer(index[0], index[1]) }
|
||||||
|
return StructureMatrix(structure)
|
||||||
|
}
|
||||||
|
|
||||||
|
override fun point(size: Int, initializer: (Int) -> T): Point<T> = bufferFactory(size, initializer)
|
||||||
|
}
|
||||||
|
|
||||||
|
data class StructureMatrix<T : Any>(
|
||||||
|
val structure: NDStructure<T>,
|
||||||
|
override val features: Set<MatrixFeature> = emptySet()
|
||||||
|
) : Matrix<T> {
|
||||||
|
|
||||||
|
init {
|
||||||
|
if (structure.shape.size != 2) {
|
||||||
|
error("Dimension mismatch for matrix structure")
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
override val rowNum: Int
|
||||||
|
get() = structure.shape[0]
|
||||||
|
override val colNum: Int
|
||||||
|
get() = structure.shape[1]
|
||||||
|
|
||||||
|
|
||||||
|
override val shape: IntArray get() = structure.shape
|
||||||
|
|
||||||
|
override fun get(index: IntArray): T = structure[index]
|
||||||
|
|
||||||
|
override fun get(i: Int, j: Int): T = structure[i, j]
|
||||||
|
|
||||||
|
override fun elements(): Sequence<Pair<IntArray, T>> = structure.elements()
|
||||||
|
}
|
@ -0,0 +1,10 @@
|
|||||||
|
package scientifik.kmath.linear
|
||||||
|
|
||||||
|
class VirtualMatrix<T : Any>(
|
||||||
|
override val rowNum: Int,
|
||||||
|
override val colNum: Int,
|
||||||
|
override val features: Set<MatrixFeature> = emptySet(),
|
||||||
|
val generator: (i: Int, j: Int) -> T
|
||||||
|
) : Matrix<T> {
|
||||||
|
override fun get(i: Int, j: Int): T = generator(i, j)
|
||||||
|
}
|
@ -33,8 +33,8 @@ interface Space<T> {
|
|||||||
operator fun T.times(k: Number) = multiply(this, k.toDouble())
|
operator fun T.times(k: Number) = multiply(this, k.toDouble())
|
||||||
operator fun T.div(k: Number) = multiply(this, 1.0 / k.toDouble())
|
operator fun T.div(k: Number) = multiply(this, 1.0 / k.toDouble())
|
||||||
operator fun Number.times(b: T) = b * this
|
operator fun Number.times(b: T) = b * this
|
||||||
fun Iterable<T>.sum(): T = fold(zero) { left, right -> add(left,right) }
|
|
||||||
|
|
||||||
|
fun Iterable<T>.sum(): T = fold(zero) { left, right -> add(left,right) }
|
||||||
fun Sequence<T>.sum(): T = fold(zero) { left, right -> add(left, right) }
|
fun Sequence<T>.sum(): T = fold(zero) { left, right -> add(left, right) }
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -95,6 +95,8 @@ inline class ListBuffer<T>(private val list: List<T>) : Buffer<T> {
|
|||||||
override fun iterator(): Iterator<T> = list.iterator()
|
override fun iterator(): Iterator<T> = list.iterator()
|
||||||
}
|
}
|
||||||
|
|
||||||
|
fun <T> List<T>.asBuffer() = ListBuffer(this)
|
||||||
|
|
||||||
inline class MutableListBuffer<T>(private val list: MutableList<T>) : MutableBuffer<T> {
|
inline class MutableListBuffer<T>(private val list: MutableList<T>) : MutableBuffer<T> {
|
||||||
|
|
||||||
override val size: Int
|
override val size: Int
|
||||||
@ -191,6 +193,25 @@ inline class ReadOnlyBuffer<T>(private val buffer: MutableBuffer<T>) : Buffer<T>
|
|||||||
override fun iterator(): Iterator<T> = buffer.iterator()
|
override fun iterator(): Iterator<T> = buffer.iterator()
|
||||||
}
|
}
|
||||||
|
|
||||||
|
/**
|
||||||
|
* A buffer with content calculated on-demand. The calculated contect is not stored, so it is recalculated on each call.
|
||||||
|
* Useful when one needs single element from the buffer.
|
||||||
|
*/
|
||||||
|
class VirtualBuffer<T>(override val size: Int, private val generator: (Int) -> T) : Buffer<T> {
|
||||||
|
override fun get(index: Int): T = generator(index)
|
||||||
|
|
||||||
|
override fun iterator(): Iterator<T> = (0 until size).asSequence().map(generator).iterator()
|
||||||
|
|
||||||
|
override fun contentEquals(other: Buffer<*>): Boolean {
|
||||||
|
return if (other is VirtualBuffer) {
|
||||||
|
this.size == other.size && this.generator == other.generator
|
||||||
|
} else {
|
||||||
|
super.contentEquals(other)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Convert this buffer to read-only buffer
|
* Convert this buffer to read-only buffer
|
||||||
*/
|
*/
|
||||||
|
@ -22,7 +22,7 @@ class MatrixTest {
|
|||||||
|
|
||||||
@Test
|
@Test
|
||||||
fun testTranspose() {
|
fun testTranspose() {
|
||||||
val matrix = MatrixSpace.real(3, 3).one
|
val matrix = MatrixContext.real(3, 3).one
|
||||||
val transposed = matrix.transpose()
|
val transposed = matrix.transpose()
|
||||||
assertEquals(matrix.context, transposed.context)
|
assertEquals(matrix.context, transposed.context)
|
||||||
assertEquals((matrix as StructureMatrix).structure, (transposed as StructureMatrix).structure)
|
assertEquals((matrix as StructureMatrix).structure, (transposed as StructureMatrix).structure)
|
||||||
|
@ -6,7 +6,7 @@ import kotlin.test.assertEquals
|
|||||||
class RealLUSolverTest {
|
class RealLUSolverTest {
|
||||||
@Test
|
@Test
|
||||||
fun testInvertOne() {
|
fun testInvertOne() {
|
||||||
val matrix = MatrixSpace.real(2, 2).one
|
val matrix = MatrixContext.real(2, 2).one
|
||||||
val inverted = RealLUSolver.inverse(matrix)
|
val inverted = RealLUSolver.inverse(matrix)
|
||||||
assertEquals(matrix, inverted)
|
assertEquals(matrix, inverted)
|
||||||
}
|
}
|
||||||
|
Loading…
Reference in New Issue
Block a user