diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt index f605f6df3..8e429087a 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt @@ -2,6 +2,7 @@ package scientifik.kmath.histogram import scientifik.kmath.structures.ArrayBuffer import scientifik.kmath.structures.Buffer +import scientifik.kmath.structures.DoubleBuffer typealias Point = Buffer @@ -49,6 +50,8 @@ interface MutableHistogram>: Histogram{ fun MutableHistogram.put(vararg point: T) = put(ArrayBuffer(point)) +fun MutableHistogram.put(vararg point: Number) = put(DoubleBuffer(point.map { it.toDouble() }.toDoubleArray())) + fun MutableHistogram.fill(sequence: Iterable>) = sequence.forEach { put(it) } /** diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt index 9ea17c730..25a1bf33a 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt @@ -4,11 +4,11 @@ import scientifik.kmath.linear.Vector import scientifik.kmath.operations.Space import scientifik.kmath.structures.NDStructure -data class BinTemplate>(val center: Vector, val sizes: Vector) { +data class BinTemplate>(val center: Vector, val sizes: Vector) { fun contains(vector: Point): Boolean { if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}") - val upper = center + sizes/2.0 - val lower = center - sizes/2.0 + val upper = center + sizes / 2.0 + val lower = center - sizes / 2.0 return vector.asSequence().mapIndexed { i, value -> value in lower[i]..upper[i] }.all { it } diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt index 0b71405c6..223e63419 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt @@ -1,5 +1,7 @@ package scientifik.kmath.linear +import scientifik.kmath.operations.DoubleField +import scientifik.kmath.operations.Field import scientifik.kmath.structures.MutableNDStructure import scientifik.kmath.structures.NDStructure import scientifik.kmath.structures.genericNdStructure @@ -9,7 +11,7 @@ import kotlin.math.absoluteValue /** * Implementation based on Apache common-maths LU-decomposition */ -abstract class LUDecomposition>(val matrix: Matrix) { +abstract class LUDecomposition, F : Field>(val matrix: Matrix) { private val field get() = matrix.context.field /** Entries of LU decomposition. */ @@ -31,7 +33,7 @@ abstract class LUDecomposition>(val matrix: Matrix) { * L is a lower-triangular matrix * @return the L matrix (or null if decomposed matrix is singular) */ - val l: Matrix by lazy { + val l: Matrix by lazy { matrix.context.produce { i, j -> when { j < i -> lu[i, j] @@ -48,7 +50,7 @@ abstract class LUDecomposition>(val matrix: Matrix) { * U is an upper-triangular matrix * @return the U matrix (or null if decomposed matrix is singular) */ - val u: Matrix by lazy { + val u: Matrix by lazy { matrix.context.produce { i, j -> if (j >= i) lu[i, j] else field.zero } @@ -64,7 +66,7 @@ abstract class LUDecomposition>(val matrix: Matrix) { * @return the P rows permutation matrix (or null if decomposed matrix is singular) * @see .getPivot */ - val p: Matrix by lazy { + 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 @@ -181,7 +183,7 @@ abstract class LUDecomposition>(val matrix: Matrix) { } -class RealLUDecomposition(matrix: Matrix, private val singularityThreshold: Double = DEFAULT_TOO_SMALL) : LUDecomposition(matrix) { +class RealLUDecomposition(matrix: RealMatrix, private val singularityThreshold: Double = DEFAULT_TOO_SMALL) : LUDecomposition(matrix) { override fun isSingular(value: Double): Boolean { return value.absoluteValue < singularityThreshold } @@ -194,12 +196,12 @@ class RealLUDecomposition(matrix: Matrix, private val singularityThresho /** Specialized solver. */ -object RealLUSolver : LinearSolver { +object RealLUSolver : LinearSolver { - fun decompose(mat: Matrix, threshold: Double = 1e-11): RealLUDecomposition = RealLUDecomposition(mat, threshold) + fun decompose(mat: Matrix, threshold: Double = 1e-11): RealLUDecomposition = RealLUDecomposition(mat, threshold) - override fun solve(a: Matrix, b: Matrix): Matrix { + override fun solve(a: RealMatrix, b: RealMatrix): RealMatrix { val decomposition = decompose(a, a.context.field.zero) if (b.rows != a.rows) { diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt index 9e8692b42..acd27dbd9 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt @@ -1,184 +1,17 @@ package scientifik.kmath.linear -import scientifik.kmath.operations.* -import scientifik.kmath.structures.ExtendedNDField -import scientifik.kmath.structures.GenericNDField -import scientifik.kmath.structures.NDField - -/** - * 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 MatrixSpace(val rows: Int, val columns: Int, val field: Field) : Space> { - - /** - * Produce the element of this space - */ - abstract fun produce(initializer: (Int, Int) -> T): Matrix - - /** - * Produce new matrix space with given dimensions. The space produced could be raised from cache since [MatrixSpace] does not have mutable elements - */ - abstract fun produceSpace(rows: Int, columns: Int): MatrixSpace - - override val zero: Matrix by lazy { - produce { _, _ -> field.zero } - } - -// val one: Matrix by lazy { -// produce { i, j -> if (i == j) field.one else field.zero } -// } - - override fun add(a: Matrix, b: Matrix): Matrix { - return produce { i, j -> with(field) { a[i, j] + b[i, j] } } - } - - override fun multiply(a: Matrix, k: Double): Matrix { - //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: Matrix, b: Matrix): Matrix { - 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 until a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) } - } - } - - override fun equals(other: Any?): Boolean { - if (this === other) return true - if (other !is MatrixSpace<*>) return false - - if (rows != other.rows) return false - if (columns != other.columns) return false - if (field != other.field) return false - - return true - } - - override fun hashCode(): Int { - var result = rows - result = 31 * result + columns - result = 31 * result + field.hashCode() - return result - } -} - -infix fun Matrix.dot(b: Matrix): Matrix = this.context.multiply(this, b) - -/** - * A matrix-like structure - */ -interface Matrix : SpaceElement, MatrixSpace> { - /** - * 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 - - override val self: Matrix - get() = this - - fun transpose(): Matrix { - return object : Matrix { - override val context: MatrixSpace = 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] - } - } - - companion object { - - /** - * Create [ArrayMatrix] with custom field - */ - fun of(rows: Int, columns: Int, field: Field, initializer: (Int, Int) -> T) = - ArrayMatrix(ArrayMatrixSpace(rows, columns, field), initializer) - - /** - * Create [ArrayMatrix] of doubles. The implementation in general should be faster than generic one due to boxing. - */ - fun ofReal(rows: Int, columns: Int, initializer: (Int, Int) -> Double) = - ArrayMatrix(ArrayMatrixSpace(rows, columns, DoubleField, realNDFieldFactory), initializer) - - /** - * Create a diagonal value matrix. By default value equals [Field.one]. - */ - fun diagonal(rows: Int, columns: Int, field: Field, values: (Int) -> T = { field.one }): Matrix { - return of(rows, columns, field) { i, j -> if (i == j) values(i) else field.zero } - } - - /** - * Equality check on two generic matrices - */ - fun equals(mat1: Matrix<*>, mat2: Matrix<*>): Boolean { - if (mat1 === mat2) return true - if (mat1.context != mat2.context) return false - for (i in 0 until mat1.rows) { - for (j in 0 until mat2.columns) { - if (mat1[i, j] != mat2[i, j]) return false - } - } - return true - } - } -} - - - - -typealias NDFieldFactory = (IntArray) -> NDField - -internal fun genericNDFieldFactory(field: Field): NDFieldFactory = { index -> GenericNDField(index, field) } -internal val realNDFieldFactory: NDFieldFactory = { index -> ExtendedNDField(index, DoubleField) } - - -/** - * NDArray-based implementation of vector space. By default uses slow [GenericNDField], but could be overridden with custom [NDField] factory. - */ -class ArrayMatrixSpace( - rows: Int, - columns: Int, - field: Field, - val ndFactory: NDFieldFactory = genericNDFieldFactory(field) -) : MatrixSpace(rows, columns, field) { - - val ndField by lazy { - ndFactory(intArrayOf(rows, columns)) - } - - override fun produce(initializer: (Int, Int) -> T): Matrix = ArrayMatrix(this, initializer) - - override fun produceSpace(rows: Int, columns: Int): ArrayMatrixSpace { - return ArrayMatrixSpace(rows, columns, field, ndFactory) - } -} - +import scientifik.kmath.operations.DoubleField +import scientifik.kmath.operations.Field +import scientifik.kmath.operations.Norm /** * 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.toMatrix()).toVector() - fun inverse(a: Matrix): Matrix = solve(a, Matrix.diagonal(a.rows, a.columns, a.context.field)) +interface LinearSolver> { + fun solve(a: Matrix, b: Matrix): Matrix + fun solve(a: Matrix, b: Vector): Vector = solve(a, b.toMatrix()).toVector() + fun inverse(a: Matrix): Matrix = solve(a, Matrix.diagonal(a.rows, a.columns, a.context.field)) } /** @@ -192,7 +25,7 @@ fun List.toVector() = Vector.ofReal(this.size) { this[it] } /** * Convert matrix to vector if it is possible */ -fun Matrix.toVector(): Vector { +fun > Matrix.toVector(): Vector { return when { this.columns == 1 -> { // if (this is ArrayMatrix) { @@ -208,7 +41,7 @@ fun Matrix.toVector(): Vector { } } -fun Vector.toMatrix(): Matrix { +fun > Vector.toMatrix(): Matrix { // return if (this is ArrayVector) { // //Reuse existing underlying array // ArrayMatrix(ArrayMatrixSpace(size, 1, context.field, context.ndFactory), array) @@ -216,11 +49,14 @@ fun Vector.toMatrix(): Matrix { // //Generic vector // matrix(size, 1, context.field) { i, j -> get(i) } // } - return Matrix.of(size, 1, context.field) { i, _ -> get(i) } + return Matrix.of(size, 1, context.space) { i, _ -> get(i) } } -object VectorL2Norm: Norm, Double> { - override fun norm(arg: Vector): Double { +object VectorL2Norm : Norm, Double> { + override fun norm(arg: Vector): Double { return kotlin.math.sqrt(arg.sumByDouble { it.toDouble() }) } } + +typealias RealVector = Vector +typealias RealMatrix = Matrix diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt new file mode 100644 index 000000000..f7ae04849 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt @@ -0,0 +1,171 @@ +package scientifik.kmath.linear + +import scientifik.kmath.operations.* +import scientifik.kmath.structures.ExtendedNDField +import scientifik.kmath.structures.GenericNDField +import scientifik.kmath.structures.NDField + +/** + * 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 MatrixSpace>(val rows: Int, val columns: Int, val field: F) : Space> { + + /** + * Produce the element of this space + */ + abstract fun produce(initializer: (Int, Int) -> T): Matrix + + /** + * Produce new matrix space with given dimensions. The space produced could be raised from cache since [MatrixSpace] does not have mutable elements + */ + abstract fun produceSpace(rows: Int, columns: Int): MatrixSpace + + override val zero: Matrix by lazy { + produce { _, _ -> field.zero } + } + +// val one: Matrix by lazy { +// produce { i, j -> if (i == j) field.one else field.zero } +// } + + override fun add(a: Matrix, b: Matrix): Matrix { + return produce { i, j -> with(field) { a[i, j] + b[i, j] } } + } + + override fun multiply(a: Matrix, k: Double): Matrix { + //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: Matrix, b: Matrix): Matrix { + 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 until a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) } + } + } + + override fun equals(other: Any?): Boolean { + if (this === other) return true + if (other !is MatrixSpace<*,*>) return false + + if (rows != other.rows) return false + if (columns != other.columns) return false + if (field != other.field) return false + + return true + } + + override fun hashCode(): Int { + var result = rows + result = 31 * result + columns + result = 31 * result + field.hashCode() + return result + } +} + +infix fun > Matrix.dot(b: Matrix): Matrix = this.context.multiply(this, b) + +/** + * A matrix-like structure + */ +interface Matrix> : SpaceElement, MatrixSpace> { + /** + * 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 + + override val self: Matrix + get() = this + + fun transpose(): Matrix { + return object : Matrix { + override val context: MatrixSpace = 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] + } + } + + companion object { + + /** + * Create [ArrayMatrix] with custom field + */ + fun > of(rows: Int, columns: Int, field: F, initializer: (Int, Int) -> T) = + ArrayMatrix(ArrayMatrixSpace(rows, columns, field), initializer) + + /** + * Create [ArrayMatrix] of doubles. The implementation in general should be faster than generic one due to boxing. + */ + fun ofReal(rows: Int, columns: Int, initializer: (Int, Int) -> Double) = + ArrayMatrix(ArrayMatrixSpace(rows, columns, DoubleField, realNDFieldFactory), initializer) + + /** + * Create a diagonal value matrix. By default value equals [Field.one]. + */ + fun > diagonal(rows: Int, columns: Int, field: F, values: (Int) -> T = { field.one }): Matrix { + return of(rows, columns, field) { i, j -> if (i == j) values(i) else field.zero } + } + + /** + * Equality check on two generic matrices + */ + fun equals(mat1: Matrix<*, *>, mat2: Matrix<*, *>): Boolean { + if (mat1 === mat2) return true + if (mat1.context != mat2.context) return false + for (i in 0 until mat1.rows) { + for (j in 0 until mat2.columns) { + if (mat1[i, j] != mat2[i, j]) return false + } + } + return true + } + } +} + + + + +typealias NDFieldFactory = (IntArray) -> NDField + +internal fun > genericNDFieldFactory(field: F): NDFieldFactory = { index -> GenericNDField(index, field) } +internal val realNDFieldFactory: NDFieldFactory = { index -> ExtendedNDField(index, DoubleField) } + + +/** + * NDArray-based implementation of vector space. By default uses slow [GenericNDField], but could be overridden with custom [NDField] factory. + */ +class ArrayMatrixSpace>( + rows: Int, + columns: Int, + field: F, + val ndFactory: NDFieldFactory = genericNDFieldFactory(field) +) : MatrixSpace(rows, columns, field) { + + val ndField by lazy { + ndFactory(intArrayOf(rows, columns)) + } + + override fun produce(initializer: (Int, Int) -> T): Matrix = ArrayMatrix(this, initializer) + + override fun produceSpace(rows: Int, columns: Int): ArrayMatrixSpace { + return ArrayMatrixSpace(rows, columns, field, ndFactory) + } +} diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt index 70480542f..3e0f91ebe 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt @@ -12,34 +12,34 @@ import scientifik.kmath.structures.get * A linear space for vectors. * Could be used on any point-like structure */ -abstract class VectorSpace(val size: Int, val field: Field) : Space> { +abstract class VectorSpace>(val size: Int, val space: S) : Space> { - abstract fun produce(initializer: (Int) -> T): Vector + abstract fun produce(initializer: (Int) -> T): Vector - override val zero: Vector by lazy { produce { field.zero } } + override val zero: Vector by lazy { produce { space.zero } } - override fun add(a: Point, b: Point): Vector = produce { with(field) { a[it] + b[it] } } + override fun add(a: Point, b: Point): Vector = produce { with(space) { a[it] + b[it] } } - override fun multiply(a: Point, k: Double): Vector = produce { with(field) { a[it] * k } } + override fun multiply(a: Point, k: Double): Vector = produce { with(space) { a[it] * k } } } /** * A point coupled to the linear space */ -interface Vector : SpaceElement, VectorSpace>, Point, Iterable { +interface Vector> : SpaceElement, VectorSpace>, Point, Iterable { override val size: Int get() = context.size - override operator fun plus(b: Point): Vector = context.add(self, b) - override operator fun minus(b: Point): Vector = context.add(self, context.multiply(b, -1.0)) - override operator fun times(k: Number): Vector = context.multiply(self, k.toDouble()) - override operator fun div(k: Number): Vector = context.multiply(self, 1.0 / k.toDouble()) + override operator fun plus(b: Point): Vector = context.add(self, b) + override operator fun minus(b: Point): Vector = context.add(self, context.multiply(b, -1.0)) + override operator fun times(k: Number): Vector = context.multiply(self, k.toDouble()) + override operator fun div(k: Number): Vector = context.multiply(self, 1.0 / k.toDouble()) companion object { /** * Create vector with custom field */ - fun of(size: Int, field: Field, initializer: (Int) -> T) = + fun > of(size: Int, field: F, initializer: (Int) -> T) = ArrayVector(ArrayVectorSpace(size, field), initializer) /** @@ -50,7 +50,7 @@ interface Vector : SpaceElement, VectorSpace>, Point, It fun ofReal(vararg point: Double) = point.toVector() - fun equals(v1: Vector<*>, v2: Vector<*>): Boolean { + fun equals(v1: Vector<*,*>, v2: Vector<*,*>): Boolean { if (v1 === v2) return true if (v1.context != v2.context) return false for (i in 0 until v2.size) { @@ -61,24 +61,24 @@ interface Vector : SpaceElement, VectorSpace>, Point, It } } -class ArrayVectorSpace( +class ArrayVectorSpace>( size: Int, - field: Field, - val ndFactory: NDFieldFactory = genericNDFieldFactory(field) -) : VectorSpace(size, field) { + field: F, + val ndFactory: NDFieldFactory = genericNDFieldFactory(field) +) : VectorSpace(size, field) { val ndField by lazy { ndFactory(intArrayOf(size)) } - override fun produce(initializer: (Int) -> T): Vector = ArrayVector(this, initializer) + override fun produce(initializer: (Int) -> T): Vector = ArrayVector(this, initializer) } /** * Member of [ArrayMatrixSpace] which wraps 2-D array */ -class ArrayMatrix internal constructor(override val context: ArrayMatrixSpace, val element: NDElement) : Matrix { +class ArrayMatrix> internal constructor(override val context: ArrayMatrixSpace, val element: NDElement) : Matrix { - constructor(context: ArrayMatrixSpace, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) }) + constructor(context: ArrayMatrixSpace, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) }) override val rows: Int get() = context.rows @@ -88,13 +88,13 @@ class ArrayMatrix internal constructor(override val context: ArrayMatri return element[i, j] } - override val self: ArrayMatrix get() = this + override val self: ArrayMatrix get() = this } -class ArrayVector internal constructor(override val context: VectorSpace, val element: NDElement) : Vector { +class ArrayVector> internal constructor(override val context: VectorSpace, val element: NDElement) : Vector { - constructor(context: ArrayVectorSpace, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) }) + constructor(context: ArrayVectorSpace, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) }) init { if (context.size != element.shape[0]) { @@ -106,13 +106,12 @@ class ArrayVector internal constructor(override val context: VectorSpac return element[index] } - override val self: ArrayVector get() = this + override val self: ArrayVector get() = this override fun iterator(): Iterator = (0 until size).map { element[it] }.iterator() - override fun copy(): ArrayVector = ArrayVector(context, element) + override fun copy(): ArrayVector = ArrayVector(context, element) override fun toString(): String = this.joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() } } -typealias RealVector = Vector \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt index 08751d934..93ffad5a3 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt @@ -9,33 +9,33 @@ import scientifik.kmath.operations.TrigonometricOperations /** * NDField that supports [ExtendedField] operations on its elements */ -class ExtendedNDField(shape: IntArray, override val field: ExtendedField) : NDField(shape, field), - TrigonometricOperations>, - PowerOperations>, - ExponentialOperations> { +class ExtendedNDField>(shape: IntArray, field: F) : NDField(shape, field), + TrigonometricOperations>, + PowerOperations>, + ExponentialOperations> { - override fun produceStructure(initializer: (IntArray) -> N): NDStructure { - return genericNdStructure(shape, initializer) + override fun produceStructure(initializer: F.(IntArray) -> N): NDStructure { + return genericNdStructure(shape) { field.initializer(it) } } - override fun power(arg: NDElement, pow: Double): NDElement { - return arg.transform { d -> with(field){power(d,pow)} } + override fun power(arg: NDElement, pow: Double): NDElement { + return arg.transform { d -> with(field) { power(d, pow) } } } - override fun exp(arg: NDElement): NDElement { - return arg.transform { d -> with(field){exp(d)} } + override fun exp(arg: NDElement): NDElement { + return arg.transform { d -> with(field) { exp(d) } } } - override fun ln(arg: NDElement): NDElement { - return arg.transform { d -> with(field){ln(d)} } + override fun ln(arg: NDElement): NDElement { + return arg.transform { d -> with(field) { ln(d) } } } - override fun sin(arg: NDElement): NDElement { - return arg.transform { d -> with(field){sin(d)} } + override fun sin(arg: NDElement): NDElement { + return arg.transform { d -> with(field) { sin(d) } } } - override fun cos(arg: NDElement): NDElement { - return arg.transform { d -> with(field){cos(d)} } + override fun cos(arg: NDElement): NDElement { + return arg.transform { d -> with(field) { cos(d) } } } } diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt index 6d3ce05a5..d8e9cc31b 100644 --- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt +++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt @@ -15,24 +15,24 @@ class ShapeMismatchException(val expected: IntArray, val actual: IntArray) : Run * @param field - operations field defined on individual array element * @param T the type of the element contained in NDArray */ -abstract class NDField(val shape: IntArray, open val field: Field) : Field> { +abstract class NDField>(val shape: IntArray, val field: F) : Field> { - abstract fun produceStructure(initializer: (IntArray) -> T): NDStructure + abstract fun produceStructure(initializer: F.(IntArray) -> T): NDStructure /** * Create new instance of NDArray using field shape and given initializer * The producer takes list of indices as argument and returns contained value */ - fun produce(initializer: (IntArray) -> T): NDElement = NDElement(this, produceStructure(initializer)) + fun produce(initializer: F.(IntArray) -> T): NDElement = NDElement(this, produceStructure(initializer)) - override val zero: NDElement by lazy { - produce { this.field.zero } + override val zero: NDElement by lazy { + produce { zero } } /** * Check the shape of given NDArray and throw exception if it does not coincide with shape of the field */ - private fun checkShape(vararg elements: NDElement) { + private fun checkShape(vararg elements: NDElement) { elements.forEach { if (!shape.contentEquals(it.shape)) { throw ShapeMismatchException(shape, it.shape) @@ -43,7 +43,7 @@ abstract class NDField(val shape: IntArray, open val field: Field) : Field /** * Element-by-element addition */ - override fun add(a: NDElement, b: NDElement): NDElement { + override fun add(a: NDElement, b: NDElement): NDElement { checkShape(a, b) return produce { with(field) { a[it] + b[it] } } } @@ -51,18 +51,18 @@ abstract class NDField(val shape: IntArray, open val field: Field) : Field /** * Multiply all elements by cinstant */ - override fun multiply(a: NDElement, k: Double): NDElement { + override fun multiply(a: NDElement, k: Double): NDElement { checkShape(a) return produce { with(field) { a[it] * k } } } - override val one: NDElement - get() = produce { this.field.one } + override val one: NDElement + get() = produce { one } /** * Element-by-element multiplication */ - override fun multiply(a: NDElement, b: NDElement): NDElement { + override fun multiply(a: NDElement, b: NDElement): NDElement { checkShape(a) return produce { with(field) { a[it] * b[it] } } } @@ -70,65 +70,65 @@ abstract class NDField(val shape: IntArray, open val field: Field) : Field /** * Element-by-element division */ - override fun divide(a: NDElement, b: NDElement): NDElement { + override fun divide(a: NDElement, b: NDElement): NDElement { checkShape(a) return produce { with(field) { a[it] / b[it] } } } - /** - * Reverse sum operation - */ - operator fun T.plus(arg: NDElement): NDElement = arg + this - - /** - * Reverse minus operation - */ - operator fun T.minus(arg: NDElement): NDElement = arg.transform { _, value -> - with(arg.context.field) { - this@minus - value - } - } - - /** - * Reverse product operation - */ - operator fun T.times(arg: NDElement): NDElement = arg * this - - /** - * Reverse division operation - */ - operator fun T.div(arg: NDElement): NDElement = arg.transform { _, value -> - with(arg.context.field) { - this@div / value - } - } +// /** +// * Reverse sum operation +// */ +// operator fun T.plus(arg: NDElement): NDElement = arg + this +// +// /** +// * Reverse minus operation +// */ +// operator fun T.minus(arg: NDElement): NDElement = arg.transform { _, value -> +// with(arg.context.field) { +// this@minus - value +// } +// } +// +// /** +// * Reverse product operation +// */ +// operator fun T.times(arg: NDElement): NDElement = arg * this +// +// /** +// * Reverse division operation +// */ +// operator fun T.div(arg: NDElement): NDElement = arg.transform { _, value -> +// with(arg.context.field) { +// this@div / value +// } +// } } /** * Immutable [NDStructure] coupled to the context. Emulates Python ndarray */ -class NDElement(override val context: NDField, private val structure: NDStructure) : FieldElement, NDField>, NDStructure by structure { +class NDElement>(override val context: NDField, private val structure: NDStructure) : FieldElement, NDField>, NDStructure by structure { //TODO ensure structure is immutable - override val self: NDElement + override val self: NDElement get() = this - inline fun transform(crossinline action: (IntArray, T) -> T): NDElement = context.produce { action(it, get(*it)) } - inline fun transform(crossinline action: (T) -> T): NDElement = context.produce { action(get(*it)) } + inline fun transform(crossinline action: (IntArray, T) -> T): NDElement = context.produce { action(it, get(*it)) } + inline fun transform(crossinline action: (T) -> T): NDElement = context.produce { action(get(*it)) } } /** * Element by element application of any operation on elements to the whole array. Just like in numpy */ -operator fun Function1.invoke(ndElement: NDElement): NDElement = ndElement.transform { _, value -> this(value) } +operator fun > Function1.invoke(ndElement: NDElement): NDElement = ndElement.transform { _, value -> this(value) } /* plus and minus */ /** * Summation operation for [NDElement] and single element */ -operator fun NDElement.plus(arg: T): NDElement = transform { _, value -> +operator fun > NDElement.plus(arg: T): NDElement = transform { _, value -> with(context.field) { arg + value } @@ -137,7 +137,7 @@ operator fun NDElement.plus(arg: T): NDElement = transform { _, value /** * Subtraction operation between [NDElement] and single element */ -operator fun NDElement.minus(arg: T): NDElement = transform { _, value -> +operator fun > NDElement.minus(arg: T): NDElement = transform { _, value -> with(context.field) { arg - value } @@ -148,7 +148,7 @@ operator fun NDElement.minus(arg: T): NDElement = transform { _, value /** * Product operation for [NDElement] and single element */ -operator fun NDElement.times(arg: T): NDElement = transform { _, value -> +operator fun > NDElement.times(arg: T): NDElement = transform { _, value -> with(context.field) { arg * value } @@ -157,14 +157,14 @@ operator fun NDElement.times(arg: T): NDElement = transform { _, value /** * Division operation between [NDElement] and single element */ -operator fun NDElement.div(arg: T): NDElement = transform { _, value -> +operator fun > NDElement.div(arg: T): NDElement = transform { _, value -> with(context.field) { arg / value } } -class GenericNDField(shape: IntArray, field: Field) : NDField(shape, field) { - override fun produceStructure(initializer: (IntArray) -> T): NDStructure = genericNdStructure(shape, initializer) +class GenericNDField>(shape: IntArray, field: F) : NDField(shape, field) { + override fun produceStructure(initializer: F.(IntArray) -> T): NDStructure = genericNdStructure(shape) { field.initializer(it) } } //typealias NDFieldFactory = (IntArray)->NDField @@ -173,23 +173,24 @@ object NDArrays { /** * Create a platform-optimized NDArray of doubles */ - fun realNDArray(shape: IntArray, initializer: (IntArray) -> Double = { 0.0 }): NDElement { + fun realNDArray(shape: IntArray, initializer: DoubleField.(IntArray) -> Double = { 0.0 }): NDElement { return ExtendedNDField(shape, DoubleField).produce(initializer) } - fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDElement { + fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDElement { return realNDArray(intArrayOf(dim)) { initializer(it[0]) } } - fun real2DArray(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }): NDElement { + fun real2DArray(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }): NDElement { return realNDArray(intArrayOf(dim1, dim2)) { initializer(it[0], it[1]) } } - fun real3DArray(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }): NDElement { + fun real3DArray(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }): NDElement { return realNDArray(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) } } - inline fun produceReal(shape: IntArray, block: ExtendedNDField.() -> NDElement) = ExtendedNDField(shape, DoubleField).run(block) + inline fun produceReal(shape: IntArray, block: ExtendedNDField.() -> NDElement) = + ExtendedNDField(shape, DoubleField).run(block) // /** // * Simple boxing NDField @@ -199,7 +200,7 @@ object NDArrays { /** * Simple boxing NDArray */ - fun create(field: Field, shape: IntArray, initializer: (IntArray) -> T): NDElement { + fun > create(field: F, shape: IntArray, initializer: (IntArray) -> T): NDElement { return GenericNDField(shape, field).produce { initializer(it) } } } diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt index cd58a63b4..990adb0c7 100644 --- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt +++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt @@ -32,7 +32,7 @@ class NumberNDFieldTest { for (i in 0..2) { for (j in 0..2) { val expected = (i * 10 + j).toDouble() - assertEquals(expected, array[i, j],"Error at index [$i, $j]") + assertEquals(expected, array[i, j], "Error at index [$i, $j]") } } } @@ -41,25 +41,25 @@ class NumberNDFieldTest { fun testExternalFunction() { val function: (Double) -> Double = { x -> x.pow(2) + 2 * x + 1 } val result = function(array1) + 1.0 - assertEquals(10.0, result[1,1]) + assertEquals(10.0, result[1, 1]) } @Test fun testLibraryFunction() { val abs: (Double) -> Double = ::abs val result = abs(array2) - assertEquals(2.0, result[0,2]) + assertEquals(2.0, result[0, 2]) } - object L2Norm: Norm, Double> { - override fun norm(arg: NDElement): Double { + object L2Norm : Norm, Double> { + override fun norm(arg: NDElement): Double { return kotlin.math.sqrt(arg.sumByDouble { it.second.toDouble() }) } } @Test - fun testInternalContext(){ - produceReal(array1.shape){ + fun testInternalContext() { + produceReal(array1.shape) { with(L2Norm) { 1 + norm(array1) + exp(array2) }