forked from kscience/kmath
Merge branch 'dev' into gsl-experiment
# Conflicts: # kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt # kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
This commit is contained in:
commit
0e8f6e29ee
@ -16,18 +16,22 @@
|
||||
- ND4J support module submitting `NDStructure` and `NDAlgebra` over `INDArray`.
|
||||
- Coroutine-deterministic Monte-Carlo scope with a random number generator.
|
||||
- Some minor utilities to `kmath-for-real`.
|
||||
- Generic operation result parameter to `MatrixContext`
|
||||
|
||||
### Changed
|
||||
- Package changed from `scientifik` to `kscience.kmath`.
|
||||
- Gradle version: 6.6 -> 6.7
|
||||
- Gradle version: 6.6 -> 6.7.1
|
||||
- Minor exceptions refactor (throwing `IllegalArgumentException` by argument checks instead of `IllegalStateException`)
|
||||
- `Polynomial` secondary constructor made function.
|
||||
- Kotlin version: 1.3.72 -> 1.4.20-M1
|
||||
- Kotlin version: 1.3.72 -> 1.4.20
|
||||
- `kmath-ast` doesn't depend on heavy `kotlin-reflect` library.
|
||||
- Full autodiff refactoring based on `Symbol`
|
||||
- `kmath-prob` renamed to `kmath-stat`
|
||||
- Grid generators moved to `kmath-for-real`
|
||||
- Use `Point<Double>` instead of specialized type in `kmath-for-real`
|
||||
- Optimized dot product for buffer matrices moved to `kmath-for-real`
|
||||
- EjmlMatrix context is an object
|
||||
- Matrix LUP `inverse` renamed to `inverseWithLUP`
|
||||
|
||||
### Deprecated
|
||||
|
||||
@ -35,6 +39,7 @@
|
||||
- `kmath-koma` module because it doesn't support Kotlin 1.4.
|
||||
- Support of `legacy` JS backend (we will support only IR)
|
||||
- `toGrid` method.
|
||||
- Public visibility of `BufferAccessor2D`
|
||||
|
||||
### Fixed
|
||||
- `symbol` method in `MstExtendedField` (https://github.com/mipt-npm/kmath/pull/140)
|
||||
|
@ -68,7 +68,7 @@ benchmark {
|
||||
// This one matches sourceSet name above
|
||||
|
||||
configurations.register("fast") {
|
||||
warmups = 5 // number of warmup iterations
|
||||
warmups = 1 // number of warmup iterations
|
||||
iterations = 3 // number of iterations
|
||||
iterationTime = 500 // time in seconds per iteration
|
||||
iterationTimeUnit = "ms" // time unity for iterationTime, default is seconds
|
||||
|
@ -1,4 +1,4 @@
|
||||
package kscience.kmath.structures
|
||||
package kscience.kmath.benchmarks
|
||||
|
||||
import org.openjdk.jmh.annotations.Benchmark
|
||||
import org.openjdk.jmh.annotations.Scope
|
@ -1,7 +1,9 @@
|
||||
package kscience.kmath.structures
|
||||
package kscience.kmath.benchmarks
|
||||
|
||||
import kscience.kmath.operations.Complex
|
||||
import kscience.kmath.operations.complex
|
||||
import kscience.kmath.structures.MutableBuffer
|
||||
import kscience.kmath.structures.RealBuffer
|
||||
import org.openjdk.jmh.annotations.Benchmark
|
||||
import org.openjdk.jmh.annotations.Scope
|
||||
import org.openjdk.jmh.annotations.State
|
@ -0,0 +1,50 @@
|
||||
package kscience.kmath.linear
|
||||
|
||||
|
||||
import kotlinx.benchmark.Benchmark
|
||||
import kscience.kmath.commons.linear.CMMatrixContext
|
||||
import kscience.kmath.commons.linear.CMMatrixContext.dot
|
||||
import kscience.kmath.commons.linear.inverse
|
||||
import kscience.kmath.commons.linear.toCM
|
||||
import kscience.kmath.ejml.EjmlMatrixContext
|
||||
import kscience.kmath.ejml.inverse
|
||||
import kscience.kmath.ejml.toEjml
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.structures.Matrix
|
||||
import org.openjdk.jmh.annotations.Scope
|
||||
import org.openjdk.jmh.annotations.State
|
||||
import kotlin.random.Random
|
||||
|
||||
@State(Scope.Benchmark)
|
||||
class LinearAlgebraBenchmark {
|
||||
companion object {
|
||||
val random = Random(1224)
|
||||
val dim = 100
|
||||
|
||||
//creating invertible matrix
|
||||
val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
|
||||
val l = Matrix.real(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 }
|
||||
val matrix = l dot u
|
||||
}
|
||||
|
||||
@Benchmark
|
||||
fun kmathLUPInversion() {
|
||||
MatrixContext.real.inverseWithLUP(matrix)
|
||||
}
|
||||
|
||||
@Benchmark
|
||||
fun cmLUPInversion() {
|
||||
CMMatrixContext {
|
||||
val cm = matrix.toCM() //avoid overhead on conversion
|
||||
inverse(cm)
|
||||
}
|
||||
}
|
||||
|
||||
@Benchmark
|
||||
fun ejmlInverse() {
|
||||
EjmlMatrixContext {
|
||||
val km = matrix.toEjml() //avoid overhead on conversion
|
||||
inverse(km)
|
||||
}
|
||||
}
|
||||
}
|
@ -0,0 +1,60 @@
|
||||
package kscience.kmath.benchmarks
|
||||
|
||||
import kotlinx.benchmark.Benchmark
|
||||
import kscience.kmath.commons.linear.CMMatrixContext
|
||||
import kscience.kmath.commons.linear.CMMatrixContext.dot
|
||||
import kscience.kmath.commons.linear.toCM
|
||||
import kscience.kmath.ejml.EjmlMatrixContext
|
||||
import kscience.kmath.ejml.toEjml
|
||||
import kscience.kmath.linear.real
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.structures.Matrix
|
||||
import org.openjdk.jmh.annotations.Scope
|
||||
import org.openjdk.jmh.annotations.State
|
||||
import kotlin.random.Random
|
||||
|
||||
@State(Scope.Benchmark)
|
||||
class MultiplicationBenchmark {
|
||||
companion object {
|
||||
val random = Random(12224)
|
||||
val dim = 1000
|
||||
|
||||
//creating invertible matrix
|
||||
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 cmMatrix1 = matrix1.toCM()
|
||||
val cmMatrix2 = matrix2.toCM()
|
||||
|
||||
val ejmlMatrix1 = matrix1.toEjml()
|
||||
val ejmlMatrix2 = matrix2.toEjml()
|
||||
}
|
||||
|
||||
@Benchmark
|
||||
fun commonsMathMultiplication() {
|
||||
CMMatrixContext.invoke {
|
||||
cmMatrix1 dot cmMatrix2
|
||||
}
|
||||
}
|
||||
|
||||
@Benchmark
|
||||
fun ejmlMultiplication() {
|
||||
EjmlMatrixContext.invoke {
|
||||
ejmlMatrix1 dot ejmlMatrix2
|
||||
}
|
||||
}
|
||||
|
||||
@Benchmark
|
||||
fun ejmlMultiplicationwithConversion() {
|
||||
val ejmlMatrix1 = matrix1.toEjml()
|
||||
val ejmlMatrix2 = matrix2.toEjml()
|
||||
EjmlMatrixContext.invoke {
|
||||
ejmlMatrix1 dot ejmlMatrix2
|
||||
}
|
||||
}
|
||||
|
||||
@Benchmark
|
||||
fun bufferedMultiplication() {
|
||||
matrix1 dot matrix2
|
||||
}
|
||||
}
|
@ -1,7 +1,8 @@
|
||||
package kscience.kmath.structures
|
||||
package kscience.kmath.benchmarks
|
||||
|
||||
import kscience.kmath.operations.RealField
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.structures.*
|
||||
import org.openjdk.jmh.annotations.Benchmark
|
||||
import org.openjdk.jmh.annotations.Scope
|
||||
import org.openjdk.jmh.annotations.State
|
@ -1,7 +1,10 @@
|
||||
package kscience.kmath.structures
|
||||
package kscience.kmath.benchmarks
|
||||
|
||||
import kscience.kmath.operations.RealField
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.structures.BufferedNDField
|
||||
import kscience.kmath.structures.NDField
|
||||
import kscience.kmath.structures.RealNDField
|
||||
import kscience.kmath.viktor.ViktorNDField
|
||||
import org.jetbrains.bio.viktor.F64Array
|
||||
import org.openjdk.jmh.annotations.Benchmark
|
||||
@ -36,7 +39,7 @@ internal class ViktorBenchmark {
|
||||
|
||||
@Benchmark
|
||||
fun rawViktor() {
|
||||
val one = F64Array.full(init = 1.0, shape = *intArrayOf(dim, dim))
|
||||
val one = F64Array.full(init = 1.0, shape = intArrayOf(dim, dim))
|
||||
var res = one
|
||||
repeat(n) { res = res + one }
|
||||
}
|
@ -1,11 +0,0 @@
|
||||
package kscience.kmath.utils
|
||||
|
||||
import kotlin.contracts.InvocationKind
|
||||
import kotlin.contracts.contract
|
||||
import kotlin.system.measureTimeMillis
|
||||
|
||||
internal inline fun measureAndPrint(title: String, block: () -> Unit) {
|
||||
contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
|
||||
val time = measureTimeMillis(block)
|
||||
println("$title completed in $time millis")
|
||||
}
|
@ -70,6 +70,7 @@ fun main() {
|
||||
//minimize the chi^2 in given starting point. Derivatives are not required, they are already included.
|
||||
val result: OptimizationResult<Double> = chi2.minimize(a to 1.5, b to 0.9, c to 1.0)
|
||||
|
||||
//display a page with plot and numerical results
|
||||
val page = Plotly.page {
|
||||
plot {
|
||||
scatter {
|
||||
|
@ -1,50 +0,0 @@
|
||||
package kscience.kmath.linear
|
||||
|
||||
import kscience.kmath.commons.linear.CMMatrixContext
|
||||
import kscience.kmath.commons.linear.inverse
|
||||
import kscience.kmath.commons.linear.toCM
|
||||
import kscience.kmath.ejml.EjmlMatrixContext
|
||||
import kscience.kmath.ejml.inverse
|
||||
import kscience.kmath.operations.RealField
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.structures.Matrix
|
||||
import kotlin.random.Random
|
||||
import kotlin.system.measureTimeMillis
|
||||
|
||||
fun main() {
|
||||
val random = Random(1224)
|
||||
val dim = 100
|
||||
//creating invertible matrix
|
||||
val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
|
||||
val l = Matrix.real(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 }
|
||||
val matrix = l dot u
|
||||
|
||||
val n = 5000 // iterations
|
||||
|
||||
MatrixContext.real {
|
||||
repeat(50) { inverse(matrix) }
|
||||
val inverseTime = measureTimeMillis { repeat(n) { inverse(matrix) } }
|
||||
println("[kmath] Inversion of $n matrices $dim x $dim finished in $inverseTime millis")
|
||||
}
|
||||
|
||||
//commons-math
|
||||
|
||||
val commonsTime = measureTimeMillis {
|
||||
CMMatrixContext {
|
||||
val cm = matrix.toCM() //avoid overhead on conversion
|
||||
repeat(n) { inverse(cm) }
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis")
|
||||
|
||||
val ejmlTime = measureTimeMillis {
|
||||
(EjmlMatrixContext(RealField)) {
|
||||
val km = matrix.toEjml() //avoid overhead on conversion
|
||||
repeat(n) { inverse(km) }
|
||||
}
|
||||
}
|
||||
|
||||
println("[ejml] Inversion of $n matrices $dim x $dim finished in $ejmlTime millis")
|
||||
}
|
@ -1,38 +0,0 @@
|
||||
package kscience.kmath.linear
|
||||
|
||||
import kscience.kmath.commons.linear.CMMatrixContext
|
||||
import kscience.kmath.commons.linear.toCM
|
||||
import kscience.kmath.ejml.EjmlMatrixContext
|
||||
import kscience.kmath.operations.RealField
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.structures.Matrix
|
||||
import kotlin.random.Random
|
||||
import kotlin.system.measureTimeMillis
|
||||
|
||||
fun main() {
|
||||
val random = Random(12224)
|
||||
val dim = 1000
|
||||
//creating invertible matrix
|
||||
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 }
|
||||
|
||||
// //warmup
|
||||
// matrix1 dot matrix2
|
||||
|
||||
CMMatrixContext {
|
||||
val cmMatrix1 = matrix1.toCM()
|
||||
val cmMatrix2 = matrix2.toCM()
|
||||
val cmTime = measureTimeMillis { cmMatrix1 dot cmMatrix2 }
|
||||
println("CM implementation time: $cmTime")
|
||||
}
|
||||
|
||||
(EjmlMatrixContext(RealField)) {
|
||||
val ejmlMatrix1 = matrix1.toEjml()
|
||||
val ejmlMatrix2 = matrix2.toEjml()
|
||||
val ejmlTime = measureTimeMillis { ejmlMatrix1 dot ejmlMatrix2 }
|
||||
println("EJML implementation time: $ejmlTime")
|
||||
}
|
||||
|
||||
val genericTime = measureTimeMillis { val res = matrix1 dot matrix2 }
|
||||
println("Generic implementation time: $genericTime")
|
||||
}
|
@ -17,3 +17,7 @@ kotlin.sourceSets {
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
readme{
|
||||
maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE
|
||||
}
|
@ -83,11 +83,10 @@ public object ArithmeticsEvaluator : Grammar<MST>() {
|
||||
}
|
||||
|
||||
/**
|
||||
* Tries to parse the string into [MST].
|
||||
* Tries to parse the string into [MST]. Returns [ParseResult] representing expression or error.
|
||||
*
|
||||
* @receiver the string to parse.
|
||||
* @return the [MST] node.
|
||||
* @author Alexander Nozik
|
||||
*/
|
||||
public fun String.tryParseMath(): ParseResult<MST> = ArithmeticsEvaluator.tryParseToEnd(this)
|
||||
|
||||
@ -96,6 +95,5 @@ public fun String.tryParseMath(): ParseResult<MST> = ArithmeticsEvaluator.tryPar
|
||||
*
|
||||
* @receiver the string to parse.
|
||||
* @return the [MST] node.
|
||||
* @author Alexander Nozik
|
||||
*/
|
||||
public fun String.parseMath(): MST = ArithmeticsEvaluator.parseToEnd(this)
|
||||
|
@ -5,8 +5,7 @@ import kscience.kmath.structures.Matrix
|
||||
import kscience.kmath.structures.NDStructure
|
||||
import org.apache.commons.math3.linear.*
|
||||
|
||||
public class CMMatrix(public val origin: RealMatrix, features: Set<MatrixFeature>? = null) :
|
||||
FeaturedMatrix<Double> {
|
||||
public class CMMatrix(public val origin: RealMatrix, features: Set<MatrixFeature>? = null) : FeaturedMatrix<Double> {
|
||||
public override val rowNum: Int get() = origin.rowDimension
|
||||
public override val colNum: Int get() = origin.columnDimension
|
||||
|
||||
@ -55,7 +54,7 @@ public fun Point<Double>.toCM(): CMVector = if (this is CMVector) this else {
|
||||
|
||||
public fun RealVector.toPoint(): CMVector = CMVector(this)
|
||||
|
||||
public object CMMatrixContext : MatrixContext<Double> {
|
||||
public object CMMatrixContext : MatrixContext<Double, CMMatrix> {
|
||||
public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): CMMatrix {
|
||||
val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } }
|
||||
return CMMatrix(Array2DRowRealMatrix(array))
|
||||
@ -79,7 +78,7 @@ public object CMMatrixContext : MatrixContext<Double> {
|
||||
public override fun multiply(a: Matrix<Double>, k: Number): CMMatrix =
|
||||
CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble()))
|
||||
|
||||
public override operator fun Matrix<Double>.times(value: Double): Matrix<Double> =
|
||||
public override operator fun Matrix<Double>.times(value: Double): CMMatrix =
|
||||
produce(rowNum, colNum) { i, j -> get(i, j) * value }
|
||||
}
|
||||
|
||||
|
@ -1,5 +1,3 @@
|
||||
import ru.mipt.npm.gradle.Maturity
|
||||
|
||||
plugins {
|
||||
id("ru.mipt.npm.mpp")
|
||||
id("ru.mipt.npm.native")
|
||||
@ -13,7 +11,7 @@ kotlin.sourceSets.commonMain {
|
||||
|
||||
readme {
|
||||
description = "Core classes, algebra definitions, basic linear algebra"
|
||||
maturity = Maturity.DEVELOPMENT
|
||||
maturity = ru.mipt.npm.gradle.Maturity.DEVELOPMENT
|
||||
propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md"))
|
||||
|
||||
feature(
|
||||
|
@ -9,8 +9,8 @@ import kscience.kmath.structures.*
|
||||
*/
|
||||
public class BufferMatrixContext<T : Any, R : Ring<T>>(
|
||||
public override val elementContext: R,
|
||||
private val bufferFactory: BufferFactory<T>
|
||||
) : GenericMatrixContext<T, R> {
|
||||
private val bufferFactory: BufferFactory<T>,
|
||||
) : GenericMatrixContext<T, R, BufferMatrix<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)
|
||||
@ -22,15 +22,15 @@ public class BufferMatrixContext<T : Any, R : Ring<T>>(
|
||||
}
|
||||
|
||||
@Suppress("OVERRIDE_BY_INLINE")
|
||||
public object RealMatrixContext : GenericMatrixContext<Double, RealField> {
|
||||
public object RealMatrixContext : GenericMatrixContext<Double, RealField, BufferMatrix<Double>> {
|
||||
public override val elementContext: RealField
|
||||
get() = RealField
|
||||
|
||||
public override inline fun produce(
|
||||
rows: Int,
|
||||
columns: Int,
|
||||
initializer: (i: Int, j: Int) -> Double
|
||||
): Matrix<Double> {
|
||||
initializer: (i: Int, j: Int) -> Double,
|
||||
): BufferMatrix<Double> {
|
||||
val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
|
||||
return BufferMatrix(rows, columns, buffer)
|
||||
}
|
||||
@ -43,15 +43,15 @@ public class BufferMatrix<T : Any>(
|
||||
public override val rowNum: Int,
|
||||
public override val colNum: Int,
|
||||
public val buffer: Buffer<out T>,
|
||||
public override val features: Set<MatrixFeature> = emptySet()
|
||||
public override val features: Set<MatrixFeature> = emptySet(),
|
||||
) : FeaturedMatrix<T> {
|
||||
override val shape: IntArray
|
||||
get() = intArrayOf(rowNum, colNum)
|
||||
|
||||
init {
|
||||
require(buffer.size == rowNum * colNum) { "Dimension mismatch for matrix structure" }
|
||||
}
|
||||
|
||||
override val shape: IntArray get() = intArrayOf(rowNum, colNum)
|
||||
|
||||
public override fun suggestFeature(vararg features: MatrixFeature): BufferMatrix<T> =
|
||||
BufferMatrix(rowNum, colNum, buffer, this.features + features)
|
||||
|
||||
@ -86,28 +86,3 @@ public class BufferMatrix<T : Any>(
|
||||
else "Matrix(rowsNum = $rowNum, colNum = $colNum, features=$features)"
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Optimized dot product for real matrices
|
||||
*/
|
||||
public infix fun BufferMatrix<Double>.dot(other: BufferMatrix<Double>): BufferMatrix<Double> {
|
||||
require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
|
||||
val array = DoubleArray(this.rowNum * other.colNum)
|
||||
|
||||
//convert to array to insure there is not memory indirection
|
||||
fun Buffer<out Double>.unsafeArray() = if (this is RealBuffer)
|
||||
array
|
||||
else
|
||||
DoubleArray(size) { get(it) }
|
||||
|
||||
val a = this.buffer.unsafeArray()
|
||||
val b = other.buffer.unsafeArray()
|
||||
|
||||
for (i in (0 until rowNum))
|
||||
for (j in (0 until other.colNum))
|
||||
for (k in (0 until colNum))
|
||||
array[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j]
|
||||
|
||||
val buffer = RealBuffer(array)
|
||||
return BufferMatrix(rowNum, other.colNum, buffer)
|
||||
}
|
||||
|
@ -27,9 +27,8 @@ public interface FeaturedMatrix<T : Any> : Matrix<T> {
|
||||
public inline fun Structure2D.Companion.real(
|
||||
rows: Int,
|
||||
columns: Int,
|
||||
initializer: (Int, Int) -> Double
|
||||
): Matrix<Double> =
|
||||
MatrixContext.real.produce(rows, columns, initializer)
|
||||
initializer: (Int, Int) -> Double,
|
||||
): BufferMatrix<Double> = MatrixContext.real.produce(rows, columns, initializer)
|
||||
|
||||
/**
|
||||
* Build a square matrix from given elements.
|
||||
@ -58,7 +57,7 @@ public inline fun <reified T : Any> Matrix<*>.getFeature(): T? =
|
||||
/**
|
||||
* Diagonal matrix of ones. The matrix is virtual no actual matrix is created
|
||||
*/
|
||||
public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.one(rows: Int, columns: Int): FeaturedMatrix<T> =
|
||||
public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R, *>.one(rows: Int, columns: Int): FeaturedMatrix<T> =
|
||||
VirtualMatrix(rows, columns, DiagonalFeature) { i, j ->
|
||||
if (i == j) elementContext.one else elementContext.zero
|
||||
}
|
||||
@ -67,7 +66,7 @@ public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.one(rows: Int, colu
|
||||
/**
|
||||
* A virtual matrix of zeroes
|
||||
*/
|
||||
public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.zero(rows: Int, columns: Int): FeaturedMatrix<T> =
|
||||
public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R, *>.zero(rows: Int, columns: Int): FeaturedMatrix<T> =
|
||||
VirtualMatrix(rows, columns) { _, _ -> elementContext.zero }
|
||||
|
||||
public class TransposedFeature<T : Any>(public val original: Matrix<T>) : MatrixFeature
|
||||
@ -81,6 +80,4 @@ public fun <T : Any> Matrix<T>.transpose(): Matrix<T> {
|
||||
rowNum,
|
||||
setOf(TransposedFeature(this))
|
||||
) { i, j -> get(j, i) }
|
||||
}
|
||||
|
||||
public infix fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = with(MatrixContext.real) { dot(other) }
|
||||
}
|
@ -1,25 +1,18 @@
|
||||
package kscience.kmath.linear
|
||||
|
||||
import kscience.kmath.operations.Field
|
||||
import kscience.kmath.operations.RealField
|
||||
import kscience.kmath.operations.Ring
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.structures.BufferAccessor2D
|
||||
import kscience.kmath.structures.Matrix
|
||||
import kscience.kmath.structures.Structure2D
|
||||
import kotlin.reflect.KClass
|
||||
import kscience.kmath.operations.*
|
||||
import kscience.kmath.structures.*
|
||||
|
||||
/**
|
||||
* Common implementation of [LUPDecompositionFeature]
|
||||
*/
|
||||
public class LUPDecomposition<T : Any>(
|
||||
public val context: GenericMatrixContext<T, out Field<T>>,
|
||||
public val context: MatrixContext<T, FeaturedMatrix<T>>,
|
||||
public val elementContext: Field<T>,
|
||||
public val lu: Structure2D<T>,
|
||||
public val pivot: IntArray,
|
||||
private val even: Boolean
|
||||
private val even: Boolean,
|
||||
) : LUPDecompositionFeature<T>, DeterminantFeature<T> {
|
||||
public val elementContext: Field<T>
|
||||
get() = context.elementContext
|
||||
|
||||
/**
|
||||
* Returns the matrix L of the decomposition.
|
||||
@ -64,23 +57,25 @@ public class LUPDecomposition<T : Any>(
|
||||
|
||||
}
|
||||
|
||||
public fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.abs(value: T): T =
|
||||
@PublishedApi
|
||||
internal fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, *>.abs(value: T): T =
|
||||
if (value > elementContext.zero) value else elementContext { -value }
|
||||
|
||||
/**
|
||||
* Create a lup decomposition of generic matrix
|
||||
* Create a lup decomposition of generic matrix.
|
||||
*/
|
||||
public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
|
||||
type: KClass<T>,
|
||||
public fun <T : Comparable<T>> MatrixContext<T, FeaturedMatrix<T>>.lup(
|
||||
factory: MutableBufferFactory<T>,
|
||||
elementContext: Field<T>,
|
||||
matrix: Matrix<T>,
|
||||
checkSingular: (T) -> Boolean
|
||||
checkSingular: (T) -> Boolean,
|
||||
): LUPDecomposition<T> {
|
||||
require(matrix.rowNum == matrix.colNum) { "LU decomposition supports only square matrices" }
|
||||
val m = matrix.colNum
|
||||
val pivot = IntArray(matrix.rowNum)
|
||||
|
||||
//TODO just waits for KEEP-176
|
||||
BufferAccessor2D(type, matrix.rowNum, matrix.colNum).run {
|
||||
BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {
|
||||
elementContext {
|
||||
val lu = create(matrix)
|
||||
|
||||
@ -112,14 +107,14 @@ public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.l
|
||||
luRow[col] = sum
|
||||
|
||||
// maintain best permutation choice
|
||||
if (this@lup.abs(sum) > largest) {
|
||||
largest = this@lup.abs(sum)
|
||||
if (abs(sum) > largest) {
|
||||
largest = abs(sum)
|
||||
max = row
|
||||
}
|
||||
}
|
||||
|
||||
// Singularity check
|
||||
check(!checkSingular(this@lup.abs(lu[max, col]))) { "The matrix is singular" }
|
||||
check(!checkSingular(abs(lu[max, col]))) { "The matrix is singular" }
|
||||
|
||||
// Pivot if necessary
|
||||
if (max != col) {
|
||||
@ -143,23 +138,23 @@ public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.l
|
||||
for (row in col + 1 until m) lu[row, col] /= luDiag
|
||||
}
|
||||
|
||||
return LUPDecomposition(this@lup, lu.collect(), pivot, even)
|
||||
return LUPDecomposition(this@lup, elementContext, lu.collect(), pivot, even)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
|
||||
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.lup(
|
||||
matrix: Matrix<T>,
|
||||
checkSingular: (T) -> Boolean
|
||||
): LUPDecomposition<T> = lup(T::class, matrix, checkSingular)
|
||||
noinline checkSingular: (T) -> Boolean,
|
||||
): LUPDecomposition<T> = lup(MutableBuffer.Companion::auto, elementContext, matrix, checkSingular)
|
||||
|
||||
public fun GenericMatrixContext<Double, RealField>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> =
|
||||
lup(Double::class, matrix) { it < 1e-11 }
|
||||
public fun MatrixContext<Double, FeaturedMatrix<Double>>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> =
|
||||
lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 }
|
||||
|
||||
public fun <T : Any> LUPDecomposition<T>.solve(type: KClass<T>, matrix: Matrix<T>): Matrix<T> {
|
||||
public fun <T : Any> LUPDecomposition<T>.solveWithLUP(factory: MutableBufferFactory<T>, matrix: Matrix<T>): FeaturedMatrix<T> {
|
||||
require(matrix.rowNum == pivot.size) { "Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}" }
|
||||
|
||||
BufferAccessor2D(type, matrix.rowNum, matrix.colNum).run {
|
||||
BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {
|
||||
elementContext {
|
||||
// Apply permutations to b
|
||||
val bp = create { _, _ -> zero }
|
||||
@ -201,27 +196,34 @@ public fun <T : Any> LUPDecomposition<T>.solve(type: KClass<T>, matrix: Matrix<T
|
||||
}
|
||||
}
|
||||
|
||||
public inline fun <reified T : Any> LUPDecomposition<T>.solve(matrix: Matrix<T>): Matrix<T> = solve(T::class, matrix)
|
||||
public inline fun <reified T : Any> LUPDecomposition<T>.solveWithLUP(matrix: Matrix<T>): Matrix<T> =
|
||||
solveWithLUP(MutableBuffer.Companion::auto, matrix)
|
||||
|
||||
/**
|
||||
* Solve a linear equation **a*x = b**
|
||||
* Solve a linear equation **a*x = b** using LUP decomposition
|
||||
*/
|
||||
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.solve(
|
||||
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.solveWithLUP(
|
||||
a: Matrix<T>,
|
||||
b: Matrix<T>,
|
||||
checkSingular: (T) -> Boolean
|
||||
): Matrix<T> {
|
||||
noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
|
||||
noinline checkSingular: (T) -> Boolean,
|
||||
): FeaturedMatrix<T> {
|
||||
// Use existing decomposition if it is provided by matrix
|
||||
val decomposition = a.getFeature() ?: lup(T::class, a, checkSingular)
|
||||
return decomposition.solve(T::class, b)
|
||||
val decomposition = a.getFeature() ?: lup(bufferFactory, elementContext, a, checkSingular)
|
||||
return decomposition.solveWithLUP(bufferFactory, b)
|
||||
}
|
||||
|
||||
public fun RealMatrixContext.solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> = solve(a, b) { it < 1e-11 }
|
||||
public fun RealMatrixContext.solveWithLUP(a: Matrix<Double>, b: Matrix<Double>): FeaturedMatrix<Double> =
|
||||
solveWithLUP(a, b) { it < 1e-11 }
|
||||
|
||||
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.inverse(
|
||||
public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.inverseWithLUP(
|
||||
matrix: Matrix<T>,
|
||||
checkSingular: (T) -> Boolean
|
||||
): Matrix<T> = solve(matrix, one(matrix.rowNum, matrix.colNum), checkSingular)
|
||||
noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
|
||||
noinline checkSingular: (T) -> Boolean,
|
||||
): FeaturedMatrix<T> = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular)
|
||||
|
||||
public fun RealMatrixContext.inverse(matrix: Matrix<Double>): Matrix<Double> =
|
||||
solve(matrix, one(matrix.rowNum, matrix.colNum)) { it < 1e-11 }
|
||||
/**
|
||||
* Inverses a square matrix using LUP decomposition. Non square matrix will throw a error.
|
||||
*/
|
||||
public fun RealMatrixContext.inverseWithLUP(matrix: Matrix<Double>): FeaturedMatrix<Double> =
|
||||
solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), Buffer.Companion::real) { it < 1e-11 }
|
||||
|
@ -12,17 +12,17 @@ import kscience.kmath.structures.asSequence
|
||||
/**
|
||||
* Basic operations on matrices. Operates on [Matrix]
|
||||
*/
|
||||
public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
|
||||
public interface MatrixContext<T : Any, out M : Matrix<T>> : SpaceOperations<Matrix<T>> {
|
||||
/**
|
||||
* Produce a matrix with this context and given dimensions
|
||||
*/
|
||||
public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T>
|
||||
public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): M
|
||||
|
||||
public override fun binaryOperation(operation: String, left: Matrix<T>, right: Matrix<T>): Matrix<T> =
|
||||
when (operation) {
|
||||
"dot" -> left dot right
|
||||
else -> super.binaryOperation(operation, left, right)
|
||||
}
|
||||
@Suppress("UNCHECKED_CAST")
|
||||
public override fun binaryOperation(operation: String, left: Matrix<T>, right: Matrix<T>): M = when (operation) {
|
||||
"dot" -> left dot right
|
||||
else -> super.binaryOperation(operation, left, right) as M
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes the dot product of this matrix and another one.
|
||||
@ -31,7 +31,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
|
||||
* @param other the multiplier.
|
||||
* @return the dot product.
|
||||
*/
|
||||
public infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T>
|
||||
public infix fun Matrix<T>.dot(other: Matrix<T>): M
|
||||
|
||||
/**
|
||||
* Computes the dot product of this matrix and a vector.
|
||||
@ -49,7 +49,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
|
||||
* @param value the multiplier.
|
||||
* @receiver the product.
|
||||
*/
|
||||
public operator fun Matrix<T>.times(value: T): Matrix<T>
|
||||
public operator fun Matrix<T>.times(value: T): M
|
||||
|
||||
/**
|
||||
* Multiplies an element by a matrix of it.
|
||||
@ -58,7 +58,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
|
||||
* @param m the multiplier.
|
||||
* @receiver the product.
|
||||
*/
|
||||
public operator fun T.times(m: Matrix<T>): Matrix<T> = m * this
|
||||
public operator fun T.times(m: Matrix<T>): M = m * this
|
||||
|
||||
public companion object {
|
||||
/**
|
||||
@ -71,18 +71,18 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
|
||||
*/
|
||||
public fun <T : Any, R : Ring<T>> buffered(
|
||||
ring: R,
|
||||
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
|
||||
): GenericMatrixContext<T, R> = BufferMatrixContext(ring, bufferFactory)
|
||||
bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
|
||||
): GenericMatrixContext<T, R, BufferMatrix<T>> = BufferMatrixContext(ring, bufferFactory)
|
||||
|
||||
/**
|
||||
* Automatic buffered matrix, unboxed if it is possible
|
||||
*/
|
||||
public inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R> =
|
||||
public inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R, BufferMatrix<T>> =
|
||||
buffered(ring, Buffer.Companion::auto)
|
||||
}
|
||||
}
|
||||
|
||||
public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
|
||||
public interface GenericMatrixContext<T : Any, R : Ring<T>, out M : Matrix<T>> : MatrixContext<T, M> {
|
||||
/**
|
||||
* The ring context for matrix elements
|
||||
*/
|
||||
@ -93,7 +93,7 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
|
||||
*/
|
||||
public fun point(size: Int, initializer: (Int) -> T): Point<T>
|
||||
|
||||
public override infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T> {
|
||||
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})" }
|
||||
|
||||
@ -114,10 +114,10 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
|
||||
}
|
||||
}
|
||||
|
||||
public override operator fun Matrix<T>.unaryMinus(): Matrix<T> =
|
||||
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>): Matrix<T> {
|
||||
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}]"
|
||||
}
|
||||
@ -125,7 +125,7 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
|
||||
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>): Matrix<T> {
|
||||
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}]"
|
||||
}
|
||||
@ -133,11 +133,11 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
|
||||
return produce(rowNum, colNum) { i, j -> elementContext { get(i, j) + b[i, j] } }
|
||||
}
|
||||
|
||||
public override fun multiply(a: Matrix<T>, k: Number): Matrix<T> =
|
||||
public override fun multiply(a: Matrix<T>, k: Number): M =
|
||||
produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] * k } }
|
||||
|
||||
public operator fun Number.times(matrix: FeaturedMatrix<T>): Matrix<T> = matrix * this
|
||||
public operator fun Number.times(matrix: FeaturedMatrix<T>): M = multiply(matrix, this)
|
||||
|
||||
public override operator fun Matrix<T>.times(value: T): Matrix<T> =
|
||||
public override operator fun Matrix<T>.times(value: T): M =
|
||||
produce(rowNum, colNum) { i, j -> elementContext { get(i, j) * value } }
|
||||
}
|
||||
|
@ -38,6 +38,11 @@ public fun <T> Space<T>.average(data: Iterable<T>): T = sum(data) / data.count()
|
||||
*/
|
||||
public fun <T> Space<T>.average(data: Sequence<T>): T = sum(data) / data.count()
|
||||
|
||||
/**
|
||||
* Absolute of the comparable [value]
|
||||
*/
|
||||
public fun <T : Comparable<T>> Space<T>.abs(value: T): T = if (value > zero) value else -value
|
||||
|
||||
/**
|
||||
* Returns the sum of all elements in the iterable in provided space.
|
||||
*
|
||||
|
@ -1,25 +1,31 @@
|
||||
package kscience.kmath.structures
|
||||
|
||||
import kotlin.reflect.KClass
|
||||
|
||||
/**
|
||||
* A context that allows to operate on a [MutableBuffer] as on 2d array
|
||||
*/
|
||||
public class BufferAccessor2D<T : Any>(public val type: KClass<T>, public val rowNum: Int, public val colNum: Int) {
|
||||
internal class BufferAccessor2D<T : Any>(
|
||||
public val rowNum: Int,
|
||||
public val colNum: Int,
|
||||
val factory: MutableBufferFactory<T>,
|
||||
) {
|
||||
public operator fun Buffer<T>.get(i: Int, j: Int): T = get(i + colNum * j)
|
||||
|
||||
public operator fun MutableBuffer<T>.set(i: Int, j: Int, value: T) {
|
||||
set(i + colNum * j, value)
|
||||
}
|
||||
|
||||
public inline fun create(init: (i: Int, j: Int) -> T): MutableBuffer<T> =
|
||||
MutableBuffer.auto(type, rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
|
||||
public inline fun create(crossinline init: (i: Int, j: Int) -> T): MutableBuffer<T> =
|
||||
factory(rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
|
||||
|
||||
public fun create(mat: Structure2D<T>): MutableBuffer<T> = create { i, j -> mat[i, j] }
|
||||
|
||||
//TODO optimize wrapper
|
||||
public fun MutableBuffer<T>.collect(): Structure2D<T> =
|
||||
NDStructure.auto(type, rowNum, colNum) { (i, j) -> get(i, j) }.as2D()
|
||||
public fun MutableBuffer<T>.collect(): Structure2D<T> = NDStructure.build(
|
||||
DefaultStrides(intArrayOf(rowNum, colNum)),
|
||||
factory
|
||||
) { (i, j) ->
|
||||
get(i, j)
|
||||
}.as2D()
|
||||
|
||||
public inner class Row(public val buffer: MutableBuffer<T>, public val rowIndex: Int) : MutableBuffer<T> {
|
||||
override val size: Int get() = colNum
|
||||
@ -30,7 +36,7 @@ public class BufferAccessor2D<T : Any>(public val type: KClass<T>, public val ro
|
||||
buffer[rowIndex, index] = value
|
||||
}
|
||||
|
||||
override fun copy(): MutableBuffer<T> = MutableBuffer.auto(type, colNum) { get(it) }
|
||||
override fun copy(): MutableBuffer<T> = factory(colNum) { get(it) }
|
||||
override operator fun iterator(): Iterator<T> = (0 until colNum).map(::get).iterator()
|
||||
|
||||
}
|
||||
|
@ -1,5 +1,6 @@
|
||||
package kscience.kmath.linear
|
||||
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.structures.Matrix
|
||||
import kscience.kmath.structures.NDStructure
|
||||
import kscience.kmath.structures.as2D
|
||||
@ -38,7 +39,7 @@ class MatrixTest {
|
||||
infix fun Matrix<Double>.pow(power: Int): Matrix<Double> {
|
||||
var res = this
|
||||
repeat(power - 1) {
|
||||
res = res dot this
|
||||
res = RealMatrixContext.invoke { res dot this@pow }
|
||||
}
|
||||
return res
|
||||
}
|
||||
|
@ -9,7 +9,7 @@ class RealLUSolverTest {
|
||||
@Test
|
||||
fun testInvertOne() {
|
||||
val matrix = MatrixContext.real.one(2, 2)
|
||||
val inverted = MatrixContext.real.inverse(matrix)
|
||||
val inverted = MatrixContext.real.inverseWithLUP(matrix)
|
||||
assertEquals(matrix, inverted)
|
||||
}
|
||||
|
||||
@ -37,7 +37,7 @@ class RealLUSolverTest {
|
||||
1.0, 3.0
|
||||
)
|
||||
|
||||
val inverted = MatrixContext.real.inverse(matrix)
|
||||
val inverted = MatrixContext.real.inverseWithLUP(matrix)
|
||||
|
||||
val expected = Matrix.square(
|
||||
0.375, -0.125,
|
||||
|
@ -18,3 +18,7 @@ kotlin.sourceSets {
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
readme{
|
||||
maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE
|
||||
}
|
||||
|
@ -42,7 +42,7 @@ public interface DMatrix<T, R : Dimension, C : Dimension> : Structure2D<T> {
|
||||
* An inline wrapper for a Matrix
|
||||
*/
|
||||
public inline class DMatrixWrapper<T, R : Dimension, C : Dimension>(
|
||||
public val structure: Structure2D<T>
|
||||
private val structure: Structure2D<T>
|
||||
) : DMatrix<T, R, C> {
|
||||
override val shape: IntArray get() = structure.shape
|
||||
override operator fun get(i: Int, j: Int): T = structure[i, j]
|
||||
@ -81,7 +81,7 @@ public inline class DPointWrapper<T, D : Dimension>(public val point: Point<T>)
|
||||
/**
|
||||
* Basic operations on dimension-safe matrices. Operates on [Matrix]
|
||||
*/
|
||||
public inline class DMatrixContext<T : Any, Ri : Ring<T>>(public val context: GenericMatrixContext<T, Ri>) {
|
||||
public inline class DMatrixContext<T : Any, Ri : Ring<T>>(public val context: GenericMatrixContext<T, Ri, Matrix<T>>) {
|
||||
public inline fun <reified R : Dimension, reified C : Dimension> Matrix<T>.coerce(): DMatrix<T, R, C> {
|
||||
require(rowNum == Dimension.dim<R>().toInt()) {
|
||||
"Row number mismatch: expected ${Dimension.dim<R>()} but found $rowNum"
|
||||
|
@ -2,22 +2,21 @@ package kscience.kmath.ejml
|
||||
|
||||
import kscience.kmath.linear.MatrixContext
|
||||
import kscience.kmath.linear.Point
|
||||
import kscience.kmath.operations.Space
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.structures.Matrix
|
||||
import org.ejml.simple.SimpleMatrix
|
||||
|
||||
/**
|
||||
* Converts this matrix to EJML one.
|
||||
*/
|
||||
public fun Matrix<Double>.toEjml(): EjmlMatrix =
|
||||
if (this is EjmlMatrix) this else EjmlMatrixContext.produce(rowNum, colNum) { i, j -> get(i, j) }
|
||||
|
||||
/**
|
||||
* Represents context of basic operations operating with [EjmlMatrix].
|
||||
*
|
||||
* @author Iaroslav Postovalov
|
||||
*/
|
||||
public class EjmlMatrixContext(private val space: Space<Double>) : MatrixContext<Double> {
|
||||
/**
|
||||
* Converts this matrix to EJML one.
|
||||
*/
|
||||
public fun Matrix<Double>.toEjml(): EjmlMatrix =
|
||||
if (this is EjmlMatrix) this else produce(rowNum, colNum) { i, j -> get(i, j) }
|
||||
public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> {
|
||||
|
||||
/**
|
||||
* Converts this vector to EJML one.
|
||||
@ -47,12 +46,10 @@ public class EjmlMatrixContext(private val space: Space<Double>) : MatrixContext
|
||||
EjmlMatrix(toEjml().origin - b.toEjml().origin)
|
||||
|
||||
public override fun multiply(a: Matrix<Double>, k: Number): EjmlMatrix =
|
||||
produce(a.rowNum, a.colNum) { i, j -> space { a[i, j] * k } }
|
||||
produce(a.rowNum, a.colNum) { i, j -> a[i, j] * k.toDouble() }
|
||||
|
||||
public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix =
|
||||
EjmlMatrix(toEjml().origin.scale(value))
|
||||
|
||||
public companion object
|
||||
}
|
||||
|
||||
/**
|
||||
|
@ -1,13 +1,14 @@
|
||||
package kscience.kmath.real
|
||||
|
||||
import kscience.kmath.linear.FeaturedMatrix
|
||||
import kscience.kmath.linear.MatrixContext
|
||||
import kscience.kmath.linear.RealMatrixContext.elementContext
|
||||
import kscience.kmath.linear.VirtualMatrix
|
||||
import kscience.kmath.linear.inverseWithLUP
|
||||
import kscience.kmath.misc.UnstableKMathAPI
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.operations.sum
|
||||
import kscience.kmath.structures.Buffer
|
||||
import kscience.kmath.structures.Matrix
|
||||
import kscience.kmath.structures.RealBuffer
|
||||
import kscience.kmath.structures.asIterable
|
||||
import kotlin.math.pow
|
||||
@ -24,7 +25,7 @@ import kotlin.math.pow
|
||||
* Functions that help create a real (Double) matrix
|
||||
*/
|
||||
|
||||
public typealias RealMatrix = Matrix<Double>
|
||||
public typealias RealMatrix = FeaturedMatrix<Double>
|
||||
|
||||
public fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double): RealMatrix =
|
||||
MatrixContext.real.produce(rowNum, colNum, initializer)
|
||||
@ -86,18 +87,6 @@ public operator fun Double.minus(matrix: RealMatrix): RealMatrix =
|
||||
// row, col -> matrix[row, col] / this
|
||||
//}
|
||||
|
||||
/*
|
||||
* Per-element (!) square and power operations
|
||||
*/
|
||||
|
||||
public fun RealMatrix.square(): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { row, col ->
|
||||
this[row, col].pow(2)
|
||||
}
|
||||
|
||||
public fun RealMatrix.pow(n: Int): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { i, j ->
|
||||
this[i, j].pow(n)
|
||||
}
|
||||
|
||||
/*
|
||||
* Operations on two matrices (per-element!)
|
||||
*/
|
||||
@ -157,3 +146,35 @@ public fun RealMatrix.sum(): Double = elements().map { (_, value) -> value }.sum
|
||||
public fun RealMatrix.min(): Double? = elements().map { (_, value) -> value }.minOrNull()
|
||||
public fun RealMatrix.max(): Double? = elements().map { (_, value) -> value }.maxOrNull()
|
||||
public fun RealMatrix.average(): Double = elements().map { (_, value) -> value }.average()
|
||||
|
||||
public inline fun RealMatrix.map(transform: (Double) -> Double): RealMatrix =
|
||||
MatrixContext.real.produce(rowNum, colNum) { i, j ->
|
||||
transform(get(i, j))
|
||||
}
|
||||
|
||||
/**
|
||||
* Inverse a square real matrix using LUP decomposition
|
||||
*/
|
||||
public fun RealMatrix.inverseWithLUP(): RealMatrix = MatrixContext.real.inverseWithLUP(this)
|
||||
|
||||
//extended operations
|
||||
|
||||
public fun RealMatrix.pow(p: Double): RealMatrix = map { it.pow(p) }
|
||||
|
||||
public fun RealMatrix.pow(p: Int): RealMatrix = map { it.pow(p) }
|
||||
|
||||
public fun exp(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.exp(it) }
|
||||
|
||||
public fun sqrt(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.sqrt(it) }
|
||||
|
||||
public fun RealMatrix.square(): RealMatrix = map { it.pow(2) }
|
||||
|
||||
public fun sin(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.sin(it) }
|
||||
|
||||
public fun cos(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.cos(it) }
|
||||
|
||||
public fun tan(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.tan(it) }
|
||||
|
||||
public fun ln(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.ln(it) }
|
||||
|
||||
public fun log10(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.log10(it) }
|
@ -3,7 +3,6 @@ package kscience.kmath.real
|
||||
import kscience.kmath.linear.Point
|
||||
import kscience.kmath.operations.Norm
|
||||
import kscience.kmath.structures.Buffer
|
||||
import kscience.kmath.structures.RealBuffer
|
||||
import kscience.kmath.structures.asBuffer
|
||||
import kscience.kmath.structures.asIterable
|
||||
import kotlin.math.pow
|
||||
@ -11,17 +10,11 @@ import kotlin.math.sqrt
|
||||
|
||||
public typealias RealVector = Point<Double>
|
||||
|
||||
public inline fun RealVector(size: Int, init: (Int) -> Double): RealVector = RealBuffer(size, init)
|
||||
public fun RealVector(vararg doubles: Double): RealVector = RealBuffer(doubles)
|
||||
|
||||
public fun DoubleArray.asVector(): RealVector = asBuffer()
|
||||
public fun List<Double>.asVector(): RealVector = asBuffer()
|
||||
|
||||
|
||||
public object VectorL2Norm : Norm<Point<out Number>, Double> {
|
||||
override fun norm(arg: Point<out Number>): Double = sqrt(arg.asIterable().sumByDouble(Number::toDouble))
|
||||
}
|
||||
|
||||
public operator fun Buffer.Companion.invoke(vararg doubles: Double): RealVector = doubles.asBuffer()
|
||||
|
||||
/**
|
||||
* Fill the vector of given [size] with given [value]
|
||||
@ -36,12 +29,6 @@ public inline fun RealVector.map(transform: (Double) -> Double): RealVector =
|
||||
public inline fun RealVector.mapIndexed(transform: (index: Int, value: Double) -> Double): RealVector =
|
||||
Buffer.real(size) { transform(it, get(it)) }
|
||||
|
||||
public fun RealVector.pow(p: Double): RealVector = map { it.pow(p) }
|
||||
|
||||
public fun RealVector.pow(p: Int): RealVector = map { it.pow(p) }
|
||||
|
||||
public fun exp(vector: RealVector): RealVector = vector.map { kotlin.math.exp(it) }
|
||||
|
||||
public operator fun RealVector.plus(other: RealVector): RealVector =
|
||||
mapIndexed { index, value -> value + other[index] }
|
||||
|
||||
@ -71,3 +58,25 @@ public operator fun RealVector.div(other: RealVector): RealVector =
|
||||
public operator fun RealVector.div(number: Number): RealVector = map { it / number.toDouble() }
|
||||
|
||||
public operator fun Number.div(vector: RealVector): RealVector = vector.map { toDouble() / it }
|
||||
|
||||
//extended operations
|
||||
|
||||
public fun RealVector.pow(p: Double): RealVector = map { it.pow(p) }
|
||||
|
||||
public fun RealVector.pow(p: Int): RealVector = map { it.pow(p) }
|
||||
|
||||
public fun exp(vector: RealVector): RealVector = vector.map { kotlin.math.exp(it) }
|
||||
|
||||
public fun sqrt(vector: RealVector): RealVector = vector.map { kotlin.math.sqrt(it) }
|
||||
|
||||
public fun RealVector.square(): RealVector = map { it.pow(2) }
|
||||
|
||||
public fun sin(vector: RealVector): RealVector = vector.map { kotlin.math.sin(it) }
|
||||
|
||||
public fun cos(vector: RealVector): RealVector = vector.map { kotlin.math.cos(it) }
|
||||
|
||||
public fun tan(vector: RealVector): RealVector = vector.map { kotlin.math.tan(it) }
|
||||
|
||||
public fun ln(vector: RealVector): RealVector = vector.map { kotlin.math.ln(it) }
|
||||
|
||||
public fun log10(vector: RealVector): RealVector = vector.map { kotlin.math.log10(it) }
|
||||
|
@ -0,0 +1,31 @@
|
||||
package kscience.kmath.real
|
||||
|
||||
import kscience.kmath.linear.BufferMatrix
|
||||
import kscience.kmath.structures.Buffer
|
||||
import kscience.kmath.structures.RealBuffer
|
||||
|
||||
|
||||
/**
|
||||
* Optimized dot product for real matrices
|
||||
*/
|
||||
public infix fun BufferMatrix<Double>.dot(other: BufferMatrix<Double>): BufferMatrix<Double> {
|
||||
require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
|
||||
val resultArray = DoubleArray(this.rowNum * other.colNum)
|
||||
|
||||
//convert to array to insure there is no memory indirection
|
||||
fun Buffer<out Double>.unsafeArray() = if (this is RealBuffer)
|
||||
this.array
|
||||
else
|
||||
DoubleArray(size) { get(it) }
|
||||
|
||||
val a = this.buffer.unsafeArray()
|
||||
val b = other.buffer.unsafeArray()
|
||||
|
||||
for (i in (0 until rowNum))
|
||||
for (j in (0 until other.colNum))
|
||||
for (k in (0 until colNum))
|
||||
resultArray[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j]
|
||||
|
||||
val buffer = RealBuffer(resultArray)
|
||||
return BufferMatrix(rowNum, other.colNum, buffer)
|
||||
}
|
@ -1,34 +1,33 @@
|
||||
package kaceince.kmath.real
|
||||
|
||||
import kscience.kmath.linear.MatrixContext
|
||||
import kscience.kmath.linear.asMatrix
|
||||
import kscience.kmath.linear.transpose
|
||||
import kscience.kmath.linear.*
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.real.RealVector
|
||||
import kscience.kmath.real.plus
|
||||
import kscience.kmath.structures.Buffer
|
||||
import kotlin.test.Test
|
||||
import kotlin.test.assertEquals
|
||||
|
||||
internal class RealVectorTest {
|
||||
@Test
|
||||
fun testSum() {
|
||||
val vector1 = RealVector(5) { it.toDouble() }
|
||||
val vector2 = RealVector(5) { 5 - it.toDouble() }
|
||||
val vector1 = Buffer.real(5) { it.toDouble() }
|
||||
val vector2 = Buffer.real(5) { 5 - it.toDouble() }
|
||||
val sum = vector1 + vector2
|
||||
assertEquals(5.0, sum[2])
|
||||
}
|
||||
|
||||
@Test
|
||||
fun testVectorToMatrix() {
|
||||
val vector = RealVector(5) { it.toDouble() }
|
||||
val vector = Buffer.real(5) { it.toDouble() }
|
||||
val matrix = vector.asMatrix()
|
||||
assertEquals(4.0, matrix[4, 0])
|
||||
}
|
||||
|
||||
@Test
|
||||
fun testDot() {
|
||||
val vector1 = RealVector(5) { it.toDouble() }
|
||||
val vector2 = RealVector(5) { 5 - it.toDouble() }
|
||||
val vector1 = Buffer.real(5) { it.toDouble() }
|
||||
val vector2 = Buffer.real(5) { 5 - it.toDouble() }
|
||||
val matrix1 = vector1.asMatrix()
|
||||
val matrix2 = vector2.asMatrix().transpose()
|
||||
val product = MatrixContext.real { matrix1 dot matrix2 }
|
||||
|
@ -3,7 +3,6 @@ package kscience.kmath.histogram
|
||||
import kscience.kmath.linear.Point
|
||||
import kscience.kmath.operations.SpaceOperations
|
||||
import kscience.kmath.operations.invoke
|
||||
import kscience.kmath.real.asVector
|
||||
import kscience.kmath.structures.*
|
||||
import kotlin.math.floor
|
||||
|
||||
@ -123,8 +122,8 @@ public class RealHistogram(
|
||||
*```
|
||||
*/
|
||||
public fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): RealHistogram = RealHistogram(
|
||||
ranges.map(ClosedFloatingPointRange<Double>::start).asVector(),
|
||||
ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asVector()
|
||||
ranges.map(ClosedFloatingPointRange<Double>::start).asBuffer(),
|
||||
ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asBuffer()
|
||||
)
|
||||
|
||||
/**
|
||||
|
@ -4,6 +4,8 @@ import kscience.kmath.histogram.RealHistogram
|
||||
import kscience.kmath.histogram.fill
|
||||
import kscience.kmath.histogram.put
|
||||
import kscience.kmath.real.RealVector
|
||||
import kscience.kmath.real.invoke
|
||||
import kscience.kmath.structures.Buffer
|
||||
import kotlin.random.Random
|
||||
import kotlin.test.*
|
||||
|
||||
|
@ -1,8 +1,8 @@
|
||||
package kscience.kmath.histogram
|
||||
|
||||
import kscience.kmath.real.RealVector
|
||||
import kscience.kmath.real.asVector
|
||||
import kscience.kmath.structures.Buffer
|
||||
import kscience.kmath.structures.asBuffer
|
||||
import java.util.*
|
||||
import kotlin.math.floor
|
||||
|
||||
@ -16,7 +16,7 @@ public class UnivariateBin(
|
||||
//TODO add weighting
|
||||
public override val value: Number get() = counter.sum()
|
||||
|
||||
public override val center: RealVector get() = doubleArrayOf(position).asVector()
|
||||
public override val center: RealVector get() = doubleArrayOf(position).asBuffer()
|
||||
public override val dimension: Int get() = 1
|
||||
|
||||
public operator fun contains(value: Double): Boolean = value in (position - size / 2)..(position + size / 2)
|
||||
|
Loading…
Reference in New Issue
Block a user