From a9c084773b3931c37ce5ec3db5d7d58abcaffc6d Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Tue, 9 Oct 2018 12:22:35 +0300 Subject: [PATCH] Linear algebra in progress --- build.gradle | 2 +- .../scientifik/kmath/commons/CommonsMatrix.kt | 3 + .../kmath/linear/LUDecomposition.kt | 321 ++++++++++++++++++ .../scientifik/kmath/linear/LinearAlgrebra.kt | 250 ++++++++++++++ .../scientifik/kmath/operations/Complex.kt | 41 +++ .../scientifik/kmath/operations/Fields.kt | 41 --- .../kmath/structures/BufferNDField.kt | 21 +- .../kmath/structures/LinearAlgrebra.kt | 169 --------- .../scientifik/kmath/structures/NDArray.kt | 20 +- .../scientifik/kmath/structures/NDArrays.kt | 56 +-- .../kmath/structures/ArrayMatrixTest.kt | 1 + .../kmath/structures/SimpleNDFieldTest.kt | 4 +- .../scientifik/kmath/structures/_NDArrays.kt | 2 +- .../scientifik/kmath/structures/_NDArrays.kt | 16 +- settings.gradle | 1 + 15 files changed, 708 insertions(+), 240 deletions(-) create mode 100644 kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt delete mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt diff --git a/build.gradle b/build.gradle index 8db73094a..04700f377 100644 --- a/build.gradle +++ b/build.gradle @@ -1,5 +1,5 @@ buildscript { - ext.kotlin_version = '1.3.0-rc-116' + ext.kotlin_version = '1.3.0-rc-146' repositories { jcenter() diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt new file mode 100644 index 000000000..a4a9a67ae --- /dev/null +++ b/kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt @@ -0,0 +1,3 @@ +package scientifik.kmath.commons + +//val solver: DecompositionSolver \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt new file mode 100644 index 000000000..6d329809f --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt @@ -0,0 +1,321 @@ +package scientifik.kmath.linear + +import scientifik.kmath.operations.Field +import scientifik.kmath.structures.MutableNDArray +import scientifik.kmath.structures.NDArray +import scientifik.kmath.structures.NDArrays + +/** + * Calculates the LUP-decomposition of a square matrix. + * + * The LUP-decomposition of a matrix A consists of three matrices L, U and + * P that satisfy: PA = LU. L is lower triangular (with unit + * diagonal terms), U is upper triangular and P is a permutation matrix. All + * matrices are mm. + * + * As shown by the presence of the P matrix, this decomposition is + * implemented using partial pivoting. + * + * This class is based on the class with similar name from the + * [JAMA](http://math.nist.gov/javanumerics/jama/) library. + * + * * a [getP][.getP] method has been added, + * * the `det` method has been renamed as [ getDeterminant][.getDeterminant], + * * the `getDoublePivot` method has been removed (but the int based + * [getPivot][.getPivot] method has been kept), + * * the `solve` and `isNonSingular` methods have been replaced + * by a [getSolver][.getSolver] method and the equivalent methods + * provided by the returned [DecompositionSolver]. + * + * + * @see [MathWorld](http://mathworld.wolfram.com/LUDecomposition.html) + * + * @see [Wikipedia](http://en.wikipedia.org/wiki/LU_decomposition) + * + * @since 2.0 (changed to concrete class in 3.0) + * + * @param matrix The matrix to decompose. + * @param singularityThreshold threshold (based on partial row norm) + * under which a matrix is considered singular + * @throws NonSquareMatrixException if matrix is not square + + */ +abstract class LUDecomposition>(val matrix: Matrix) { + + private val field get() = matrix.context.field + /** Entries of LU decomposition. */ + internal val lu: NDArray + /** 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 { + val pair = matrix.context.field.calculateLU() + lu = pair.first + pivot = pair.second + } + + /** + * Returns the matrix L of the decomposition. + * + * L is a lower-triangular matrix + * @return the L matrix (or null if decomposed matrix is singular) + */ + val l: Matrix by lazy { + matrix.context.produce { i, j -> + when { + j < i -> lu[i, j] + j == i -> matrix.context.field.one + else -> matrix.context.field.zero + } + } + } + + + /** + * Returns the matrix U of the decomposition. + * + * U is an upper-triangular matrix + * @return the U matrix (or null if decomposed matrix is singular) + */ + val u: Matrix by lazy { + matrix.context.produce { i, j -> + if (j >= i) lu[i, j] else field.zero + } + } + + /** + * Returns the P rows permutation matrix. + * + * P is a sparse matrix with exactly one element set to 1.0 in + * each row and each column, all other elements being set to 0.0. + * + * The positions of the 1 elements are given by the [ pivot permutation vector][.getPivot]. + * @return the P rows permutation matrix (or null if decomposed matrix is singular) + * @see .getPivot + */ + val p: Matrix by lazy { + matrix.context.produce { i, j -> + //TODO ineffective. Need sparse matrix for that + if (j == pivot[i]) field.one else field.zero + } + } + + /** + * Return the determinant of the matrix + * @return determinant of the matrix + */ + val determinant: T + get() { + with(matrix.context.field) { + var determinant = if (even) one else -one + for (i in 0 until matrix.rows) { + determinant *= lu[i, i] + } + return determinant + } + } + +// /** +// * Get a solver for finding the A X = B solution in exact linear +// * sense. +// * @return a solver +// */ +// val solver: DecompositionSolver +// get() = Solver(lu, pivot, singular) + + /** + * In-place transformation for [MutableNDArray], using given transformation for each element + */ + operator fun MutableNDArray.set(i: Int, j: Int, value: T) { + this[listOf(i, j)] = value + } + + abstract fun isSingular(value: T): Boolean + + private fun Field.calculateLU(): Pair, IntArray> { + if (matrix.rows != matrix.columns) { + error("LU decomposition supports only square matrices") + } + + fun T.abs() = if (this > zero) this else -this + + val m = matrix.columns + val pivot = IntArray(matrix.rows) + //TODO fix performance + val lu: MutableNDArray = NDArrays.createMutable(matrix.context.field, listOf(matrix.rows, matrix.columns)) { index -> matrix[index[0], index[1]] } + + // Initialize permutation array and parity + for (row in 0 until m) { + pivot[row] = row + } + even = true + + // Loop over columns + for (col in 0 until m) { + + // upper + for (row in 0 until col) { + var sum = lu[row, col] + for (i in 0 until row) { + sum -= lu[row, i] * lu[i, col] + } + lu[row, col] = sum + } + + // lower + val max = (col until m).maxBy { row -> + var sum = lu[row, col] + for (i in 0 until col) { + sum -= lu[row, i] * lu[i, col] + } + //luRow[col] = sum + lu[row, col] = sum + + sum.abs() + } ?: col + + // Singularity check + if (isSingular(lu[max, col].abs())) { + error("Singular matrix") + } + + // Pivot if necessary + if (max != col) { + //var tmp = zero + //val luMax = lu[max] + //val luCol = lu[col] + for (i in 0 until m) { + lu[max, i] = lu[col, i] + lu[col, i] = lu[max, i] + } + val temp = pivot[max] + pivot[max] = pivot[col] + pivot[col] = temp + even = !even + } + + // Divide the lower elements by the "winning" diagonal elt. + val luDiag = lu[col, col] + for (row in col + 1 until m) { + lu[row, col] /= luDiag + } + } + return Pair(lu, pivot) + } + + /** + * Returns the pivot permutation vector. + * @return the pivot permutation vector + * @see .getP + */ + fun getPivot(): IntArray { + return pivot.copyOf() + } + + companion object { + /** Default bound to determine effective singularity in LU decomposition. */ + private const val DEFAULT_TOO_SMALL = 1e-11 + } +} + +class RealLUDecomposition(matrix: Matrix, private val singularityThreshold: Double = 1e-11) : LUDecomposition(matrix) { + override fun isSingular(value: Double): Boolean { + return value < singularityThreshold + } +} + + +/** Specialized solver. */ +class RealLUSolver : LinearSolver { + +// +// /** {@inheritDoc} */ +// override fun solve(b: RealVector): RealVector { +// val m = pivot.size +// if (b.getDimension() != m) { +// throw DimensionMismatchException(b.getDimension(), m) +// } +// if (singular) { +// throw SingularMatrixException() +// } +// +// val bp = DoubleArray(m) +// +// // Apply permutations to b +// for (row in 0 until m) { +// bp[row] = b.getEntry(pivot[row]) +// } +// +// // Solve LY = b +// for (col in 0 until m) { +// val bpCol = bp[col] +// for (i in col + 1 until m) { +// bp[i] -= bpCol * lu[i][col] +// } +// } +// +// // Solve UX = Y +// for (col in m - 1 downTo 0) { +// bp[col] /= lu[col][col] +// val bpCol = bp[col] +// for (i in 0 until col) { +// bp[i] -= bpCol * lu[i][col] +// } +// } +// +// return ArrayRealVector(bp, false) +// } + + + fun decompose(mat: Matrix, threshold: Double = 1e-11): RealLUDecomposition = RealLUDecomposition(mat, threshold) + + override fun solve(a: Matrix, b: Matrix): Matrix { + val decomposition = decompose(a, a.context.field.zero) + + if (b.rows != a.rows) { + error("Matrix dimension mismatch expected ${a.rows}, but got ${b.rows}") + } + + // Apply permutations to b + val bp = Array(a.rows) { DoubleArray(b.columns) } + for (row in 0 until a.rows) { + val bpRow = bp[row] + val pRow = decomposition.pivot[row] + for (col in 0 until b.columns) { + bpRow[col] = b[pRow, col] + } + } + + // Solve LY = b + for (col in 0 until a.rows) { + val bpCol = bp[col] + for (i in col + 1 until a.rows) { + val bpI = bp[i] + val luICol = decomposition.lu[i, col] + for (j in 0 until b.columns) { + bpI[j] -= bpCol[j] * luICol + } + } + } + + // Solve UX = Y + for (col in a.rows - 1 downTo 0) { + val bpCol = bp[col] + val luDiag = decomposition.lu[col, col] + for (j in 0 until b.columns) { + bpCol[j] /= luDiag + } + for (i in 0 until col) { + val bpI = bp[i] + val luICol = decomposition.lu[i, col] + for (j in 0 until b.columns) { + bpI[j] -= bpCol[j] * luICol + } + } + } + + return a.context.produce { i, j -> bp[i][j] } + } +} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt new file mode 100644 index 000000000..07295e0aa --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt @@ -0,0 +1,250 @@ +package scientifik.kmath.linear + +import scientifik.kmath.operations.DoubleField +import scientifik.kmath.operations.Field +import scientifik.kmath.operations.Space +import scientifik.kmath.operations.SpaceElement +import scientifik.kmath.structures.NDArray +import scientifik.kmath.structures.NDArrays.createFactory +import scientifik.kmath.structures.NDFieldFactory +import scientifik.kmath.structures.realNDFieldFactory + +/** + * The space for linear elements. Supports scalar product alongside with standard linear operations. + * @param T type of individual element of the vector or matrix + * @param V the type of vector space element + */ +abstract class LinearSpace>(val rows: Int, val columns: Int, val field: Field) : Space { + + /** + * Produce the element of this space + */ + abstract fun produce(initializer: (Int, Int) -> T): V + + /** + * Produce new linear space with given dimensions. The space produced could be raised from cache since [LinearSpace] does not have mutable elements + */ + abstract fun produceSpace(rows: Int, columns: Int): LinearSpace + + override val zero: V by lazy { + produce { _, _ -> field.zero } + } + + val one: V by lazy { + produce { i, j -> if (i == j) field.one else field.zero } + } + + override fun add(a: V, b: V): V { + return produce { i, j -> with(field) { a[i, j] + b[i, j] } } + } + + override fun multiply(a: V, k: Double): V { + //TODO it is possible to implement scalable linear elements which normed values and adjustable scale to save memory and processing poser + return produce { i, j -> with(field) { a[i, j] * k } } + } + + /** + * Dot product. Throws exception on dimension mismatch + */ + fun multiply(a: V, b: V): V { + if (a.rows != b.columns) { + //TODO replace by specific exception + error("Dimension mismatch in linear structure dot product: [${a.rows},${a.columns}]*[${b.rows},${b.columns}]") + } + return produceSpace(a.rows, b.columns).produce { i, j -> + (0..a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) } + } + } + + infix fun V.dot(b: V): V = multiply(this, b) +} + +/** + * A specialized [LinearSpace] which works with vectors + */ +abstract class VectorSpace>(size: Int, field: Field) : LinearSpace(size, 1, field) + +/** + * A matrix-like structure + */ +interface Matrix { + val context: LinearSpace> + /** + * Number of rows + */ + val rows: Int + /** + * Number of columns + */ + val columns: Int + + /** + * Get element in row [i] and column [j]. Throws error in case of call ounside structure dimensions + */ + operator fun get(i: Int, j: Int): T + + fun transpose(): Matrix { + return object : Matrix { + override val context: LinearSpace> = this@Matrix.context + override val rows: Int = this@Matrix.columns + override val columns: Int = this@Matrix.rows + override fun get(i: Int, j: Int): T = this@Matrix[j, i] + } + } +} + +interface Vector : Matrix { + override val context: VectorSpace> + override val columns: Int + get() = 1 + + operator fun get(i: Int) = get(i, 0) +} + + +/** + * NDArray-based implementation of vector space. By default uses slow [SimpleNDField], but could be overridden with custom [NDField] factory. + */ +class ArraySpace( + rows: Int, + columns: Int, + field: Field, + val ndFactory: NDFieldFactory = createFactory(field) +) : LinearSpace>(rows, columns, field) { + + val ndField by lazy { + ndFactory(listOf(rows, columns)) + } + + override fun produce(initializer: (Int, Int) -> T): Matrix = ArrayMatrix(this, initializer) + + override fun produceSpace(rows: Int, columns: Int): ArraySpace { + return ArraySpace(rows, columns, field, ndFactory) + } +} + +class ArrayVectorSpace( + size: Int, + field: Field, + val ndFactory: NDFieldFactory = createFactory(field) +) : VectorSpace>(size, field) { + val ndField by lazy { + ndFactory(listOf(size)) + } + + override fun produce(initializer: (Int, Int) -> T): Vector = produceVector { i -> initializer(i, 0) } + + fun produceVector(initializer: (Int) -> T): Vector = ArrayVector(this, initializer) + + override fun produceSpace(rows: Int, columns: Int): LinearSpace> { + TODO("not implemented") //To change body of created functions use File | Settings | File Templates. + } +} + +/** + * Member of [ArraySpace] which wraps 2-D array + */ +class ArrayMatrix internal constructor(override val context: ArraySpace, val array: NDArray) : Matrix, SpaceElement, ArraySpace> { + + constructor(context: ArraySpace, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) }) + + override val rows: Int get() = context.rows + + override val columns: Int get() = context.columns + + override fun get(i: Int, j: Int): T { + return array[i, j] + } + + override val self: ArrayMatrix get() = this +} + + +class ArrayVector internal constructor(override val context: ArrayVectorSpace, val array: NDArray) : Vector, SpaceElement, ArrayVectorSpace> { + + constructor(context: ArrayVectorSpace, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) }) + + init { + if (context.columns != 1) { + error("Vector must have single column") + } + if (context.rows != array.shape[0]) { + error("Array dimension mismatch") + } + } + + //private val array = context.ndField.produce { list -> initializer(list[0]) } + + + override val rows: Int get() = context.rows + + override val columns: Int = 1 + + override fun get(i: Int, j: Int): T { + return array[i] + } + + override val self: ArrayVector get() = this + +} + +/** + * A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors + */ +interface LinearSolver { + fun solve(a: Matrix, b: Matrix): Matrix + fun solve(a: Matrix, b: Vector): Vector = solve(a, b as Matrix).toVector() + fun inverse(a: Matrix): Matrix = solve(a, a.context.one) +} + +/** + * Create vector with custom field + */ +fun vector(size: Int, field: Field, initializer: (Int) -> T) = + ArrayVector(ArrayVectorSpace(size, field), initializer) + +/** + * Create vector of [Double] + */ +fun realVector(size: Int, initializer: (Int) -> Double) = + ArrayVector(ArrayVectorSpace(size, DoubleField, realNDFieldFactory), initializer) + +/** + * Convert vector to array (copying content of array) + */ +fun Array.asVector(field: Field) = vector(size, field) { this[it] } + +fun DoubleArray.asVector() = realVector(this.size) { this[it] } + +/** + * Create [ArrayMatrix] with custom field + */ +fun matrix(rows: Int, columns: Int, field: Field, initializer: (Int, Int) -> T) = + ArrayMatrix(ArraySpace(rows, columns, field), initializer) + +/** + * Create [ArrayMatrix] of doubles. + */ +fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) = + ArrayMatrix(ArraySpace(rows, columns, DoubleField, realNDFieldFactory), initializer) + + +/** + * Convert matrix to vector if it is possible + */ +fun Matrix.toVector(): Vector { + return when { + this is Vector -> return this + this.columns == 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) } + } + } + else -> error("Can't convert matrix with more than one column to vector") + } +} + diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt new file mode 100644 index 000000000..dd3be256e --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt @@ -0,0 +1,41 @@ +package scientifik.kmath.operations + +/** + * A field for complex numbers + */ +object ComplexField : Field { + override val zero: Complex = Complex(0.0, 0.0) + + override fun add(a: Complex, b: Complex): Complex = Complex(a.re + b.re, a.im + b.im) + + override fun multiply(a: Complex, k: Double): Complex = Complex(a.re * k, a.im * k) + + override val one: Complex = Complex(1.0, 0.0) + + override fun multiply(a: Complex, b: Complex): Complex = Complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re) + + override fun divide(a: Complex, b: Complex): Complex = Complex(a.re * b.re + a.im * b.im, a.re * b.im - a.im * b.re) / b.square + +} + +/** + * Complex number class + */ +data class Complex(val re: Double, val im: Double) : FieldElement { + override val self: Complex get() = this + override val context: ComplexField + get() = ComplexField + + /** + * A complex conjugate + */ + val conjugate: Complex + get() = Complex(re, -im) + + val square: Double + get() = re * re + im * im + + val abs: Double + get() = kotlin.math.sqrt(square) + +} \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt index 0d2600d82..eeba76222 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt @@ -1,7 +1,6 @@ package scientifik.kmath.operations import kotlin.math.pow -import kotlin.math.sqrt /** * Field for real values @@ -51,46 +50,6 @@ data class Real(val value: Double) : Number(), FieldElement { } -/** - * A field for complex numbers - */ -object ComplexField : Field { - override val zero: Complex = Complex(0.0, 0.0) - - override fun add(a: Complex, b: Complex): Complex = Complex(a.re + b.re, a.im + b.im) - - override fun multiply(a: Complex, k: Double): Complex = Complex(a.re * k, a.im * k) - - override val one: Complex = Complex(1.0, 0.0) - - override fun multiply(a: Complex, b: Complex): Complex = Complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re) - - override fun divide(a: Complex, b: Complex): Complex = Complex(a.re * b.re + a.im * b.im, a.re * b.im - a.im * b.re) / b.square - -} - -/** - * Complex number class - */ -data class Complex(val re: Double, val im: Double) : FieldElement { - override val self: Complex get() = this - override val context: ComplexField - get() = ComplexField - - /** - * A complex conjugate - */ - val conjugate: Complex - get() = Complex(re, -im) - - val square: Double - get() = re * re + im * im - - val module: Double - get() = sqrt(square) - -} - /** * A field for double without boxing. Does not produce appropriate field element */ diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt index 9fe7efe87..aab0f1b97 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt @@ -9,6 +9,11 @@ import scientifik.kmath.operations.Field interface Buffer { operator fun get(index: Int): T operator fun set(index: Int, value: T) + + /** + * A shallow copy of the buffer + */ + fun copy(): Buffer } /** @@ -63,8 +68,16 @@ abstract class BufferNDField(shape: List, field: Field) : NDField( return BufferNDArray(this, buffer) } + /** + * Produce mutable NDArray instance + */ + fun produceMutable(initializer: (List) -> T): MutableNDArray { + val buffer = createBuffer(capacity) { initializer(index(it)) } + return MutableBufferedNDArray(this, buffer) + } - class BufferNDArray(override val context: BufferNDField, val data: Buffer) : NDArray { + + private open class BufferNDArray(override val context: BufferNDField, val data: Buffer) : NDArray { override fun get(vararg index: Int): T { return data[context.offset(index.asList())] @@ -88,6 +101,12 @@ abstract class BufferNDField(shape: List, field: Field) : NDField( override val self: NDArray get() = this } + + private class MutableBufferedNDArray(context: BufferNDField, data: Buffer): BufferNDArray(context,data), MutableNDArray{ + override operator fun set(index: List, value: T){ + data[context.offset(index)] = value + } + } } diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt deleted file mode 100644 index a5e6f2ce1..000000000 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt +++ /dev/null @@ -1,169 +0,0 @@ -package scientifik.kmath.structures - -import scientifik.kmath.operations.DoubleField -import scientifik.kmath.operations.Field -import scientifik.kmath.operations.Space -import scientifik.kmath.operations.SpaceElement -import scientifik.kmath.structures.NDArrays.createSimpleNDFieldFactory - -/** - * The space for linear elements. Supports scalar product alongside with standard linear operations. - * @param T type of individual element of the vector or matrix - * @param V the type of vector space element - */ -abstract class LinearSpace>(val rows: Int, val columns: Int, val field: Field) : Space { - - /** - * Produce the element of this space - */ - abstract fun produce(initializer: (Int, Int) -> T): V - - /** - * Produce new linear space with given dimensions. The space produced could be raised from cache since [LinearSpace] does not have mutable elements - */ - abstract fun produceSpace(rows: Int, columns: Int): LinearSpace - - override val zero: V by lazy { - produce { _, _ -> field.zero } - } - - override fun add(a: V, b: V): V { - return produce { i, j -> with(field) { a[i, j] + b[i, j] } } - } - - override fun multiply(a: V, k: Double): V { - //TODO it is possible to implement scalable linear elements which normed values and adjustable scale to save memory and processing poser - return produce { i, j -> with(field) { a[i, j] * k } } - } - - /** - * Dot product. Throws exception on dimension mismatch - */ - fun multiply(a: V, b: V): V { - if (a.rows != b.columns) { - //TODO replace by specific exception - error("Dimension mismatch in linear structure dot product: [${a.rows},${a.columns}]*[${b.rows},${b.columns}]") - } - return produceSpace(a.rows, b.columns).produce { i, j -> - (0..a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) } - } - } - - infix fun V.dot(b: V): V = multiply(this, b) -} - -/** - * A matrix-like structure that is not dependent on specific space implementation - */ -interface LinearStructure { - /** - * Number of rows - */ - val rows: Int - /** - * Number of columns - */ - val columns: Int - - /** - * Get element in row [i] and column [j]. Throws error in case of call ounside structure dimensions - */ - operator fun get(i: Int, j: Int): T - - fun transpose(): LinearStructure { - return object : LinearStructure { - override val rows: Int = this@LinearStructure.columns - override val columns: Int = this@LinearStructure.rows - override fun get(i: Int, j: Int): T = this@LinearStructure[j, i] - } - } -} - -interface Vector : LinearStructure { - override val columns: Int - get() = 1 - - operator fun get(i: Int) = get(i, 0) -} - - -/** - * NDArray-based implementation of vector space. By default uses slow [SimpleNDField], but could be overridden with custom [NDField] factory. - */ -class ArraySpace( - rows: Int, - columns: Int, - field: Field, - val ndFactory: NDFieldFactory = createSimpleNDFieldFactory(field) -) : LinearSpace>(rows, columns, field) { - - val ndField by lazy { - ndFactory(listOf(rows, columns)) - } - - override fun produce(initializer: (Int, Int) -> T): LinearStructure = ArrayMatrix(this, initializer) - - override fun produceSpace(rows: Int, columns: Int): LinearSpace> { - return ArraySpace(rows, columns, field, ndFactory) - } -} - -/** - * Member of [ArraySpace] which wraps 2-D array - */ -class ArrayMatrix(override val context: ArraySpace, initializer: (Int, Int) -> T) : LinearStructure, SpaceElement, ArraySpace> { - - private val array = context.ndField.produce { list -> initializer(list[0], list[1]) } - - //val list: List> = (0 until rows).map { i -> (0 until columns).map { j -> initializer(i, j) } } - - override val rows: Int get() = context.rows - - override val columns: Int get() = context.columns - - override fun get(i: Int, j: Int): T { - return array[i, j] - } - - override val self: ArrayMatrix get() = this -} - - -class ArrayVector(override val context: ArraySpace, initializer: (Int) -> T) : Vector, SpaceElement, ArraySpace> { - - init { - if (context.columns != 1) { - error("Vector must have single column") - } - } - - private val array = context.ndField.produce { list -> initializer(list[0]) } - - - override val rows: Int get() = context.rows - - override val columns: Int = 1 - - override fun get(i: Int, j: Int): T { - return array[i] - } - - override val self: ArrayVector get() = this - -} - -fun vector(size: Int, field: Field, initializer: (Int) -> T) = - ArrayVector(ArraySpace(size, 1, field), initializer) - -fun realVector(size: Int, initializer: (Int) -> Double) = - ArrayVector(ArraySpace(size, 1, DoubleField, realNDFieldFactory), initializer) - -fun Array.asVector(field: Field) = vector(size, field) { this[it] } - -fun DoubleArray.asVector() = realVector(this.size) { this[it] } - -fun matrix(rows: Int, columns: Int, field: Field, initializer: (Int, Int) -> T) = - ArrayMatrix(ArraySpace(rows, columns, field), initializer) - -fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) = - ArrayMatrix(ArraySpace(rows, columns, DoubleField, realNDFieldFactory), initializer) diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArray.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArray.kt index b1ceed5b1..85f60704f 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArray.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArray.kt @@ -73,7 +73,9 @@ abstract class NDField(val shape: List, val field: Field) : Field : FieldElement, NDField> { /** @@ -125,6 +127,22 @@ interface NDArray : FieldElement, NDField> { } } +/** + * In-place mutable [NDArray] + */ +interface MutableNDArray : NDArray { + operator fun set(index: List, value: T) +} + +/** + * In-place transformation for [MutableNDArray], using given transformation for each element + */ +fun MutableNDArray.transformInPlace(action: (List, T) -> T) { + for ((index, oldValue) in this) { + this[index] = action(index, oldValue) + } +} + /** * Element by element application of any operation on elements to the whole array. Just like in numpy */ diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt index ee5ff0a74..739e91832 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt @@ -9,6 +9,28 @@ typealias NDFieldFactory = (shape: List) -> NDField */ expect val realNDFieldFactory: NDFieldFactory + +class SimpleNDField(field: Field, shape: List) : BufferNDField(shape, field) { + override fun createBuffer(capacity: Int, initializer: (Int) -> T): Buffer { + val array = ArrayList(capacity) + (0 until capacity).forEach { + array.add(initializer(it)) + } + + return BufferOfObjects(array) + } + + private class BufferOfObjects(val array: ArrayList) : Buffer { + override fun get(index: Int): T = array[index] + + override fun set(index: Int, value: T) { + array[index] = value + } + + override fun copy(): Buffer = BufferOfObjects(ArrayList(array)) + } +} + object NDArrays { /** * Create a platform-optimized NDArray of doubles @@ -29,26 +51,22 @@ object NDArrays { return realNDArray(listOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) } } - class SimpleNDField(field: Field, shape: List) : BufferNDField(shape, field) { - override fun createBuffer(capacity: Int, initializer: (Int) -> T): Buffer { - val array = ArrayList(capacity) - (0 until capacity).forEach { - array.add(initializer(it)) - } + /** + * Simple boxing NDField + */ + fun createFactory(field: Field): NDFieldFactory = { shape -> SimpleNDField(field, shape) } - return object : Buffer { - override fun get(index: Int): T = array[index] - - override fun set(index: Int, value: T) { - array[index] = initializer(index) - } - } - } - } - - fun createSimpleNDFieldFactory(field: Field): NDFieldFactory = { list -> SimpleNDField(field, list) } - - fun simpleNDArray(field: Field, shape: List, initializer: (List) -> T): NDArray { + /** + * Simple boxing NDArray + */ + fun create(field: Field, shape: List, initializer: (List) -> T): NDArray { return SimpleNDField(field, shape).produce { initializer(it) } } + + /** + * Mutable boxing NDArray + */ + fun createMutable(field: Field, shape: List, initializer: (List) -> T): MutableNDArray { + return SimpleNDField(field, shape).produceMutable { initializer(it) } + } } \ No newline at end of file diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt index 10d383a22..77a436626 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt @@ -1,5 +1,6 @@ package scientifik.kmath.structures +import scientifik.kmath.linear.realVector import kotlin.test.Test import kotlin.test.assertEquals diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt index 21431dad5..94c1e15cc 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt @@ -1,7 +1,7 @@ package scientifik.kmath.structures import scientifik.kmath.operations.DoubleField -import scientifik.kmath.structures.NDArrays.simpleNDArray +import scientifik.kmath.structures.NDArrays.create import kotlin.test.Test import kotlin.test.assertEquals @@ -9,7 +9,7 @@ import kotlin.test.assertEquals class SimpleNDFieldTest{ @Test fun testStrides(){ - val ndArray = simpleNDArray(DoubleField, listOf(10,10)){(it[0]+it[1]).toDouble()} + val ndArray = create(DoubleField, listOf(10,10)){(it[0]+it[1]).toDouble()} assertEquals(ndArray[5,5], 10.0) } diff --git a/kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt b/kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt index 29e64a575..32d66a04b 100644 --- a/kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt +++ b/kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt @@ -5,4 +5,4 @@ import scientifik.kmath.operations.DoubleField /** * Using boxing implementation for js */ -actual val realNDFieldFactory: NDFieldFactory = NDArrays.createSimpleNDFieldFactory(DoubleField) \ No newline at end of file +actual val realNDFieldFactory: NDFieldFactory = NDArrays.createFactory(DoubleField) \ No newline at end of file diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt index c46d4210a..19d6c4041 100644 --- a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt +++ b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt @@ -7,12 +7,18 @@ private class RealNDField(shape: List) : BufferNDField(shape, Doubl override fun createBuffer(capacity: Int, initializer: (Int) -> Double): Buffer { val array = DoubleArray(capacity, initializer) val buffer = DoubleBuffer.wrap(array) - return object : Buffer { - override fun get(index: Int): Double = buffer.get(index) + return BufferOfDoubles(buffer) + } - override fun set(index: Int, value: Double) { - buffer.put(index, value) - } + private class BufferOfDoubles(val buffer: DoubleBuffer): Buffer{ + override fun get(index: Int): Double = buffer.get(index) + + override fun set(index: Int, value: Double) { + buffer.put(index, value) + } + + override fun copy(): Buffer { + return BufferOfDoubles(buffer) } } } diff --git a/settings.gradle b/settings.gradle index 03760221a..988a94843 100644 --- a/settings.gradle +++ b/settings.gradle @@ -10,4 +10,5 @@ enableFeaturePreview('GRADLE_METADATA') rootProject.name = 'kmath' include ':kmath-core' +include ':kmath-commons'