Merge branch 'dev' into gsl-experiment

# Conflicts:
#	kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt

kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt

@@ -7,10 +7,16 @@ import kscience.kmath.structures.asBuffer
 import kotlin.math.sqrt
 
 /**
- * A 2d structure plus optional matrix-specific features
+ * A [Matrix] that holds [MatrixFeature] objects.
+ *
+ * @param T the type of items.
  */
 public interface FeaturedMatrix<T : Any> : Matrix<T> {
-    override val shape: IntArray get() = intArrayOf(rowNum, colNum)
+    public override val shape: IntArray get() = intArrayOf(rowNum, colNum)
+
+    /**
+     * The set of features this matrix possesses.
+     */
     public val features: Set<MatrixFeature>
 
     /**

kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt

@@ -4,16 +4,15 @@ import kscience.kmath.operations.*
 import kscience.kmath.structures.*
 
 /**
- * Common implementation of [LUPDecompositionFeature]
+ * Common implementation of [LupDecompositionFeature].
  */
 public class LUPDecomposition<T : Any>(
     public val context: MatrixContext<T, FeaturedMatrix<T>>,
     public val elementContext: Field<T>,
-    public val lu: Structure2D<T>,
+    public val lu: Matrix<T>,
     public val pivot: IntArray,
     private val even: Boolean,
-) : LUPDecompositionFeature<T>, DeterminantFeature<T> {
-
+) : LupDecompositionFeature<T>, DeterminantFeature<T> {
     /**
      * Returns the matrix L of the decomposition.
      *
@@ -151,7 +150,10 @@ public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext
 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>.solveWithLUP(factory: MutableBufferFactory<T>, matrix: Matrix<T>): FeaturedMatrix<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(matrix.rowNum, matrix.colNum, factory).run {
@@ -217,4 +219,4 @@ public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext
     matrix: Matrix<T>,
     noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
     noinline checkSingular: (T) -> Boolean,
-): FeaturedMatrix<T> = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular)
+): FeaturedMatrix<T> = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular)

kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixFeatures.kt

@@ -1,62 +1,154 @@
 package kscience.kmath.linear
 
 /**
- * A marker interface representing some matrix feature like diagonal, sparse, zero, etc. Features used to optimize matrix
- * operations performance in some cases.
+ * A marker interface representing some properties of matrices or additional transformations of them. Features are used
+ * to optimize matrix operations performance in some cases or retrieve the APIs.
  */
 public interface MatrixFeature
 
 /**
- * The matrix with this feature is considered to have only diagonal non-null elements
+ * Matrices with this feature are considered to have only diagonal non-null elements.
  */
 public object DiagonalFeature : MatrixFeature
 
 /**
- * Matrix with this feature has all zero elements
+ * Matrices with this feature have all zero elements.
  */
 public object ZeroFeature : MatrixFeature
 
 /**
- * Matrix with this feature have unit elements on diagonal and zero elements in all other places
+ * Matrices with this feature have unit elements on diagonal and zero elements in all other places.
  */
 public object UnitFeature : MatrixFeature
 
 /**
- * Inverted matrix feature
+ * Matrices with this feature can be inverted: [inverse] = `a`<sup>-1</sup> where `a` is the owning matrix.
+ *
+ * @param T the type of matrices' items.
  */
 public interface InverseMatrixFeature<T : Any> : MatrixFeature {
+    /**
+     * The inverse matrix of the matrix that owns this feature.
+     */
     public val inverse: FeaturedMatrix<T>
 }
 
 /**
- * A determinant container
+ * Matrices with this feature can compute their determinant.
  */
 public interface DeterminantFeature<T : Any> : MatrixFeature {
+    /**
+     * The determinant of the matrix that owns this feature.
+     */
     public val determinant: T
 }
 
+/**
+ * Produces a [DeterminantFeature] where the [DeterminantFeature.determinant] is [determinant].
+ *
+ * @param determinant the value of determinant.
+ * @return a new [DeterminantFeature].
+ */
 @Suppress("FunctionName")
 public fun <T : Any> DeterminantFeature(determinant: T): DeterminantFeature<T> = object : DeterminantFeature<T> {
     override val determinant: T = determinant
 }
 
 /**
- * Lower triangular matrix
+ * Matrices with this feature are lower triangular ones.
  */
 public object LFeature : MatrixFeature
 
 /**
- * Upper triangular feature
+ * Matrices with this feature are upper triangular ones.
  */
 public object UFeature : MatrixFeature
 
 /**
- * TODO add documentation
+ * Matrices with this feature support LU factorization with partial pivoting: *[p] &middot; a = [l] &middot; [u]* where
+ * *a* is the owning matrix.
+ *
+ * @param T the type of matrices' items.
  */
-public interface LUPDecompositionFeature<T : Any> : MatrixFeature {
+public interface LupDecompositionFeature<T : Any> : MatrixFeature {
+    /**
+     * The lower triangular matrix in this decomposition. It may have [LFeature].
+     */
     public val l: FeaturedMatrix<T>
+
+    /**
+     * The upper triangular matrix in this decomposition. It may have [UFeature].
+     */
     public val u: FeaturedMatrix<T>
+
+    /**
+     * The permutation matrix in this decomposition.
+     */
     public val p: FeaturedMatrix<T>
 }
 
+/**
+ * Matrices with this feature are orthogonal ones: *a &middot; a<sup>T</sup> = u* where *a* is the owning matrix, *u*
+ * is the unit matrix ([UnitFeature]).
+ */
+public object OrthogonalFeature : MatrixFeature
+
+/**
+ * Matrices with this feature support QR factorization: *a = [q] &middot; [r]* where *a* is the owning matrix.
+ *
+ * @param T the type of matrices' items.
+ */
+public interface QRDecompositionFeature<T : Any> : MatrixFeature {
+    /**
+     * The orthogonal matrix in this decomposition. It may have [OrthogonalFeature].
+     */
+    public val q: FeaturedMatrix<T>
+
+    /**
+     * The upper triangular matrix in this decomposition. It may have [UFeature].
+     */
+    public val r: FeaturedMatrix<T>
+}
+
+/**
+ * Matrices with this feature support Cholesky factorization: *a = [l] &middot; [l]<sup>H</sup>* where *a* is the
+ * owning matrix.
+ *
+ * @param T the type of matrices' items.
+ */
+public interface CholeskyDecompositionFeature<T : Any> : MatrixFeature {
+    /**
+     * The triangular matrix in this decomposition. It may have either [UFeature] or [LFeature].
+     */
+    public val l: FeaturedMatrix<T>
+}
+
+/**
+ * Matrices with this feature support SVD: *a = [u] &middot; [s] &middot; [v]<sup>H</sup>* where *a* is the owning
+ * matrix.
+ *
+ * @param T the type of matrices' items.
+ */
+public interface SingularValueDecompositionFeature<T : Any> : MatrixFeature {
+    /**
+     * The matrix in this decomposition. It is unitary, and it consists from left singular vectors.
+     */
+    public val u: FeaturedMatrix<T>
+
+    /**
+     * The matrix in this decomposition. Its main diagonal elements are singular values.
+     */
+    public val s: FeaturedMatrix<T>
+
+    /**
+     * The matrix in this decomposition. It is unitary, and it consists from right singular vectors.
+     */
+    public val v: FeaturedMatrix<T>
+
+    /**
+     * The buffer of singular values of this SVD.
+     */
+    public val singularValues: Point<T>
+}
+
 //TODO add sparse matrix feature

kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure2D.kt

@@ -1,12 +1,40 @@
 package kscience.kmath.structures
 
 /**
- * A structure that is guaranteed to be two-dimensional
+ * A structure that is guaranteed to be two-dimensional.
+ *
+ * @param T the type of items.
  */
 public interface Structure2D<T> : NDStructure<T> {
+    /**
+     * The number of rows in this structure.
+     */
     public val rowNum: Int get() = shape[0]
+
+    /**
+     * The number of columns in this structure.
+     */
     public val colNum: Int get() = shape[1]
 
+    /**
+     * The buffer of rows of this structure. It gets elements from the structure dynamically.
+     */
+    public val rows: Buffer<Buffer<T>>
+        get() = VirtualBuffer(rowNum) { i -> VirtualBuffer(colNum) { j -> get(i, j) } }
+
+    /**
+     * The buffer of columns of this structure. It gets elements from the structure dynamically.
+     */
+    public val columns: Buffer<Buffer<T>>
+        get() = VirtualBuffer(colNum) { j -> VirtualBuffer(rowNum) { i -> get(i, j) } }
+
+    /**
+     * Retrieves an element from the structure by two indices.
+     *
+     * @param i the first index.
+     * @param j the second index.
+     * @return an element.
+     */
     public operator fun get(i: Int, j: Int): T
 
     override operator fun get(index: IntArray): T {
@@ -14,15 +42,9 @@ public interface Structure2D<T> : NDStructure<T> {
         return get(index[0], index[1])
     }
 
-    public val rows: Buffer<Buffer<T>>
-        get() = VirtualBuffer(rowNum) { i -> VirtualBuffer(colNum) { j -> get(i, j) } }
-
-    public val columns: Buffer<Buffer<T>>
-        get() = VirtualBuffer(colNum) { j -> VirtualBuffer(rowNum) { i -> get(i, j) } }
-
     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))
+        for (i in 0 until rowNum)
+            for (j in 0 until colNum) yield(intArrayOf(i, j) to get(i, j))
     }
 
     public companion object
@@ -47,4 +69,9 @@ public fun <T> NDStructure<T>.as2D(): Structure2D<T> = if (shape.size == 2)
 else
     error("Can't create 2d-structure from ${shape.size}d-structure")
 
+/**
+ * Alias for [Structure2D] with more familiar name.
+ *
+ * @param T the type of items.
+ */
 public typealias Matrix<T> = Structure2D<T>

kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt

@@ -1,10 +1,8 @@
 package kscience.kmath.ejml
 
-import kscience.kmath.linear.DeterminantFeature
-import kscience.kmath.linear.FeaturedMatrix
-import kscience.kmath.linear.LUPDecompositionFeature
-import kscience.kmath.linear.MatrixFeature
+import kscience.kmath.linear.*
 import kscience.kmath.structures.NDStructure
+import kscience.kmath.structures.RealBuffer
 import org.ejml.dense.row.factory.DecompositionFactory_DDRM
 import org.ejml.simple.SimpleMatrix
 
@@ -14,7 +12,7 @@ import org.ejml.simple.SimpleMatrix
  * @property origin the underlying [SimpleMatrix].
  * @author Iaroslav Postovalov
  */
-public class EjmlMatrix(public val origin: SimpleMatrix, features: Set<MatrixFeature>? = null) :
+public class EjmlMatrix(public val origin: SimpleMatrix, features: Set<MatrixFeature> = emptySet()) :
     FeaturedMatrix<Double> {
     public override val rowNum: Int
         get() = origin.numRows()
@@ -22,32 +20,63 @@ public class EjmlMatrix(public val origin: SimpleMatrix, features: Set<MatrixFea
     public override val colNum: Int
         get() = origin.numCols()
 
-    public override val features: Set<MatrixFeature> = setOf(
-        object : LUPDecompositionFeature<Double>, DeterminantFeature<Double> {
-            override val determinant: Double
-                get() = origin.determinant()
+    public override val shape: IntArray by lazy { intArrayOf(rowNum, colNum) }
 
-            private val lup by lazy {
-                val ludecompositionF64 = DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols())
-                    .also { it.decompose(origin.ddrm.copy()) }
+    public override val features: Set<MatrixFeature> = hashSetOf(
+        object : InverseMatrixFeature<Double> {
+            override val inverse: FeaturedMatrix<Double> by lazy { EjmlMatrix(origin.invert()) }
+        },
 
-                Triple(
-                    EjmlMatrix(SimpleMatrix(ludecompositionF64.getRowPivot(null))),
-                    EjmlMatrix(SimpleMatrix(ludecompositionF64.getLower(null))),
-                    EjmlMatrix(SimpleMatrix(ludecompositionF64.getUpper(null))),
-                )
+        object : DeterminantFeature<Double> {
+            override val determinant: Double by lazy(origin::determinant)
+        },
+
+        object : SingularValueDecompositionFeature<Double> {
+            private val svd by lazy {
+                DecompositionFactory_DDRM.svd(origin.numRows(), origin.numCols(), true, true, false)
+                    .apply { decompose(origin.ddrm.copy()) }
             }
 
-            override val l: FeaturedMatrix<Double>
-                get() = lup.second
+            override val u: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getU(null, false))) }
+            override val s: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getW(null))) }
+            override val v: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(svd.getV(null, false))) }
+            override val singularValues: Point<Double> by lazy { RealBuffer(svd.singularValues) }
+        },
 
