Refactor/linear algebra #241

Merged
altavir merged 11 commits from refactor/linear-algebra into dev 2021-03-14 21:49:47 +03:00
10 changed files with 41 additions and 126 deletions
Showing only changes of commit e01f90b5e0 - Show all commits

View File

@ -9,6 +9,7 @@
- Exponential operations merged with hyperbolic functions - Exponential operations merged with hyperbolic functions
- Space is replaced by Group. Space is reserved for vector spaces. - Space is replaced by Group. Space is reserved for vector spaces.
- VectorSpace is now a vector space - VectorSpace is now a vector space
-
### Deprecated ### Deprecated

View File

@ -4,15 +4,15 @@ import space.kscience.kmath.real.*
import space.kscience.kmath.structures.RealBuffer import space.kscience.kmath.structures.RealBuffer
fun main() { fun main() {
val x0 = Vector(0.0, 0.0, 0.0) val x0 = Point(0.0, 0.0, 0.0)
val sigma = Vector(1.0, 1.0, 1.0) val sigma = Point(1.0, 1.0, 1.0)
val gaussian: (Vector<Double>) -> Double = { x -> val gaussian: (Point<Double>) -> Double = { x ->
require(x.size == x0.size) require(x.size == x0.size)
kotlin.math.exp(-((x - x0) / sigma).square().sum()) kotlin.math.exp(-((x - x0) / sigma).square().sum())
} }
fun ((Vector<Double>) -> Double).grad(x: Vector<Double>): Vector<Double> { fun ((Point<Double>) -> Double).grad(x: Point<Double>): Point<Double> {
require(x.size == x0.size) require(x.size == x0.size)
return RealBuffer(x.size) { i -> return RealBuffer(x.size) { i ->
val h = sigma[i] / 5 val h = sigma[i] / 5

View File

@ -62,7 +62,7 @@ public class CMMatrix(public val origin: RealMatrix) : Matrix<Double> {
override fun hashCode(): Int = origin.hashCode() override fun hashCode(): Int = origin.hashCode()
} }
public inline class CMVector(public val origin: RealVector) : Vector<Double> { public inline class CMVector(public val origin: RealVector) : Point<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)
@ -94,7 +94,7 @@ public object CMLinearSpace : LinearSpace<Double, RealField> {
} }
} }
public fun Vector<Double>.toCM(): CMVector = if (this is CMVector) this else { public fun Point<Double>.toCM(): CMVector = if (this is CMVector) this else {
val array = DoubleArray(size) { this[it] } val array = DoubleArray(size) { this[it] }
ArrayRealVector(array).wrap() ArrayRealVector(array).wrap()
} }
@ -102,22 +102,22 @@ public object CMLinearSpace : LinearSpace<Double, RealField> {
internal fun RealMatrix.wrap(): CMMatrix = CMMatrix(this) internal fun RealMatrix.wrap(): CMMatrix = CMMatrix(this)
internal fun RealVector.wrap(): CMVector = CMVector(this) internal fun RealVector.wrap(): CMVector = CMVector(this)
override fun buildVector(size: Int, initializer: RealField.(Int) -> Double): Vector<Double> = override fun buildVector(size: Int, initializer: RealField.(Int) -> Double): Point<Double> =
ArrayRealVector(DoubleArray(size) { RealField.initializer(it) }).wrap() ArrayRealVector(DoubleArray(size) { RealField.initializer(it) }).wrap()
override fun Matrix<Double>.plus(other: Matrix<Double>): CMMatrix = override fun Matrix<Double>.plus(other: Matrix<Double>): CMMatrix =
toCM().origin.add(other.toCM().origin).wrap() toCM().origin.add(other.toCM().origin).wrap()
override fun Vector<Double>.plus(other: Vector<Double>): CMVector = override fun Point<Double>.plus(other: Point<Double>): CMVector =
toCM().origin.add(other.toCM().origin).wrap() toCM().origin.add(other.toCM().origin).wrap()
override fun Vector<Double>.minus(other: Vector<Double>): CMVector = override fun Point<Double>.minus(other: Point<Double>): CMVector =
toCM().origin.subtract(other.toCM().origin).wrap() toCM().origin.subtract(other.toCM().origin).wrap()
public override fun Matrix<Double>.dot(other: Matrix<Double>): CMMatrix = public override fun Matrix<Double>.dot(other: Matrix<Double>): CMMatrix =
toCM().origin.multiply(other.toCM().origin).wrap() toCM().origin.multiply(other.toCM().origin).wrap()
public override fun Matrix<Double>.dot(vector: Vector<Double>): CMVector = public override fun Matrix<Double>.dot(vector: Point<Double>): CMVector =
toCM().origin.preMultiply(vector.toCM().origin).wrap() toCM().origin.preMultiply(vector.toCM().origin).wrap()
public override operator fun Matrix<Double>.minus(other: Matrix<Double>): CMMatrix = public override operator fun Matrix<Double>.minus(other: Matrix<Double>): CMMatrix =
@ -129,10 +129,10 @@ public object CMLinearSpace : LinearSpace<Double, RealField> {
override fun Double.times(m: Matrix<Double>): CMMatrix = override fun Double.times(m: Matrix<Double>): CMMatrix =
m * this m * this
override fun Vector<Double>.times(value: Double): CMVector = override fun Point<Double>.times(value: Double): CMVector =
toCM().origin.mapMultiply(value).wrap() toCM().origin.mapMultiply(value).wrap()
override fun Double.times(v: Vector<Double>): CMVector = override fun Double.times(v: Point<Double>): CMVector =
v * this v * this
} }

View File

@ -474,11 +474,6 @@ public final class space/kscience/kmath/linear/LFeature : space/kscience/kmath/l
public static final field INSTANCE Lspace/kscience/kmath/linear/LFeature; public static final field INSTANCE Lspace/kscience/kmath/linear/LFeature;
} }
public final class space/kscience/kmath/linear/LinearAlgebraKt {
public static final fun asMatrix (Lspace/kscience/kmath/structures/Buffer;)Lspace/kscience/kmath/linear/VirtualMatrix;
public static final fun asVector (Lspace/kscience/kmath/nd/Structure2D;)Lspace/kscience/kmath/structures/Buffer;
}
public abstract interface class space/kscience/kmath/linear/LinearSolver { public abstract interface class space/kscience/kmath/linear/LinearSolver {
public abstract fun inverse (Lspace/kscience/kmath/nd/Structure2D;)Lspace/kscience/kmath/nd/Structure2D; public abstract fun inverse (Lspace/kscience/kmath/nd/Structure2D;)Lspace/kscience/kmath/nd/Structure2D;
public abstract fun solve (Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/nd/Structure2D;)Lspace/kscience/kmath/nd/Structure2D; public abstract fun solve (Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/nd/Structure2D;)Lspace/kscience/kmath/nd/Structure2D;
@ -489,6 +484,11 @@ public final class space/kscience/kmath/linear/LinearSolver$DefaultImpls {
public static fun solve (Lspace/kscience/kmath/linear/LinearSolver;Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/structures/Buffer;)Lspace/kscience/kmath/structures/Buffer; public static fun solve (Lspace/kscience/kmath/linear/LinearSolver;Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/structures/Buffer;)Lspace/kscience/kmath/structures/Buffer;
} }
public final class space/kscience/kmath/linear/LinearSolverKt {
public static final fun asMatrix (Lspace/kscience/kmath/structures/Buffer;)Lspace/kscience/kmath/linear/VirtualMatrix;
public static final fun asVector (Lspace/kscience/kmath/nd/Structure2D;)Lspace/kscience/kmath/structures/Buffer;
}
public abstract interface class space/kscience/kmath/linear/LinearSpace { public abstract interface class space/kscience/kmath/linear/LinearSpace {
public static final field Companion Lspace/kscience/kmath/linear/LinearSpace$Companion; public static final field Companion Lspace/kscience/kmath/linear/LinearSpace$Companion;
public abstract fun buildMatrix (IILkotlin/jvm/functions/Function3;)Lspace/kscience/kmath/nd/Structure2D; public abstract fun buildMatrix (IILkotlin/jvm/functions/Function3;)Lspace/kscience/kmath/nd/Structure2D;

View File

@ -22,7 +22,7 @@ public class BufferLinearSpace<T : Any, A : Ring<T>>(
override fun buildMatrix(rows: Int, columns: Int, initializer: A.(i: Int, j: Int) -> T): Matrix<T> = override fun buildMatrix(rows: Int, columns: Int, initializer: A.(i: Int, j: Int) -> T): Matrix<T> =
ndRing(rows, columns).produce { (i, j) -> elementAlgebra.initializer(i, j) }.as2D() ndRing(rows, columns).produce { (i, j) -> elementAlgebra.initializer(i, j) }.as2D()
override fun buildVector(size: Int, initializer: A.(Int) -> T): Vector<T> = override fun buildVector(size: Int, initializer: A.(Int) -> T): Point<T> =
bufferFactory(size) { elementAlgebra.initializer(it) } bufferFactory(size) { elementAlgebra.initializer(it) }
override fun Matrix<T>.unaryMinus(): Matrix<T> = ndRing(rowNum, colNum).run { override fun Matrix<T>.unaryMinus(): Matrix<T> = ndRing(rowNum, colNum).run {
@ -62,7 +62,7 @@ public class BufferLinearSpace<T : Any, A : Ring<T>>(
} }
} }
override fun Matrix<T>.dot(vector: Vector<T>): Vector<T> { override fun Matrix<T>.dot(vector: Point<T>): Point<T> {
require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" } require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
return elementAlgebra { return elementAlgebra {
val rows = this@dot.rows.map { it.linearize() } val rows = this@dot.rows.map { it.linearize() }

View File

@ -1,9 +1,6 @@
package space.kscience.kmath.linear package space.kscience.kmath.linear
import space.kscience.kmath.nd.as1D import space.kscience.kmath.nd.as1D
import space.kscience.kmath.structures.Buffer
public typealias Point<T> = Buffer<T>
/** /**
* 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
@ -17,7 +14,7 @@ public interface LinearSolver<T : Any> {
/** /**
* Convert matrix to vector if it is possible * Convert matrix to vector if it is possible
*/ */
public fun <T : Any> Matrix<T>.asVector(): Vector<T> = public fun <T : Any> Matrix<T>.asVector(): Point<T> =
if (this.colNum == 1) if (this.colNum == 1)
as1D() as1D()
else else

View File

@ -14,7 +14,10 @@ import kotlin.reflect.KClass
*/ */
public typealias Matrix<T> = Structure2D<T> public typealias Matrix<T> = Structure2D<T>
public typealias Vector<T> = Point<T> /**
* Alias or using [Buffer] as a point/vector in a many-dimensional space.
*/
public typealias Point<T> = Buffer<T>
/** /**
* Basic operations on matrices and vectors. Operates on [Matrix]. * Basic operations on matrices and vectors. Operates on [Matrix].
@ -33,13 +36,13 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
/** /**
* Produces a point compatible with matrix space (and possibly optimized for it). * Produces a point compatible with matrix space (and possibly optimized for it).
*/ */
public fun buildVector(size: Int, initializer: A.(Int) -> T): Vector<T> public fun buildVector(size: Int, initializer: A.(Int) -> T): Point<T>
public operator fun Matrix<T>.unaryMinus(): Matrix<T> = buildMatrix(rowNum, colNum) { i, j -> public operator fun Matrix<T>.unaryMinus(): Matrix<T> = buildMatrix(rowNum, colNum) { i, j ->
-get(i, j) -get(i, j)
} }
public operator fun Vector<T>.unaryMinus(): Vector<T> = buildVector(size) { public operator fun Point<T>.unaryMinus(): Point<T> = buildVector(size) {
-get(it) -get(it)
} }
@ -54,7 +57,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
/** /**
* Vector sum * Vector sum
*/ */
public operator fun Vector<T>.plus(other: Vector<T>): Vector<T> = buildVector(size) { public operator fun Point<T>.plus(other: Point<T>): Point<T> = buildVector(size) {
get(it) + other[it] get(it) + other[it]
} }
@ -68,7 +71,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
/** /**
* Vector subtraction * Vector subtraction
*/ */
public operator fun Vector<T>.minus(other: Vector<T>): Vector<T> = buildVector(size) { public operator fun Point<T>.minus(other: Point<T>): Point<T> = buildVector(size) {
get(it) - other[it] get(it) - other[it]
} }
@ -100,7 +103,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
* @param vector the multiplier. * @param vector the multiplier.
* @return the dot product. * @return the dot product.
*/ */
public infix fun Matrix<T>.dot(vector: Vector<T>): Vector<T> { public 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})" } require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
return elementAlgebra { return elementAlgebra {
buildVector(rowNum) { i -> buildVector(rowNum) { i ->
@ -139,7 +142,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
* @param value the multiplier. * @param value the multiplier.
* @receiver the product. * @receiver the product.
*/ */
public operator fun Vector<T>.times(value: T): Vector<T> = public operator fun Point<T>.times(value: T): Point<T> =
buildVector(size) { i -> get(i) * value } buildVector(size) { i -> get(i) * value }
/** /**
@ -149,7 +152,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
* @param v the multiplier. * @param v the multiplier.
* @receiver the product. * @receiver the product.
*/ */
public operator fun T.times(v: Vector<T>): Vector<T> = v * this public operator fun T.times(v: Point<T>): Point<T> = v * this
/** /**
* Gets a feature from the matrix. This function may return some additional features to * Gets a feature from the matrix. This function may return some additional features to

View File

@ -24,7 +24,7 @@ public fun <T : Any, A : Ring<T>> LinearSpace<T, A>.matrix(rows: Int, columns: I
MatrixBuilder(this, rows, columns) MatrixBuilder(this, rows, columns)
@UnstableKMathAPI @UnstableKMathAPI
public fun <T : Any> LinearSpace<T, Ring<T>>.vector(vararg elements: T): Vector<T> { public fun <T : Any> LinearSpace<T, Ring<T>>.vector(vararg elements: T): Point<T> {
return buildVector(elements.size) { elements[it] } return buildVector(elements.size) { elements[it] }
} }

View File

@ -1,86 +0,0 @@
package space.kscience.kmath.linear
//public object RealLinearSpace:
//public object RealLinearSpace : LinearSpace<Double, RealField>, ScaleOperations<Matrix<Double>> {
//
// override val elementAlgebra: RealField get() = RealField
//
// public override fun buildMatrix(
// rows: Int,
// columns: Int,
// initializer: (i: Int, j: Int) -> Double,
// ): Matrix<Double> {
// val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
// return BufferMatrix(rows, columns, buffer)
// }
//
// public fun Matrix<Double>.toBufferMatrix(): BufferMatrix<Double> = if (this is BufferMatrix) this else {
// buildMatrix(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
//
// 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> {
// require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
// val bufferMatrix = toBufferMatrix()
// val otherBufferMatrix = other.toBufferMatrix()
// return buildMatrix(rowNum, other.colNum) { i, j ->
// var res = 0.0
// for (l in 0 until colNum) {
// res += bufferMatrix[i, l] * otherBufferMatrix[l, j]
// }
// res
// }
// }
//
// public override infix fun Matrix<Double>.dot(vector: Point<Double>): Point<Double> {
// require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
// val bufferMatrix = toBufferMatrix()
// return RealBuffer(rowNum) { i ->
// var res = 0.0
// for (j in 0 until colNum) {
// res += bufferMatrix[i, j] * vector[j]
// }
// res
// }
// }
//
// override fun add(a: Matrix<Double>, b: Matrix<Double>): BufferMatrix<Double> {
// 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 buildMatrix(a.rowNum, a.colNum) { i, j ->
// aBufferMatrix[i, j] + bBufferMatrix[i, j]
// }
// }
//
// override fun scale(a: Matrix<Double>, value: Double): BufferMatrix<Double> {
// val bufferMatrix = a.toBufferMatrix()
// 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 multiply(a: Matrix<Double>, k: Number): BufferMatrix<Double> {
//// val aBufferMatrix = a.toBufferMatrix()
//// return produce(a.rowNum, a.colNum) { i, j -> aBufferMatrix[i, j] * k.toDouble() }
//// }
////
//// override fun divide(a: Matrix<Double>, k: Number): BufferMatrix<Double> {
//// val aBufferMatrix = a.toBufferMatrix()
//// return produce(a.rowNum, a.colNum) { i, j -> aBufferMatrix[i, j] / k.toDouble() }
//// }
//}
///**
// * Partially optimized real-valued matrix
// */
//public val LinearSpace.Companion.real: RealLinearSpace get() = RealLinearSpace

View File

@ -27,7 +27,7 @@ public object EjmlLinearSpace : LinearSpace<Double, RealField> {
/** /**
* Converts this vector to EJML one. * Converts this vector to EJML one.
*/ */
public fun Vector<Double>.toEjml(): EjmlVector = when (this) { public fun Point<Double>.toEjml(): EjmlVector = when (this) {
is EjmlVector -> this is EjmlVector -> this
else -> EjmlVector(SimpleMatrix(size, 1).also { else -> EjmlVector(SimpleMatrix(size, 1).also {
(0 until it.numRows()).forEach { row -> it[row, 0] = get(row) } (0 until it.numRows()).forEach { row -> it[row, 0] = get(row) }
@ -41,7 +41,7 @@ public object EjmlLinearSpace : LinearSpace<Double, RealField> {
} }
}) })
override fun buildVector(size: Int, initializer: RealField.(Int) -> Double): Vector<Double> = override fun buildVector(size: Int, initializer: RealField.(Int) -> Double): Point<Double> =
EjmlVector(SimpleMatrix(size, 1).also { EjmlVector(SimpleMatrix(size, 1).also {
(0 until it.numRows()).forEach { row -> it[row, 0] = RealField.initializer(row) } (0 until it.numRows()).forEach { row -> it[row, 0] = RealField.initializer(row) }
}) })
@ -54,7 +54,7 @@ public object EjmlLinearSpace : LinearSpace<Double, RealField> {
public override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlMatrix = public override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlMatrix =
EjmlMatrix(toEjml().origin.mult(other.toEjml().origin)) EjmlMatrix(toEjml().origin.mult(other.toEjml().origin))
public override fun Matrix<Double>.dot(vector: Vector<Double>): EjmlVector = public override fun Matrix<Double>.dot(vector: Point<Double>): EjmlVector =
EjmlVector(toEjml().origin.mult(vector.toEjml().origin)) EjmlVector(toEjml().origin.mult(vector.toEjml().origin))
public override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlMatrix = public override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlMatrix =
@ -63,25 +63,25 @@ public object EjmlLinearSpace : LinearSpace<Double, RealField> {
public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix = public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix =
toEjml().origin.scale(value).wrapMatrix() toEjml().origin.scale(value).wrapMatrix()
override fun Vector<Double>.unaryMinus(): EjmlVector = override fun Point<Double>.unaryMinus(): EjmlVector =
toEjml().origin.negative().wrapVector() toEjml().origin.negative().wrapVector()
override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlMatrix = override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlMatrix =
(toEjml().origin + other.toEjml().origin).wrapMatrix() (toEjml().origin + other.toEjml().origin).wrapMatrix()
override fun Vector<Double>.plus(other: Vector<Double>): EjmlVector = override fun Point<Double>.plus(other: Point<Double>): EjmlVector =
(toEjml().origin + other.toEjml().origin).wrapVector() (toEjml().origin + other.toEjml().origin).wrapVector()
override fun Vector<Double>.minus(other: Vector<Double>): EjmlVector = override fun Point<Double>.minus(other: Point<Double>): EjmlVector =
(toEjml().origin - other.toEjml().origin).wrapVector() (toEjml().origin - other.toEjml().origin).wrapVector()
override fun Double.times(m: Matrix<Double>): EjmlMatrix = override fun Double.times(m: Matrix<Double>): EjmlMatrix =
m.toEjml().origin.scale(this).wrapMatrix() m.toEjml().origin.scale(this).wrapMatrix()
override fun Vector<Double>.times(value: Double): EjmlVector = override fun Point<Double>.times(value: Double): EjmlVector =
toEjml().origin.scale(value).wrapVector() toEjml().origin.scale(value).wrapVector()
override fun Double.times(v: Vector<Double>): EjmlVector = override fun Double.times(v: Point<Double>): EjmlVector =
v.toEjml().origin.scale(this).wrapVector() v.toEjml().origin.scale(this).wrapVector()
} }