-            override val u: FeaturedMatrix<Double>
-                get() = lup.third
+        object : QRDecompositionFeature<Double> {
+            private val qr by lazy {
+                DecompositionFactory_DDRM.qr().apply { decompose(origin.ddrm.copy()) }
+            }
 
-            override val p: FeaturedMatrix<Double>
-                get() = lup.first
-        }
-    ) union features.orEmpty()
+            override val q: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(qr.getQ(null, false))) }
+            override val r: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(qr.getR(null, false))) }
+        },
+
+        object : CholeskyDecompositionFeature<Double> {
+            override val l: FeaturedMatrix<Double> by lazy {
+                val cholesky =
+                    DecompositionFactory_DDRM.chol(rowNum, true).apply { decompose(origin.ddrm.copy()) }
+
+                EjmlMatrix(SimpleMatrix(cholesky.getT(null)), setOf(LFeature))
+            }
+        },
+
+        object : LupDecompositionFeature<Double> {
+            private val lup by lazy {
+                DecompositionFactory_DDRM.lu(origin.numRows(), origin.numCols()).apply { decompose(origin.ddrm.copy()) }
+            }
+
+            override val l: FeaturedMatrix<Double> by lazy {
+                EjmlMatrix(SimpleMatrix(lup.getLower(null)), setOf(LFeature))
+            }
+
+            override val u: FeaturedMatrix<Double> by lazy {
+                EjmlMatrix(SimpleMatrix(lup.getUpper(null)), setOf(UFeature))
+            }
+
+            override val p: FeaturedMatrix<Double> by lazy { EjmlMatrix(SimpleMatrix(lup.getRowPivot(null))) }
+        },
+    ) union features
 
     public override fun suggestFeature(vararg features: MatrixFeature): EjmlMatrix =
         EjmlMatrix(origin, this.features + features)

kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt

@@ -17,7 +17,6 @@ public fun Matrix<Double>.toEjml(): EjmlMatrix =
  * @author Iaroslav Postovalov
  */
 public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> {
-
     /**
      * Converts this vector to EJML one.
      */
@@ -33,6 +32,11 @@ public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> {
             }
         })
 
+    override fun point(size: Int, initializer: (Int) -> Double): Point<Double> =
+        EjmlVector(SimpleMatrix(size, 1).also {
+            (0 until it.numRows()).forEach { row -> it[row, 0] = initializer(row) }
+        })
+
     public override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlMatrix =
         EjmlMatrix(toEjml().origin.mult(other.toEjml().origin))
 
@@ -73,12 +77,3 @@ public fun EjmlMatrixContext.solve(a: Matrix<Double>, b: Matrix<Double>): EjmlMa
  */
 public fun EjmlMatrixContext.solve(a: Matrix<Double>, b: Point<Double>): EjmlVector =
     EjmlVector(a.toEjml().origin.solve(b.toEjml().origin))
-
-/**
- * Returns the inverse of given matrix: b = a^(-1).
- *
- * @param a the matrix.
- * @return the inverse of this matrix.
- * @author Iaroslav Postovalov
- */
-public fun EjmlMatrixContext.inverse(a: Matrix<Double>): EjmlMatrix = EjmlMatrix(a.toEjml().origin.invert())

kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt

@@ -1,7 +1,7 @@
 package kscience.kmath.ejml
 
 import kscience.kmath.linear.DeterminantFeature
-import kscience.kmath.linear.LUPDecompositionFeature
+import kscience.kmath.linear.LupDecompositionFeature
 import kscience.kmath.linear.MatrixFeature
 import kscience.kmath.linear.getFeature
 import org.ejml.dense.row.factory.DecompositionFactory_DDRM
@@ -44,7 +44,7 @@ internal class EjmlMatrixTest {
         val w = EjmlMatrix(m)
         val det = w.getFeature<DeterminantFeature<Double>>() ?: fail()
         assertEquals(m.determinant(), det.determinant)
-        val lup = w.getFeature<LUPDecompositionFeature<Double>>() ?: fail()
+        val lup = w.getFeature<LupDecompositionFeature<Double>>() ?: fail()
 
         val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows(), m.numCols())
             .also { it.decompose(m.ddrm.copy()) }