diff --git a/CHANGELOG.md b/CHANGELOG.md index 4e1d9f530..bbfb68501 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,9 +3,11 @@ ## Unreleased ### Added -- Explicit `mutableStructureND` builders for mutable stucures +- New Attributes-kt module that could be used as stand-alone. It declares type-safe attributes containers. +- Explicit `mutableStructureND` builders for mutable structures ### Changed +- Features replaced with Attributes. ### Deprecated diff --git a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attribute.kt b/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attribute.kt index 4b9e27655..6fa142180 100644 --- a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attribute.kt +++ b/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attribute.kt @@ -9,20 +9,31 @@ import kotlin.reflect.KType public interface Attribute +/** + * An attribute that could be either present or absent + */ +public interface FlagAttribute : Attribute +/** + * An attribute with a default value + */ public interface AttributeWithDefault : Attribute { public val default: T } +/** + * Attribute containing a set of values + */ public interface SetAttribute : Attribute> /** * An attribute that has a type parameter for value + * @param type parameter-type */ -public abstract class PolymorphicAttribute(public val type: KType) : Attribute { - override fun equals(other: Any?): Boolean = (other as? PolymorphicAttribute<*>)?.type == this.type +public abstract class PolymorphicAttribute(public val type: SafeType) : Attribute { + override fun equals(other: Any?): Boolean = other != null && + (this::class == other::class) && + (other as? PolymorphicAttribute<*>)?.type == this.type - override fun hashCode(): Int { - return type.hashCode() - } + override fun hashCode(): Int = this::class.hashCode() + type.hashCode() } diff --git a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attributes.kt b/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attributes.kt index 25767a686..f93b6446a 100644 --- a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attributes.kt +++ b/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attributes.kt @@ -24,18 +24,45 @@ public value class Attributes internal constructor(public val content: Map Attributes.getOrDefault(attribute: AttributeWithDefault): T = get(attribute) ?: attribute.default -public fun > Attributes.withAttribute( +/** + * Check if there is an attribute that matches given key by type and adheres to [predicate]. + */ +@Suppress("UNCHECKED_CAST") +public inline fun > Attributes.any(predicate: (value: T) -> Boolean): Boolean = + content.any { (mapKey, mapValue) -> mapKey is A && predicate(mapValue as T) } + +/** + * Check if there is an attribute of given type (subtypes included) + */ +public inline fun > Attributes.any(): Boolean = + content.any { (mapKey, _) -> mapKey is A } + +/** + * Check if [Attributes] contains a flag. Multiple keys that are instances of a flag could be present + */ +public inline fun Attributes.has(): Boolean = + content.keys.any { it is A } + +/** + * Create [Attributes] with an added or replaced attribute key. + */ +public fun > Attributes.withAttribute( attribute: A, - attrValue: T?, -): Attributes = Attributes( - if (attrValue == null) { - content - attribute - } else { - content + (attribute to attrValue) - } -) + attrValue: T, +): Attributes = Attributes(content + (attribute to attrValue)) + +public fun > Attributes.withAttribute(attribute: A): Attributes = + withAttribute(attribute, Unit) + +/** + * Create new [Attributes] by removing [attribute] key + */ +public fun Attributes.withoutAttribute(attribute: Attribute<*>): Attributes = Attributes(content.minus(attribute)) /** * Add an element to a [SetAttribute] @@ -63,6 +90,9 @@ public fun > Attributes.withoutAttributeElement( ) } +/** + * Create [Attributes] with a single key + */ public fun > Attributes( attribute: A, attrValue: T, diff --git a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/SafeType.kt b/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/SafeType.kt new file mode 100644 index 000000000..44ead89cb --- /dev/null +++ b/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/SafeType.kt @@ -0,0 +1,26 @@ +/* + * Copyright 2018-2023 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.attributes + +import kotlin.reflect.KClass +import kotlin.reflect.KType +import kotlin.reflect.typeOf + +/** + * Safe variant ok Kotlin [KType] that ensures that the type parameter is of the same type ask [kType] + * + * @param kType raw [KType] + */ +public class SafeType @PublishedApi internal constructor(public val kType: KType) + +public inline fun safeTypeOf(): SafeType = SafeType(typeOf()) + +/** + * Derive Kotlin [KClass] from this type and fail if the type is not a class (should not happen) + */ +@Suppress("UNCHECKED_CAST") +@UnstableAttributesAPI +public val SafeType.kClass: KClass get() = kType.classifier as KClass \ No newline at end of file diff --git a/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/annotations.kt b/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/annotations.kt new file mode 100644 index 000000000..80d93daaf --- /dev/null +++ b/attributes-kt/src/commonMain/kotlin/space/kscience/attributes/annotations.kt @@ -0,0 +1,17 @@ +/* + * Copyright 2018-2023 KMath contributors. + * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. + */ + +package space.kscience.attributes + +/** + * Marks declarations that are still experimental in the Attributes-kt APIs, which means that the design of the corresponding + * declarations has open issues that may (or may not) lead to their changes in the future. Roughly speaking, there is + * a chance of those declarations will be deprecated in the future or the semantics of their behavior may change + * in some way that may break some code. + */ +@MustBeDocumented +@Retention(value = AnnotationRetention.BINARY) +@RequiresOptIn("This API is unstable and could change in future", RequiresOptIn.Level.WARNING) +public annotation class UnstableAttributesAPI \ No newline at end of file diff --git a/buildSrc/settings.gradle.kts b/buildSrc/settings.gradle.kts index 38eebeecc..53c47a814 100644 --- a/buildSrc/settings.gradle.kts +++ b/buildSrc/settings.gradle.kts @@ -1,4 +1,4 @@ -enableFeaturePreview("TYPESAFE_PROJECT_ACCESSORS") +rootProject.name = "buildSrc" dependencyResolutionManagement { val projectProperties = java.util.Properties() @@ -13,6 +13,7 @@ dependencyResolutionManagement { val toolsVersion: String = projectProperties["toolsVersion"].toString() + @Suppress("UnstableApiUsage") repositories { mavenLocal() maven("https://repo.kotlin.link") diff --git a/buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt b/buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt index d973ebae4..22e86c2a2 100644 --- a/buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt +++ b/buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt @@ -56,6 +56,8 @@ public object EjmlLinearSpace${ops} : EjmlLinearSpace<${type}, ${kmathAlgebra}, */ override val elementAlgebra: $kmathAlgebra get() = $kmathAlgebra + override val elementType: KType get() = typeOf<$type>() + @Suppress("UNCHECKED_CAST") override fun Matrix<${type}>.toEjml(): Ejml${type}Matrix<${ejmlMatrixType}> = when { this is Ejml${type}Matrix<*> && origin is $ejmlMatrixType -> this as Ejml${type}Matrix<${ejmlMatrixType}> diff --git a/gradle.properties b/gradle.properties index 162720167..2f1a2a030 100644 --- a/gradle.properties +++ b/gradle.properties @@ -11,5 +11,4 @@ toolsVersion=0.14.9-kotlin-1.8.20 org.gradle.parallel=true org.gradle.workers.max=4 org.gradle.configureondemand=true -org.gradle.jvmargs=-Xmx4096m - +org.gradle.jvmargs=-Xmx4096m \ No newline at end of file diff --git a/gradle/wrapper/gradle-wrapper.properties b/gradle/wrapper/gradle-wrapper.properties index fae08049a..15de90249 100644 --- a/gradle/wrapper/gradle-wrapper.properties +++ b/gradle/wrapper/gradle-wrapper.properties @@ -1,5 +1,5 @@ distributionBase=GRADLE_USER_HOME distributionPath=wrapper/dists -distributionUrl=https\://services.gradle.org/distributions/gradle-8.1.1-bin.zip +distributionUrl=https\://services.gradle.org/distributions/gradle-8.2-bin.zip zipStoreBase=GRADLE_USER_HOME zipStorePath=wrapper/dists diff --git a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt index e369effdf..85ff881ef 100644 --- a/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt +++ b/kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt @@ -8,12 +8,14 @@ package space.kscience.kmath.commons.linear import org.apache.commons.math3.linear.* import space.kscience.kmath.UnstableKMathAPI import space.kscience.kmath.linear.* -import space.kscience.kmath.nd.StructureFeature +import space.kscience.kmath.nd.StructureAttribute import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.DoubleBuffer import kotlin.reflect.KClass +import kotlin.reflect.KType import kotlin.reflect.cast +import kotlin.reflect.typeOf public class CMMatrix(public val origin: RealMatrix) : Matrix { override val rowNum: Int get() = origin.rowDimension @@ -38,6 +40,8 @@ public fun RealVector.toPoint(): CMVector = CMVector(this) public object CMLinearSpace : LinearSpace { override val elementAlgebra: DoubleField get() = DoubleField + override val elementType: KType = typeOf() + override fun buildMatrix( rows: Int, columns: Int, @@ -99,7 +103,7 @@ public object CMLinearSpace : LinearSpace { v * this @UnstableKMathAPI - override fun computeFeature(structure: Matrix, type: KClass): F? { + override fun computeFeature(structure: Matrix, type: KClass): F? { //Return the feature if it is intrinsic to the structure structure.getFeature(type)?.let { return it } @@ -108,36 +112,37 @@ public object CMLinearSpace : LinearSpace { return when (type) { IsDiagonal::class -> if (origin is DiagonalMatrix) IsDiagonal else null - DeterminantFeature::class, LupDecompositionFeature::class -> object : - DeterminantFeature, - LupDecompositionFeature { + Determinant::class, LupDecompositionAttribute::class -> object : + Determinant, + LupDecompositionAttribute { private val lup by lazy { LUDecomposition(origin) } override val determinant: Double by lazy { lup.determinant } - override val l: Matrix by lazy> { CMMatrix(lup.l).withFeature(LFeature) } - override val u: Matrix by lazy> { CMMatrix(lup.u).withFeature(UFeature) } + override val l: Matrix by lazy> { CMMatrix(lup.l).withAttribute(LowerTriangular) } + override val u: Matrix by lazy> { CMMatrix(lup.u).withAttribute(UpperTriangular) } override val p: Matrix by lazy { CMMatrix(lup.p) } } - CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature { + CholeskyDecompositionAttribute::class -> object : CholeskyDecompositionAttribute { override val l: Matrix by lazy> { val cholesky = CholeskyDecomposition(origin) - CMMatrix(cholesky.l).withFeature(LFeature) + CMMatrix(cholesky.l).withAttribute(LowerTriangular) } } - QRDecompositionFeature::class -> object : QRDecompositionFeature { + QRDecompositionAttribute::class -> object : QRDecompositionAttribute { private val qr by lazy { QRDecomposition(origin) } - override val q: Matrix by lazy> { CMMatrix(qr.q).withFeature(OrthogonalFeature) } - override val r: Matrix by lazy> { CMMatrix(qr.r).withFeature(UFeature) } + override val q: Matrix by lazy> { CMMatrix(qr.q).withAttribute(OrthogonalAttribute) } + override val r: Matrix by lazy> { CMMatrix(qr.r).withAttribute(UpperTriangular) } } - SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature { + SingularValueDecompositionAttribute::class -> object : SingularValueDecompositionAttribute { private val sv by lazy { SingularValueDecomposition(origin) } override val u: Matrix by lazy { CMMatrix(sv.u) } override val s: Matrix by lazy { CMMatrix(sv.s) } override val v: Matrix by lazy { CMMatrix(sv.v) } override val singularValues: Point by lazy { DoubleBuffer(sv.singularValues) } } + else -> null }?.let(type::cast) } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSAlgebra.kt index 8c7cb0cf1..28546ab98 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSAlgebra.kt @@ -13,6 +13,8 @@ import space.kscience.kmath.structures.MutableBufferFactory import space.kscience.kmath.structures.asBuffer import kotlin.math.max import kotlin.math.min +import kotlin.reflect.KType +import kotlin.reflect.typeOf /** * Class representing both the value and the differentials of a function. diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/BufferedLinearSpace.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/BufferedLinearSpace.kt index 4bba47a91..6bfe7581b 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/BufferedLinearSpace.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/BufferedLinearSpace.kt @@ -5,12 +5,15 @@ package space.kscience.kmath.linear +import space.kscience.attributes.SafeType import space.kscience.kmath.PerformancePitfall import space.kscience.kmath.nd.* import space.kscience.kmath.operations.* import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.VirtualBuffer import space.kscience.kmath.structures.indices +import kotlin.reflect.KType +import kotlin.reflect.typeOf public class BufferedLinearSpace>( diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/DoubleLinearSpace.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/DoubleLinearSpace.kt index 940af4a86..226ededce 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/DoubleLinearSpace.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/DoubleLinearSpace.kt @@ -12,6 +12,8 @@ import space.kscience.kmath.operations.DoubleField import space.kscience.kmath.operations.invoke import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.DoubleBuffer +import kotlin.reflect.KType +import kotlin.reflect.typeOf public object DoubleLinearSpace : LinearSpace { @@ -20,7 +22,7 @@ public object DoubleLinearSpace : LinearSpace { override fun buildMatrix( rows: Int, columns: Int, - initializer: DoubleField.(i: Int, j: Int) -> Double + initializer: DoubleField.(i: Int, j: Int) -> Double, ): Matrix = DoubleFieldOpsND.structureND(ShapeND(rows, columns)) { (i, j) -> DoubleField.initializer(i, j) }.as2D() diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt index a82bafe57..b83316e84 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt @@ -5,16 +5,15 @@ package space.kscience.kmath.linear +import space.kscience.attributes.SafeType import space.kscience.kmath.UnstableKMathAPI -import space.kscience.kmath.nd.MutableStructure2D -import space.kscience.kmath.nd.Structure2D -import space.kscience.kmath.nd.StructureFeature -import space.kscience.kmath.nd.as1D +import space.kscience.kmath.nd.* import space.kscience.kmath.operations.BufferRingOps import space.kscience.kmath.operations.Ring import space.kscience.kmath.operations.invoke import space.kscience.kmath.structures.Buffer -import kotlin.reflect.KClass +import kotlin.reflect.KType +import kotlin.reflect.typeOf /** * Alias for [Structure2D] with more familiar name. @@ -31,13 +30,19 @@ public typealias MutableMatrix = MutableStructure2D */ public typealias Point = Buffer +/** + * A marker interface for algebras that operate on matrices + * @param T type of matrix element + */ +public interface MatrixOperations + /** * Basic operations on matrices and vectors. * * @param T the type of items in the matrices. * @param A the type of ring over [T]. */ -public interface LinearSpace> { +public interface LinearSpace> : MatrixOperations { public val elementAlgebra: A /** @@ -167,16 +172,16 @@ public interface LinearSpace> { public operator fun T.times(v: Point): Point = v * this /** - * Compute a feature of the structure in this scope. Structure features take precedence other context features. + * Get an attribute value for the structure in this scope. Structure features take precedence other context features. * - * @param F the type of feature. + * @param A the type of feature. * @param structure the structure. - * @param type the [KClass] instance of [F]. + * @param attribute to be computed * @return a feature object or `null` if it isn't present. */ @UnstableKMathAPI - public fun computeFeature(structure: Matrix, type: KClass): F? = - structure.getFeature(type) + public fun > attributeFor(structure: StructureND<*>, attribute: A): T? = + structure.attributes[attribute] public companion object { @@ -184,23 +189,12 @@ public interface LinearSpace> { * A structured matrix with custom buffer */ public fun > buffered( - algebra: A + algebra: A, ): LinearSpace = BufferedLinearSpace(BufferRingOps(algebra)) } } -/** - * Get a feature of the structure in this scope. Structure features take precedence other context features. - * - * @param T the type of items in the matrices. - * @param F the type of feature. - * @return a feature object or `null` if it isn't present. - */ -@UnstableKMathAPI -public inline fun LinearSpace.computeFeature(structure: Matrix): F? = - computeFeature(structure, F::class) - public inline operator fun , R> LS.invoke(block: LS.() -> R): R = run(block) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt index 650e7be5c..c44564231 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt @@ -3,67 +3,92 @@ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. */ +@file:Suppress("UnusedReceiverParameter") + package space.kscience.kmath.linear +import space.kscience.attributes.PolymorphicAttribute +import space.kscience.attributes.SafeType +import space.kscience.attributes.safeTypeOf import space.kscience.kmath.UnstableKMathAPI import space.kscience.kmath.operations.* -import space.kscience.kmath.structures.BufferAccessor2D -import space.kscience.kmath.structures.DoubleBuffer -import space.kscience.kmath.structures.MutableBuffer -import space.kscience.kmath.structures.MutableBufferFactory +import space.kscience.kmath.structures.* /** - * Common implementation of [LupDecompositionFeature]. + * Matrices with this feature support LU factorization with partial pivoting: *[p] · a = [l] · [u]* where + * *a* is the owning matrix. + * + * @param T the type of matrices' items. + * @param l The lower triangular matrix in this decomposition. It may have [LowerTriangular]. + * @param u The upper triangular matrix in this decomposition. It may have [UpperTriangular]. + * @param p he permutation matrix in this decomposition. May have [Determinant] attribute */ -public class LupDecomposition( - public val context: LinearSpace, - public val elementContext: Field, - public val lu: Matrix, - public val pivot: IntArray, - private val even: Boolean, -) : LupDecompositionFeature, DeterminantFeature { - /** - * Returns the matrix L of the decomposition. - * - * L is a lower-triangular matrix with [Ring.one] in diagonal - */ - override val l: Matrix = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j -> - when { - j < i -> lu[i, j] - j == i -> elementContext.one - else -> elementContext.zero - } - }.withFeature(LFeature) +public data class LupDecomposition( + public val l: Matrix, + public val u: Matrix, + public val p: Matrix, +) - /** - * Returns the matrix U of the decomposition. - * - * U is an upper-triangular matrix including the diagonal - */ - override val u: Matrix = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j -> - if (j >= i) lu[i, j] else elementContext.zero - }.withFeature(UFeature) +public class LupDecompositionAttribute(type: SafeType>) : + PolymorphicAttribute>(type), + MatrixAttribute> - /** - * Returns the P rows permutation matrix. - * - * P is a sparse matrix with exactly one element set to [Ring.one] in - * each row and each column, all other elements being set to [Ring.zero]. - */ - override val p: Matrix = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j -> - if (j == pivot[i]) elementContext.one else elementContext.zero - } +public val MatrixOperations.LUP: LupDecompositionAttribute + get() = LupDecompositionAttribute(safeTypeOf()) - /** - * Return the determinant of the matrix - * @return determinant of the matrix - */ - override val determinant: T by lazy { - elementContext { (0 until lu.shape[0]).fold(if (even) one else -one) { value, i -> value * lu[i, i] } } - } -} +///** +// * Common implementation of [LupDecomposition]. +// */ +//private class LupDecompositionImpl( +// public val elementContext: Field, +// public val lu: Matrix, +// public val pivot: IntBuffer, +// private val even: Boolean, +//) : LupDecomposition { +// /** +// * Returns the matrix L of the decomposition. +// * +// * L is a lower-triangular matrix with [Ring.one] in diagonal +// */ +// override val l: Matrix = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j -> +// when { +// j < i -> lu[i, j] +// j == i -> elementContext.one +// else -> elementContext.zero +// } +// }.withFeature(LowerTriangular) +// +// +// /** +// * Returns the matrix U of the decomposition. +// * +// * U is an upper-triangular matrix including the diagonal +// */ +// override val u: Matrix = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j -> +// if (j >= i) lu[i, j] else elementContext.zero +// }.withFeature(UpperTriangular) +// +// /** +// * Returns the P rows permutation matrix. +// * +// * P is a sparse matrix with exactly one element set to [Ring.one] in +// * each row and each column, all other elements being set to [Ring.zero]. +// */ +// override val p: Matrix = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j -> +// if (j == pivot[i]) elementContext.one else elementContext.zero +// } +// +// /** +// * Return the determinant of the matrix +// * @return determinant of the matrix +// */ +// override val determinant: T by lazy { +// elementContext { (0 until lu.shape[0]).fold(if(even) one else -one) { value, i -> value * lu[i, i] } } +// } +// +//} @PublishedApi internal fun > LinearSpace>.abs(value: T): T = @@ -73,7 +98,6 @@ internal fun > LinearSpace>.abs(value: T): T = * Create a lup decomposition of generic matrix. */ public fun > LinearSpace>.lup( - factory: MutableBufferFactory, matrix: Matrix, checkSingular: (T) -> Boolean, ): LupDecomposition { @@ -82,15 +106,15 @@ public fun > LinearSpace>.lup( val pivot = IntArray(matrix.rowNum) //TODO just waits for multi-receivers - BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run { + BufferAccessor2D(matrix.rowNum, matrix.colNum, elementAlgebra.bufferFactory).run { elementAlgebra { val lu = create(matrix) - // Initialize permutation array and parity + // Initialize the permutation array and parity for (row in 0 until m) pivot[row] = row var even = true - // Initialize permutation array and parity + // Initialize the permutation array and parity for (row in 0 until m) pivot[row] = row // Loop over columns @@ -145,46 +169,57 @@ public fun > LinearSpace>.lup( for (row in col + 1 until m) lu[row, col] /= luDiag } - return LupDecomposition(this@lup, elementAlgebra, lu.collect(), pivot, even) + val l: MatrixWrapper = VirtualMatrix(rowNum, colNum) { i, j -> + when { + j < i -> lu[i, j] + j == i -> one + else -> zero + } + }.withAttribute(LowerTriangular) + + val u = VirtualMatrix(rowNum, colNum) { i, j -> + if (j >= i) lu[i, j] else zero + }.withAttribute(UpperTriangular) + + val p = VirtualMatrix(rowNum, colNum) { i, j -> + if (j == pivot[i]) one else zero + }.withAttribute(Determinant, if (even) one else -one) + + return LupDecomposition(l, u, p) } } } -public inline fun > LinearSpace>.lup( - matrix: Matrix, - noinline checkSingular: (T) -> Boolean, -): LupDecomposition = lup(MutableBuffer.Companion::auto, matrix, checkSingular) public fun LinearSpace.lup( matrix: Matrix, singularityThreshold: Double = 1e-11, -): LupDecomposition = - lup(::DoubleBuffer, matrix) { it < singularityThreshold } +): LupDecomposition = lup(matrix) { it < singularityThreshold } -internal fun LupDecomposition.solve( - factory: MutableBufferFactory, +internal fun > LinearSpace.solve( + lup: LupDecomposition, matrix: Matrix, ): Matrix { - require(matrix.rowNum == pivot.size) { "Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}" } + require(matrix.rowNum == lup.l.rowNum) { "Matrix dimension mismatch. Expected ${lup.l.rowNum}, but got ${matrix.colNum}" } - BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run { - elementContext { + BufferAccessor2D(matrix.rowNum, matrix.colNum, elementAlgebra.bufferFactory).run { + elementAlgebra { // Apply permutations to b val bp = create { _, _ -> zero } - for (row in pivot.indices) { + for (row in 0 until rowNum) { val bpRow = bp.row(row) val pRow = pivot[row] for (col in 0 until matrix.colNum) bpRow[col] = matrix[pRow, col] } // Solve LY = b - for (col in pivot.indices) { + for (col in 0 until colNum) { val bpCol = bp.row(col) - for (i in col + 1 until pivot.size) { + for (i in col + 1 until colNum) { val bpI = bp.row(i) - val luICol = lu[i, col] + val luICol = lup.l[i, col] for (j in 0 until matrix.colNum) { bpI[j] -= bpCol[j] * luICol } @@ -192,19 +227,19 @@ internal fun LupDecomposition.solve( } // Solve UX = Y - for (col in pivot.size - 1 downTo 0) { + for (col in colNum - 1 downTo 0) { val bpCol = bp.row(col) - val luDiag = lu[col, col] + val luDiag = lup.u[col, col] for (j in 0 until matrix.colNum) bpCol[j] /= luDiag for (i in 0 until col) { val bpI = bp.row(i) - val luICol = lu[i, col] + val luICol = lup.u[i, col] for (j in 0 until matrix.colNum) bpI[j] -= bpCol[j] * luICol } } - return context.buildMatrix(pivot.size, matrix.colNum) { i, j -> bp[i, j] } + return buildMatrix(matrix.rowNum, matrix.colNum) { i, j -> bp[i, j] } } } } @@ -214,17 +249,16 @@ internal fun LupDecomposition.solve( */ @OptIn(UnstableKMathAPI::class) public fun , F : Field> LinearSpace.lupSolver( - bufferFactory: MutableBufferFactory, singularityCheck: (T) -> Boolean, ): LinearSolver = object : LinearSolver { override fun solve(a: Matrix, b: Matrix): Matrix { // Use existing decomposition if it is provided by matrix - val decomposition = computeFeature(a) ?: lup(bufferFactory, a, singularityCheck) - return decomposition.solve(bufferFactory, b) + val decomposition = attributeFor(a, LUP) ?: lup(a, singularityCheck) + return solve(decomposition, b) } override fun inverse(matrix: Matrix): Matrix = solve(matrix, one(matrix.rowNum, matrix.colNum)) } public fun LinearSpace.lupSolver(singularityThreshold: Double = 1e-11): LinearSolver = - lupSolver(::DoubleBuffer) { it < singularityThreshold } + lupSolver { it < singularityThreshold } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt index 4d2f01e68..7b6c36dfc 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixBuilder.kt @@ -5,6 +5,7 @@ package space.kscience.kmath.linear +import space.kscience.attributes.FlagAttribute import space.kscience.kmath.UnstableKMathAPI import space.kscience.kmath.operations.Ring import space.kscience.kmath.structures.BufferAccessor2D @@ -49,10 +50,10 @@ public inline fun LinearSpace>.column( public fun LinearSpace>.column(vararg values: T): Matrix = column(values.size, values::get) -public object SymmetricMatrixFeature : MatrixFeature +public object Symmetric : MatrixAttribute, FlagAttribute /** - * Naive implementation of a symmetric matrix builder, that adds a [SymmetricMatrixFeature] tag. The resulting matrix contains + * Naive implementation of a symmetric matrix builder, that adds a [Symmetric] tag. The resulting matrix contains * full `size^2` number of elements, but caches elements during calls to save [builder] calls. [builder] is always called in the * upper triangle region meaning that `i <= j` */ @@ -72,6 +73,6 @@ public fun > MatrixBuilder.symmetric( } else { cached } - }.withFeature(SymmetricMatrixFeature) + }.withAttribute(Symmetric) } } \ No newline at end of file diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt index 62325b39d..b494625c7 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt @@ -3,20 +3,27 @@ * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file. */ +@file:OptIn(UnstableKMathAPI::class) +@file:Suppress("UnusedReceiverParameter") + package space.kscience.kmath.linear -import space.kscience.attributes.Attribute +import space.kscience.attributes.* +import space.kscience.kmath.UnstableKMathAPI +import space.kscience.kmath.nd.StructureAttribute +import kotlin.reflect.KType +import kotlin.reflect.typeOf /** * 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 : Attribute +public interface MatrixAttribute : StructureAttribute /** * Matrices with this feature are considered to have only diagonal non-zero elements. */ -public interface IsDiagonal : MatrixFeature { +public interface IsDiagonal : MatrixAttribute, FlagAttribute { public companion object : IsDiagonal } @@ -31,130 +38,110 @@ public object IsZero : IsDiagonal public object IsUnit : IsDiagonal /** - * Matrices with this feature can be inverted: *[inverse] = a−1* where *a* is the owning matrix. + * Matrices with this feature can be inverted. * * @param T the type of matrices' items. */ -public class Inverted private constructor() : MatrixFeature> { - internal val instance: Inverted = Inverted() -} +public class Inverted(type: SafeType>) : + PolymorphicAttribute>(type), + MatrixAttribute> -@Suppress("UNCHECKED_CAST") -public val LinearSpace.Inverted: Inverted get() = Inverted.instance as Inverted +public val MatrixOperations.Inverted: Inverted get() = Inverted(safeTypeOf()) /** * Matrices with this feature can compute their determinant. * * @param T the type of matrices' items. */ -public class DeterminantFeature : MatrixFeature +public class Determinant(type: SafeType) : + PolymorphicAttribute(type), + MatrixAttribute -/** - * Produces a [DeterminantFeature] where the [DeterminantFeature.determinant] is [determinant]. - * - * @param determinant the value of determinant. - * @return a new [DeterminantFeature]. - */ -@Suppress("FunctionName") -public fun DeterminantFeature(determinant: T): DeterminantFeature = object : DeterminantFeature { - override val determinant: T = determinant -} +public inline val MatrixOperations.Determinant: Determinant get() = Determinant(safeTypeOf()) /** * Matrices with this feature are lower triangular ones. */ -public object LFeature : MatrixFeature +public object LowerTriangular : MatrixAttribute, FlagAttribute /** * Matrices with this feature are upper triangular ones. */ -public object UFeature : MatrixFeature +public object UpperTriangular : MatrixAttribute, FlagAttribute + +/** + * Matrices with this feature support LU factorization: *a = [l] · [u]* where *a* is the owning matrix. + * @param l The lower triangular matrix in this decomposition. It may have [LowerTriangular]. + * @param u The upper triangular matrix in this decomposition. It may have [UpperTriangular]. + */ +public data class LUDecomposition(val l: Matrix, val u: Matrix) /** * Matrices with this feature support LU factorization: *a = [l] · [u]* where *a* is the owning matrix. * * @param T the type of matrices' items. */ -public interface LUDecompositionFeature : MatrixFeature { - /** - * The lower triangular matrix in this decomposition. It may have [LFeature]. - */ - public val l: Matrix +public class LuDecompositionAttribute(type: SafeType>) : + PolymorphicAttribute>(type), + MatrixAttribute> - /** - * The upper triangular matrix in this decomposition. It may have [UFeature]. - */ - public val u: Matrix -} +public val MatrixOperations.LU: LuDecompositionAttribute get() = LuDecompositionAttribute(safeTypeOf()) -/** - * Matrices with this feature support LU factorization with partial pivoting: *[p] · a = [l] · [u]* where - * *a* is the owning matrix. - * - * @param T the type of matrices' items. - */ -public interface LupDecompositionFeature : MatrixFeature { - /** - * The lower triangular matrix in this decomposition. It may have [LFeature]. - */ - public val l: Matrix - - /** - * The upper triangular matrix in this decomposition. It may have [UFeature]. - */ - public val u: Matrix - - /** - * The permutation matrix in this decomposition. - */ - public val p: Matrix -} /** * Matrices with this feature are orthogonal ones: *a · aT = u* where *a* is the owning matrix, *u* * is the unit matrix ([IsUnit]). */ -public object OrthogonalFeature : MatrixFeature +public object OrthogonalAttribute : MatrixAttribute, FlagAttribute -/** - * Matrices with this feature support QR factorization: *a = [q] · [r]* where *a* is the owning matrix. - * - * @param T the type of matrices' items. - */ -public interface QRDecompositionFeature : MatrixFeature { + +public interface QRDecomposition { /** - * The orthogonal matrix in this decomposition. It may have [OrthogonalFeature]. + * The orthogonal matrix in this decomposition. It may have [OrthogonalAttribute]. */ public val q: Matrix /** - * The upper triangular matrix in this decomposition. It may have [UFeature]. + * The upper triangular matrix in this decomposition. It may have [UpperTriangular]. */ public val r: Matrix } +/** + * Matrices with this feature support QR factorization: *a = [QR.q] · [QR.r]* where *a* is the owning matrix. + * + * @param T the type of matrices' items. + */ +public class QRDecompositionAttribute(type: SafeType>) : + PolymorphicAttribute>(type), + MatrixAttribute> + +public val MatrixOperations.QR: QRDecompositionAttribute + get() = QRDecompositionAttribute(safeTypeOf()) + +public interface CholeskyDecomposition { + /** + * The triangular matrix in this decomposition. It may have either [UpperTriangular] or [LowerTriangular]. + */ + public val l: Matrix +} + /** * Matrices with this feature support Cholesky factorization: *a = [l] · [l]H* where *a* is the * owning matrix. * * @param T the type of matrices' items. */ -public interface CholeskyDecompositionFeature : MatrixFeature { - /** - * The triangular matrix in this decomposition. It may have either [UFeature] or [LFeature]. - */ - public val l: Matrix -} +public class CholeskyDecompositionAttribute(type: SafeType>) : + PolymorphicAttribute>(type), + MatrixAttribute> -/** - * Matrices with this feature support SVD: *a = [u] · [s] · [v]H* where *a* is the owning - * matrix. - * - * @param T the type of matrices' items. - */ -public interface SingularValueDecompositionFeature : MatrixFeature { +public val MatrixOperations.Cholesky: CholeskyDecompositionAttribute + get() = CholeskyDecompositionAttribute(safeTypeOf()) + +public interface SingularValueDecomposition { /** - * The matrix in this decomposition. It is unitary, and it consists from left singular vectors. + * The matrix in this decomposition. It is unitary, and it consists of left singular vectors. */ public val u: Matrix @@ -164,14 +151,27 @@ public interface SingularValueDecompositionFeature : MatrixFeature public val s: Matrix /** - * The matrix in this decomposition. It is unitary, and it consists from right singular vectors. + * The matrix in this decomposition. It is unitary, and it consists of right singular vectors. */ public val v: Matrix /** - * The buffer of singular values of this SVD. + * The buffer of singular values for this SVD. */ public val singularValues: Point } +/** + * Matrices with this feature support SVD: *a = [u] · [s] · [v]H* where *a* is the owning + * matrix. + * + * @param T the type of matrices' items. + */ +public class SingularValueDecompositionAttribute(type: SafeType>) : + PolymorphicAttribute>(type), + MatrixAttribute> + +public val MatrixOperations.SVD: SingularValueDecompositionAttribute + get() = SingularValueDecompositionAttribute(safeTypeOf()) + //TODO add sparse matrix feature diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt index b510cc697..0f0c9ce3f 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt @@ -5,19 +5,20 @@ package space.kscience.kmath.linear +import space.kscience.attributes.Attribute import space.kscience.attributes.Attributes +import space.kscience.attributes.withAttribute import space.kscience.kmath.UnstableKMathAPI -import space.kscience.kmath.misc.FeatureSet import space.kscience.kmath.operations.Ring /** - * A [Matrix] that holds [MatrixFeature] objects. + * A [Matrix] that holds [MatrixAttribute] objects. * * @param T the type of items. */ public class MatrixWrapper internal constructor( public val origin: Matrix, - public val attributes: Attributes, + override val attributes: Attributes, ) : Matrix by origin { override fun toString(): String = "MatrixWrapper(matrix=$origin, features=$attributes)" @@ -34,23 +35,31 @@ public val Matrix.origin: Matrix /** * Add a single feature to a [Matrix] */ -public fun Matrix.withFeature(newFeature: MatrixFeature): MatrixWrapper = if (this is MatrixWrapper) { - MatrixWrapper(origin, attributes.with(newFeature)) +public fun > Matrix.withAttribute( + attribute: A, + attrValue: T, +): MatrixWrapper = if (this is MatrixWrapper) { + MatrixWrapper(origin, attributes.withAttribute(attribute,attrValue)) } else { - MatrixWrapper(this, FeatureSet.of(newFeature)) + MatrixWrapper(this, Attributes(attribute, attrValue)) } -@Deprecated("To be replaced by withFeature") -public operator fun Matrix.plus(newFeature: MatrixFeature): MatrixWrapper = withFeature(newFeature) +public fun > Matrix.withAttribute( + attribute: A, +): MatrixWrapper = if (this is MatrixWrapper) { + MatrixWrapper(origin, attributes.withAttribute(attribute)) +} else { + MatrixWrapper(this, Attributes(attribute, Unit)) +} /** - * Add a collection of features to a [Matrix] + * Modify matrix attributes */ -public fun Matrix.withFeatures(newFeatures: Iterable): MatrixWrapper = +public fun Matrix.modifyAttributes(modifier: (Attributes) -> Attributes): MatrixWrapper = if (this is MatrixWrapper) { - MatrixWrapper(origin, attributes.with(newFeatures)) + MatrixWrapper(origin, modifier(attributes)) } else { - MatrixWrapper(this, FeatureSet.of(newFeatures)) + MatrixWrapper(this, modifier(Attributes.EMPTY)) } /** @@ -59,9 +68,9 @@ public fun Matrix.withFeatures(newFeatures: Iterable public fun LinearSpace>.one( rows: Int, columns: Int, -): Matrix = VirtualMatrix(rows, columns) { i, j -> +): MatrixWrapper = VirtualMatrix(rows, columns) { i, j -> if (i == j) elementAlgebra.one else elementAlgebra.zero -}.withFeature(IsUnit) +}.withAttribute(IsUnit) /** @@ -70,16 +79,16 @@ public fun LinearSpace>.one( public fun LinearSpace>.zero( rows: Int, columns: Int, -): Matrix = VirtualMatrix(rows, columns) { _, _ -> +): MatrixWrapper = VirtualMatrix(rows, columns) { _, _ -> elementAlgebra.zero -}.withFeature(IsZero) +}.withAttribute(IsZero) -public class TransposedFeature(public val original: Matrix) : MatrixFeature +public class TransposedAttribute(public val original: Matrix) : MatrixAttribute /** * Create a virtual transposed matrix without copying anything. `A.transpose().transpose() === A` */ @Suppress("UNCHECKED_CAST") @OptIn(UnstableKMathAPI::class) -public fun Matrix.transpose(): Matrix = getFeature(TransposedFeature::class)?.original as? Matrix - ?: VirtualMatrix(colNum, rowNum) { i, j -> get(j, i) }.withFeature(TransposedFeature(this)) +public fun Matrix.transpose(): Matrix = getFeature(TransposedAttribute::class)?.original as? Matrix + ?: VirtualMatrix(colNum, rowNum) { i, j -> get(j, i) }.withAttribute(TransposedAttribute(this)) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/VirtualMatrix.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/VirtualMatrix.kt index 55b970f4a..7927e4dba 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/VirtualMatrix.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/VirtualMatrix.kt @@ -5,6 +5,7 @@ package space.kscience.kmath.linear +import space.kscience.attributes.Attributes import space.kscience.kmath.nd.ShapeND @@ -16,6 +17,7 @@ import space.kscience.kmath.nd.ShapeND public class VirtualMatrix( override val rowNum: Int, override val colNum: Int, + override val attributes: Attributes = Attributes.EMPTY, public val generator: (i: Int, j: Int) -> T, ) : Matrix { @@ -24,5 +26,8 @@ public class VirtualMatrix( override operator fun get(i: Int, j: Int): T = generator(i, j) } -public fun MatrixBuilder.virtual(generator: (i: Int, j: Int) -> T): VirtualMatrix = - VirtualMatrix(rows, columns, generator) +public fun MatrixBuilder.virtual( + attributes: Attributes = Attributes.EMPTY, + generator: (i: Int, j: Int) -> T, +): VirtualMatrix = + VirtualMatrix(rows, columns, attributes, generator) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt index a75807bad..3a9a1833a 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt @@ -79,7 +79,7 @@ public interface AlgebraND> : Algebra> { * @return a feature object or `null` if it isn't present. */ @UnstableKMathAPI - public fun getFeature(structure: StructureND, type: KClass): F? = + public fun getFeature(structure: StructureND, type: KClass): F? = structure.getFeature(type) public companion object @@ -93,7 +93,7 @@ public interface AlgebraND> : Algebra> { * @return a feature object or `null` if it isn't present. */ @UnstableKMathAPI -public inline fun AlgebraND.getFeature(structure: StructureND): F? = +public inline fun AlgebraND.getFeature(structure: StructureND): F? = getFeature(structure, F::class) /** diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt index 9217f6fdc..72ecc9ef0 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt @@ -34,9 +34,9 @@ public open class BufferND( /** * Create a generic [BufferND] using provided [initializer] */ -public fun BufferND( +public inline fun BufferND( shape: ShapeND, - bufferFactory: BufferFactory = BufferFactory.boxing(), + bufferFactory: BufferFactory = BufferFactory.auto(), initializer: (IntArray) -> T, ): BufferND { val strides = Strides(shape) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt index e006d09eb..931730399 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt @@ -110,7 +110,7 @@ private value class Structure2DWrapper(val structure: StructureND) : S @PerformancePitfall override operator fun get(i: Int, j: Int): T = structure[i, j] - override fun getFeature(type: KClass): F? = structure.getFeature(type) + override fun getFeature(type: KClass): F? = structure.getFeature(type) @PerformancePitfall override fun elements(): Sequence> = structure.elements() diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt index e9814acbf..724cc6f69 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt @@ -7,6 +7,7 @@ package space.kscience.kmath.nd import space.kscience.attributes.Attribute import space.kscience.attributes.AttributeContainer +import space.kscience.attributes.SafeType import space.kscience.kmath.PerformancePitfall import space.kscience.kmath.linear.LinearSpace import space.kscience.kmath.operations.Ring @@ -15,9 +16,8 @@ import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.BufferFactory import kotlin.jvm.JvmName import kotlin.math.abs -import kotlin.reflect.KClass -public interface StructureFeature : Attribute +public interface StructureAttribute : Attribute /** * Represents n-dimensional structure i.e., multidimensional container of items of the same type and size. The number @@ -122,9 +122,14 @@ public interface StructureND : AttributeContainer, WithShape { */ public fun buffered( strides: Strides, - bufferFactory: BufferFactory = BufferFactory.boxing(), initializer: (IntArray) -> T, - ): BufferND = BufferND(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) }) + ): BufferND = BufferND(strides, Buffer.boxing(strides.linearSize) { i -> initializer(strides.index(i)) }) + + + public fun buffered( + shape: ShapeND, + initializer: (IntArray) -> T, + ): BufferND = buffered(ColumnStrides(shape), initializer) /** * Inline create NDStructure with non-boxing buffer implementation if it is possible @@ -135,16 +140,11 @@ public interface StructureND : AttributeContainer, WithShape { ): BufferND = BufferND(strides, Buffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) }) public inline fun auto( - type: KClass, + type: SafeType, strides: Strides, crossinline initializer: (IntArray) -> T, ): BufferND = BufferND(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) }) - public fun buffered( - shape: ShapeND, - bufferFactory: BufferFactory = BufferFactory.boxing(), - initializer: (IntArray) -> T, - ): BufferND = buffered(ColumnStrides(shape), bufferFactory, initializer) public inline fun auto( shape: ShapeND, @@ -159,7 +159,7 @@ public interface StructureND : AttributeContainer, WithShape { auto(ColumnStrides(ShapeND(shape)), initializer) public inline fun auto( - type: KClass, + type: SafeType, vararg shape: Int, crossinline initializer: (IntArray) -> T, ): BufferND = auto(type, ColumnStrides(ShapeND(shape)), initializer) diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt index 0960ab023..e590dcc3d 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt @@ -9,14 +9,14 @@ import space.kscience.kmath.UnstableKMathAPI import space.kscience.kmath.expressions.Symbol import space.kscience.kmath.operations.Ring.Companion.optimizedPower import space.kscience.kmath.structures.MutableBufferFactory +import kotlin.reflect.KType /** * Represents an algebraic structure. * - * @param T the type of element of this structure. + * @param T the type of element which Algebra operates on. */ public interface Algebra { - /** * Provide a factory for buffers, associated with this [Algebra] */ @@ -67,12 +67,12 @@ public interface Algebra { * * @param operation the name of operation. * @param arg the argument of operation. - * @return a result of operation. + * @return the result of the operation. */ public fun unaryOperation(operation: String, arg: T): T = unaryOperationFunction(operation)(arg) /** - * Dynamically dispatches a binary operation with the certain name. + * Dynamically dispatches a binary operation with a certain name. * * Implementations must fulfil the following requirements: * @@ -87,7 +87,7 @@ public interface Algebra { error("Binary operation '$operation' not defined in $this") /** - * Dynamically invokes a binary operation with the certain name. + * Dynamically invokes a binary operation with a certain name. * * Implementations must fulfil the following requirements: * diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt index 26a4b0a53..9047feb29 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt @@ -8,6 +8,7 @@ package space.kscience.kmath.operations import space.kscience.kmath.structures.Buffer import space.kscience.kmath.structures.MutableBuffer import space.kscience.kmath.structures.MutableBufferFactory +import kotlin.reflect.KType public interface WithSize { public val size: Int @@ -18,6 +19,9 @@ public interface WithSize { */ public interface BufferAlgebra> : Algebra> { public val elementAlgebra: A + + public val elementType: KType + public val elementBufferFactory: MutableBufferFactory get() = elementAlgebra.bufferFactory public fun buffer(size: Int, vararg elements: T): Buffer { diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt index cef8d1d4d..d9ee25321 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt @@ -5,10 +5,12 @@ package space.kscience.kmath.structures +import space.kscience.attributes.SafeType +import space.kscience.attributes.safeTypeOf import space.kscience.kmath.operations.WithSize import space.kscience.kmath.operations.asSequence import kotlin.jvm.JvmInline -import kotlin.reflect.KClass +import kotlin.reflect.typeOf /** * Function that produces [Buffer] from its size and function that supplies values. @@ -99,13 +101,13 @@ public interface Buffer : WithSize { * The [size] is specified, and each element is calculated by calling the specified [initializer] function. */ @Suppress("UNCHECKED_CAST") - public inline fun auto(type: KClass, size: Int, initializer: (Int) -> T): Buffer = - when (type) { - Double::class -> MutableBuffer.double(size) { initializer(it) as Double } as Buffer - Short::class -> MutableBuffer.short(size) { initializer(it) as Short } as Buffer - Int::class -> MutableBuffer.int(size) { initializer(it) as Int } as Buffer - Long::class -> MutableBuffer.long(size) { initializer(it) as Long } as Buffer - Float::class -> MutableBuffer.float(size) { initializer(it) as Float } as Buffer + public inline fun auto(type: SafeType, size: Int, initializer: (Int) -> T): Buffer = + when (type.kType) { + typeOf() -> MutableBuffer.double(size) { initializer(it) as Double } as Buffer + typeOf() -> MutableBuffer.short(size) { initializer(it) as Short } as Buffer + typeOf() -> MutableBuffer.int(size) { initializer(it) as Int } as Buffer + typeOf() -> MutableBuffer.long(size) { initializer(it) as Long } as Buffer + typeOf() -> MutableBuffer.float(size) { initializer(it) as Float } as Buffer else -> boxing(size, initializer) } @@ -115,8 +117,8 @@ public interface Buffer : WithSize { * * The [size] is specified, and each element is calculated by calling the specified [initializer] function. */ - public inline fun auto(size: Int, initializer: (Int) -> T): Buffer = - auto(T::class, size, initializer) + public inline fun auto(size: Int, initializer: (Int) -> T): Buffer = + auto(safeTypeOf(), size, initializer) } } diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt index 48f3e919b..3edd7d3b1 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt @@ -27,12 +27,9 @@ internal class BufferAccessor2D( fun create(mat: Structure2D): MutableBuffer = create { i, j -> mat[i, j] } //TODO optimize wrapper - fun MutableBuffer.collect(): Structure2D = StructureND.buffered( - ColumnStrides(ShapeND(rowNum, colNum)), - factory - ) { (i, j) -> - get(i, j) - }.as2D() + fun MutableBuffer.toStructure2D(): Structure2D = StructureND.buffered( + ColumnStrides(ShapeND(rowNum, colNum)) + ) { (i, j) -> get(i, j) }.as2D() inner class Row(val buffer: MutableBuffer, val rowIndex: Int) : MutableBuffer { override val size: Int get() = colNum diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MutableBuffer.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MutableBuffer.kt index c0bfc6ecc..48b05e412 100644 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MutableBuffer.kt +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/MutableBuffer.kt @@ -5,7 +5,9 @@ package space.kscience.kmath.structures -import kotlin.reflect.KClass +import space.kscience.attributes.SafeType +import space.kscience.attributes.safeTypeOf +import kotlin.reflect.typeOf /** * A generic mutable random-access structure for both primitives and objects. @@ -74,13 +76,13 @@ public interface MutableBuffer : Buffer { * The [size] is specified, and each element is calculated by calling the specified [initializer] function. */ @Suppress("UNCHECKED_CAST") - public inline fun auto(type: KClass, size: Int, initializer: (Int) -> T): MutableBuffer = - when (type) { - Double::class -> double(size) { initializer(it) as Double } as MutableBuffer - Short::class -> short(size) { initializer(it) as Short } as MutableBuffer - Int::class -> int(size) { initializer(it) as Int } as MutableBuffer - Float::class -> float(size) { initializer(it) as Float } as MutableBuffer - Long::class -> long(size) { initializer(it) as Long } as MutableBuffer + public inline fun auto(type: SafeType, size: Int, initializer: (Int) -> T): MutableBuffer = + when (type.kType) { + typeOf() -> double(size) { initializer(it) as Double } as MutableBuffer + typeOf() -> short(size) { initializer(it) as Short } as MutableBuffer + typeOf() -> int(size) { initializer(it) as Int } as MutableBuffer + typeOf() -> float(size) { initializer(it) as Float } as MutableBuffer + typeOf() -> long(size) { initializer(it) as Long } as MutableBuffer else -> boxing(size, initializer) } @@ -90,11 +92,10 @@ public interface MutableBuffer : Buffer { * * The [size] is specified, and each element is calculated by calling the specified [initializer] function. */ - @Suppress("UNCHECKED_CAST") - public inline fun auto(size: Int, initializer: (Int) -> T): MutableBuffer = - auto(T::class, size, initializer) + public inline fun auto(size: Int, initializer: (Int) -> T): MutableBuffer = + auto(safeTypeOf(), size, initializer) } } -public sealed interface PrimitiveBuffer: MutableBuffer \ No newline at end of file +public sealed interface PrimitiveBuffer : MutableBuffer \ No newline at end of file diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt index 8925fb045..ec1883128 100644 --- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt +++ b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt @@ -12,6 +12,8 @@ import space.kscience.kmath.linear.Matrix import space.kscience.kmath.linear.Point import space.kscience.kmath.nd.Structure2D import space.kscience.kmath.operations.Ring +import kotlin.reflect.KType +import kotlin.reflect.typeOf /** * [LinearSpace] implementation specialized for a certain EJML type. @@ -43,5 +45,5 @@ public abstract class EjmlLinearSpace, out M : org.ejml @Suppress("UNCHECKED_CAST") @UnstableKMathAPI public fun EjmlMatrix.inverse(): Structure2D = - computeFeature(this, InverseMatrixFeature::class)?.inverse as Structure2D + attributeFor(this, InverseMatrixFeature::class)?.inverse as Structure2D } diff --git a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt b/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt deleted file mode 100644 index c56583fa8..000000000 --- a/kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt +++ /dev/null @@ -1,1003 +0,0 @@ -/* - * Copyright 2018-2021 KMath contributors. - * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file. - */ - -/* This file is generated with buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt */ - -package space.kscience.kmath.ejml - -import org.ejml.data.* -import org.ejml.dense.row.CommonOps_DDRM -import org.ejml.dense.row.CommonOps_FDRM -import org.ejml.dense.row.factory.DecompositionFactory_DDRM -import org.ejml.dense.row.factory.DecompositionFactory_FDRM -import org.ejml.sparse.FillReducing -import org.ejml.sparse.csc.CommonOps_DSCC -import org.ejml.sparse.csc.CommonOps_FSCC -import org.ejml.sparse.csc.factory.DecompositionFactory_DSCC -import org.ejml.sparse.csc.factory.DecompositionFactory_FSCC -import org.ejml.sparse.csc.factory.LinearSolverFactory_DSCC -import org.ejml.sparse.csc.factory.LinearSolverFactory_FSCC -import space.kscience.kmath.UnstableKMathAPI -import space.kscience.kmath.linear.* -import space.kscience.kmath.linear.Matrix -import space.kscience.kmath.nd.StructureFeature -import space.kscience.kmath.operations.DoubleField -import space.kscience.kmath.operations.FloatField -import space.kscience.kmath.operations.invoke -import space.kscience.kmath.structures.DoubleBuffer -import space.kscience.kmath.structures.FloatBuffer -import kotlin.reflect.KClass -import kotlin.reflect.cast - -/** - * [EjmlVector] specialization for [Double]. - */ -public class EjmlDoubleVector(override val origin: M) : EjmlVector(origin) { - init { - require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" } - } - - override operator fun get(index: Int): Double = origin[0, index] -} - -/** - * [EjmlVector] specialization for [Float]. - */ -public class EjmlFloatVector(override val origin: M) : EjmlVector(origin) { - init { - require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" } - } - - override operator fun get(index: Int): Float = origin[0, index] -} - -/** - * [EjmlMatrix] specialization for [Double]. - */ -public class EjmlDoubleMatrix(override val origin: M) : EjmlMatrix(origin) { - override operator fun get(i: Int, j: Int): Double = origin[i, j] -} - -/** - * [EjmlMatrix] specialization for [Float]. - */ -public class EjmlFloatMatrix(override val origin: M) : EjmlMatrix(origin) { - override operator fun get(i: Int, j: Int): Float = origin[i, j] -} - -/** - * [EjmlLinearSpace] implementation based on [CommonOps_DDRM], [DecompositionFactory_DDRM] operations and - * [DMatrixRMaj] matrices. - */ -public object EjmlLinearSpaceDDRM : EjmlLinearSpace() { - /** - * The [DoubleField] reference. - */ - override val elementAlgebra: DoubleField get() = DoubleField - - @Suppress("UNCHECKED_CAST") - override fun Matrix.toEjml(): EjmlDoubleMatrix = when { - this is EjmlDoubleMatrix<*> && origin is DMatrixRMaj -> this as EjmlDoubleMatrix - else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) } - } - - @Suppress("UNCHECKED_CAST") - override fun Point.toEjml(): EjmlDoubleVector = when { - this is EjmlDoubleVector<*> && origin is DMatrixRMaj -> this as EjmlDoubleVector - else -> EjmlDoubleVector(DMatrixRMaj(size, 1).also { - (0 until it.numRows).forEach { row -> it[row, 0] = get(row) } - }) - } - - override fun buildMatrix( - rows: Int, - columns: Int, - initializer: DoubleField.(i: Int, j: Int) -> Double, - ): EjmlDoubleMatrix = DMatrixRMaj(rows, columns).also { - (0 until rows).forEach { row -> - (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) } - } - }.wrapMatrix() - - override fun buildVector( - size: Int, - initializer: DoubleField.(Int) -> Double, - ): EjmlDoubleVector = EjmlDoubleVector(DMatrixRMaj(size, 1).also { - (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) } - }) - - private fun T.wrapMatrix() = EjmlDoubleMatrix(this) - private fun T.wrapVector() = EjmlDoubleVector(this) - - override fun Matrix.unaryMinus(): Matrix = this * elementAlgebra { -one } - - override fun Matrix.dot(other: Matrix): EjmlDoubleMatrix { - val out = DMatrixRMaj(1, 1) - CommonOps_DDRM.mult(toEjml().origin, other.toEjml().origin, out) - return out.wrapMatrix() - } - - override fun Matrix.dot(vector: Point): EjmlDoubleVector { - val out = DMatrixRMaj(1, 1) - CommonOps_DDRM.mult(toEjml().origin, vector.toEjml().origin, out) - return out.wrapVector() - } - - override operator fun Matrix.minus(other: Matrix): EjmlDoubleMatrix { - val out = DMatrixRMaj(1, 1) - - CommonOps_DDRM.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra { -one }, - other.toEjml().origin, - out, - ) - - return out.wrapMatrix() - } - - override operator fun Matrix.times(value: Double): EjmlDoubleMatrix { - val res = DMatrixRMaj(1, 1) - CommonOps_DDRM.scale(value, toEjml().origin, res) - return res.wrapMatrix() - } - - override fun Point.unaryMinus(): EjmlDoubleVector { - val res = DMatrixRMaj(1, 1) - CommonOps_DDRM.changeSign(toEjml().origin, res) - return res.wrapVector() - } - - override fun Matrix.plus(other: Matrix): EjmlDoubleMatrix { - val out = DMatrixRMaj(1, 1) - - CommonOps_DDRM.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra.one, - other.toEjml().origin, - out, - ) - - return out.wrapMatrix() - } - - override fun Point.plus(other: Point): EjmlDoubleVector { - val out = DMatrixRMaj(1, 1) - - CommonOps_DDRM.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra.one, - other.toEjml().origin, - out, - ) - - return out.wrapVector() - } - - override fun Point.minus(other: Point): EjmlDoubleVector { - val out = DMatrixRMaj(1, 1) - - CommonOps_DDRM.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra { -one }, - other.toEjml().origin, - out, - ) - - return out.wrapVector() - } - - override fun Double.times(m: Matrix): EjmlDoubleMatrix = m * this - - override fun Point.times(value: Double): EjmlDoubleVector { - val res = DMatrixRMaj(1, 1) - CommonOps_DDRM.scale(value, toEjml().origin, res) - return res.wrapVector() - } - - override fun Double.times(v: Point): EjmlDoubleVector = v * this - - @UnstableKMathAPI - override fun computeFeature(structure: Matrix, type: KClass): F? { - structure.getFeature(type)?.let { return it } - val origin = structure.toEjml().origin - - return when (type) { - InverseMatrixFeature::class -> object : InverseMatrixFeature { - override val inverse: Matrix by lazy { - val res = origin.copy() - CommonOps_DDRM.invert(res) - res.wrapMatrix() - } - } - - DeterminantFeature::class -> object : DeterminantFeature { - override val determinant: Double by lazy { CommonOps_DDRM.det(origin) } - } - - SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature { - private val svd by lazy { - DecompositionFactory_DDRM.svd(origin.numRows, origin.numCols, true, true, false) - .apply { decompose(origin.copy()) } - } - - override val u: Matrix by lazy { svd.getU(null, false).wrapMatrix() } - override val s: Matrix by lazy { svd.getW(null).wrapMatrix() } - override val v: Matrix by lazy { svd.getV(null, false).wrapMatrix() } - override val singularValues: Point by lazy { DoubleBuffer(svd.singularValues) } - } - - QRDecompositionFeature::class -> object : QRDecompositionFeature { - private val qr by lazy { - DecompositionFactory_DDRM.qr().apply { decompose(origin.copy()) } - } - - override val q: Matrix by lazy { - qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature) - } - - override val r: Matrix by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) } - } - - CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature { - override val l: Matrix by lazy { - val cholesky = - DecompositionFactory_DDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) } - - cholesky.getT(null).wrapMatrix().withFeature(LFeature) - } - } - - LupDecompositionFeature::class -> object : LupDecompositionFeature { - private val lup by lazy { - DecompositionFactory_DDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) } - } - - override val l: Matrix by lazy { - lup.getLower(null).wrapMatrix().withFeature(LFeature) - } - - override val u: Matrix by lazy { - lup.getUpper(null).wrapMatrix().withFeature(UFeature) - } - - override val p: Matrix by lazy { lup.getRowPivot(null).wrapMatrix() } - } - - else -> null - }?.let{ - type.cast(it) - } - } - - /** - * Solves for *x* in the following equation: *x = [a] -1 · [b]*. - * - * @param a the base matrix. - * @param b n by p matrix. - * @return the solution for *x* that is n by p. - */ - public fun solve(a: Matrix, b: Matrix): EjmlDoubleMatrix { - val res = DMatrixRMaj(1, 1) - CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res) - return res.wrapMatrix() - } - - /** - * Solves for *x* in the following equation: *x = [a] -1 · [b]*. - * - * @param a the base matrix. - * @param b n by p vector. - * @return the solution for *x* that is n by p. - */ - public fun solve(a: Matrix, b: Point): EjmlDoubleVector { - val res = DMatrixRMaj(1, 1) - CommonOps_DDRM.solve(DMatrixRMaj(a.toEjml().origin), DMatrixRMaj(b.toEjml().origin), res) - return EjmlDoubleVector(res) - } -} - -/** - * [EjmlLinearSpace] implementation based on [CommonOps_FDRM], [DecompositionFactory_FDRM] operations and - * [FMatrixRMaj] matrices. - */ -public object EjmlLinearSpaceFDRM : EjmlLinearSpace() { - /** - * The [FloatField] reference. - */ - override val elementAlgebra: FloatField get() = FloatField - - @Suppress("UNCHECKED_CAST") - override fun Matrix.toEjml(): EjmlFloatMatrix = when { - this is EjmlFloatMatrix<*> && origin is FMatrixRMaj -> this as EjmlFloatMatrix - else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) } - } - - @Suppress("UNCHECKED_CAST") - override fun Point.toEjml(): EjmlFloatVector = when { - this is EjmlFloatVector<*> && origin is FMatrixRMaj -> this as EjmlFloatVector - else -> EjmlFloatVector(FMatrixRMaj(size, 1).also { - (0 until it.numRows).forEach { row -> it[row, 0] = get(row) } - }) - } - - override fun buildMatrix( - rows: Int, - columns: Int, - initializer: FloatField.(i: Int, j: Int) -> Float, - ): EjmlFloatMatrix = FMatrixRMaj(rows, columns).also { - (0 until rows).forEach { row -> - (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) } - } - }.wrapMatrix() - - override fun buildVector( - size: Int, - initializer: FloatField.(Int) -> Float, - ): EjmlFloatVector = EjmlFloatVector(FMatrixRMaj(size, 1).also { - (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) } - }) - - private fun T.wrapMatrix() = EjmlFloatMatrix(this) - private fun T.wrapVector() = EjmlFloatVector(this) - - override fun Matrix.unaryMinus(): Matrix = this * elementAlgebra { -one } - - override fun Matrix.dot(other: Matrix): EjmlFloatMatrix { - val out = FMatrixRMaj(1, 1) - CommonOps_FDRM.mult(toEjml().origin, other.toEjml().origin, out) - return out.wrapMatrix() - } - - override fun Matrix.dot(vector: Point): EjmlFloatVector { - val out = FMatrixRMaj(1, 1) - CommonOps_FDRM.mult(toEjml().origin, vector.toEjml().origin, out) - return out.wrapVector() - } - - override operator fun Matrix.minus(other: Matrix): EjmlFloatMatrix { - val out = FMatrixRMaj(1, 1) - - CommonOps_FDRM.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra { -one }, - other.toEjml().origin, - out, - ) - - return out.wrapMatrix() - } - - override operator fun Matrix.times(value: Float): EjmlFloatMatrix { - val res = FMatrixRMaj(1, 1) - CommonOps_FDRM.scale(value, toEjml().origin, res) - return res.wrapMatrix() - } - - override fun Point.unaryMinus(): EjmlFloatVector { - val res = FMatrixRMaj(1, 1) - CommonOps_FDRM.changeSign(toEjml().origin, res) - return res.wrapVector() - } - - override fun Matrix.plus(other: Matrix): EjmlFloatMatrix { - val out = FMatrixRMaj(1, 1) - - CommonOps_FDRM.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra.one, - other.toEjml().origin, - out, - ) - - return out.wrapMatrix() - } - - override fun Point.plus(other: Point): EjmlFloatVector { - val out = FMatrixRMaj(1, 1) - - CommonOps_FDRM.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra.one, - other.toEjml().origin, - out, - ) - - return out.wrapVector() - } - - override fun Point.minus(other: Point): EjmlFloatVector { - val out = FMatrixRMaj(1, 1) - - CommonOps_FDRM.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra { -one }, - other.toEjml().origin, - out, - ) - - return out.wrapVector() - } - - override fun Float.times(m: Matrix): EjmlFloatMatrix = m * this - - override fun Point.times(value: Float): EjmlFloatVector { - val res = FMatrixRMaj(1, 1) - CommonOps_FDRM.scale(value, toEjml().origin, res) - return res.wrapVector() - } - - override fun Float.times(v: Point): EjmlFloatVector = v * this - - @UnstableKMathAPI - override fun computeFeature(structure: Matrix, type: KClass): F? { - structure.getFeature(type)?.let { return it } - val origin = structure.toEjml().origin - - return when (type) { - InverseMatrixFeature::class -> object : InverseMatrixFeature { - override val inverse: Matrix by lazy { - val res = origin.copy() - CommonOps_FDRM.invert(res) - res.wrapMatrix() - } - } - - DeterminantFeature::class -> object : DeterminantFeature { - override val determinant: Float by lazy { CommonOps_FDRM.det(origin) } - } - - SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature { - private val svd by lazy { - DecompositionFactory_FDRM.svd(origin.numRows, origin.numCols, true, true, false) - .apply { decompose(origin.copy()) } - } - - override val u: Matrix by lazy { svd.getU(null, false).wrapMatrix() } - override val s: Matrix by lazy { svd.getW(null).wrapMatrix() } - override val v: Matrix by lazy { svd.getV(null, false).wrapMatrix() } - override val singularValues: Point by lazy { FloatBuffer(svd.singularValues) } - } - - QRDecompositionFeature::class -> object : QRDecompositionFeature { - private val qr by lazy { - DecompositionFactory_FDRM.qr().apply { decompose(origin.copy()) } - } - - override val q: Matrix by lazy { - qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature) - } - - override val r: Matrix by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) } - } - - CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature { - override val l: Matrix by lazy { - val cholesky = - DecompositionFactory_FDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) } - - cholesky.getT(null).wrapMatrix().withFeature(LFeature) - } - } - - LupDecompositionFeature::class -> object : LupDecompositionFeature { - private val lup by lazy { - DecompositionFactory_FDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) } - } - - override val l: Matrix by lazy { - lup.getLower(null).wrapMatrix().withFeature(LFeature) - } - - override val u: Matrix by lazy { - lup.getUpper(null).wrapMatrix().withFeature(UFeature) - } - - override val p: Matrix by lazy { lup.getRowPivot(null).wrapMatrix() } - } - - else -> null - }?.let{ - type.cast(it) - } - } - - /** - * Solves for *x* in the following equation: *x = [a] -1 · [b]*. - * - * @param a the base matrix. - * @param b n by p matrix. - * @return the solution for *x* that is n by p. - */ - public fun solve(a: Matrix, b: Matrix): EjmlFloatMatrix { - val res = FMatrixRMaj(1, 1) - CommonOps_FDRM.solve(FMatrixRMaj(a.toEjml().origin), FMatrixRMaj(b.toEjml().origin), res) - return res.wrapMatrix() - } - - /** - * Solves for *x* in the following equation: *x = [a] -1 · [b]*. - * - * @param a the base matrix. - * @param b n by p vector. - * @return the solution for *x* that is n by p. - */ - public fun solve(a: Matrix, b: Point): EjmlFloatVector { - val res = FMatrixRMaj(1, 1) - CommonOps_FDRM.solve(FMatrixRMaj(a.toEjml().origin), FMatrixRMaj(b.toEjml().origin), res) - return EjmlFloatVector(res) - } -} - -/** - * [EjmlLinearSpace] implementation based on [CommonOps_DSCC], [DecompositionFactory_DSCC] operations and - * [DMatrixSparseCSC] matrices. - */ -public object EjmlLinearSpaceDSCC : EjmlLinearSpace() { - /** - * The [DoubleField] reference. - */ - override val elementAlgebra: DoubleField get() = DoubleField - - @Suppress("UNCHECKED_CAST") - override fun Matrix.toEjml(): EjmlDoubleMatrix = when { - this is EjmlDoubleMatrix<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleMatrix - else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) } - } - - @Suppress("UNCHECKED_CAST") - override fun Point.toEjml(): EjmlDoubleVector = when { - this is EjmlDoubleVector<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleVector - else -> EjmlDoubleVector(DMatrixSparseCSC(size, 1).also { - (0 until it.numRows).forEach { row -> it[row, 0] = get(row) } - }) - } - - override fun buildMatrix( - rows: Int, - columns: Int, - initializer: DoubleField.(i: Int, j: Int) -> Double, - ): EjmlDoubleMatrix = DMatrixSparseCSC(rows, columns).also { - (0 until rows).forEach { row -> - (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) } - } - }.wrapMatrix() - - override fun buildVector( - size: Int, - initializer: DoubleField.(Int) -> Double, - ): EjmlDoubleVector = EjmlDoubleVector(DMatrixSparseCSC(size, 1).also { - (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) } - }) - - private fun T.wrapMatrix() = EjmlDoubleMatrix(this) - private fun T.wrapVector() = EjmlDoubleVector(this) - - override fun Matrix.unaryMinus(): Matrix = this * elementAlgebra { -one } - - override fun Matrix.dot(other: Matrix): EjmlDoubleMatrix { - val out = DMatrixSparseCSC(1, 1) - CommonOps_DSCC.mult(toEjml().origin, other.toEjml().origin, out) - return out.wrapMatrix() - } - - override fun Matrix.dot(vector: Point): EjmlDoubleVector { - val out = DMatrixSparseCSC(1, 1) - CommonOps_DSCC.mult(toEjml().origin, vector.toEjml().origin, out) - return out.wrapVector() - } - - override operator fun Matrix.minus(other: Matrix): EjmlDoubleMatrix { - val out = DMatrixSparseCSC(1, 1) - - CommonOps_DSCC.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra { -one }, - other.toEjml().origin, - out, - null, - null, - ) - - return out.wrapMatrix() - } - - override operator fun Matrix.times(value: Double): EjmlDoubleMatrix { - val res = DMatrixSparseCSC(1, 1) - CommonOps_DSCC.scale(value, toEjml().origin, res) - return res.wrapMatrix() - } - - override fun Point.unaryMinus(): EjmlDoubleVector { - val res = DMatrixSparseCSC(1, 1) - CommonOps_DSCC.changeSign(toEjml().origin, res) - return res.wrapVector() - } - - override fun Matrix.plus(other: Matrix): EjmlDoubleMatrix { - val out = DMatrixSparseCSC(1, 1) - - CommonOps_DSCC.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra.one, - other.toEjml().origin, - out, - null, - null, - ) - - return out.wrapMatrix() - } - - override fun Point.plus(other: Point): EjmlDoubleVector { - val out = DMatrixSparseCSC(1, 1) - - CommonOps_DSCC.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra.one, - other.toEjml().origin, - out, - null, - null, - ) - - return out.wrapVector() - } - - override fun Point.minus(other: Point): EjmlDoubleVector { - val out = DMatrixSparseCSC(1, 1) - - CommonOps_DSCC.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra { -one }, - other.toEjml().origin, - out, - null, - null, - ) - - return out.wrapVector() - } - - override fun Double.times(m: Matrix): EjmlDoubleMatrix = m * this - - override fun Point.times(value: Double): EjmlDoubleVector { - val res = DMatrixSparseCSC(1, 1) - CommonOps_DSCC.scale(value, toEjml().origin, res) - return res.wrapVector() - } - - override fun Double.times(v: Point): EjmlDoubleVector = v * this - - @UnstableKMathAPI - override fun computeFeature(structure: Matrix, type: KClass): F? { - structure.getFeature(type)?.let { return it } - val origin = structure.toEjml().origin - - return when (type) { - QRDecompositionFeature::class -> object : QRDecompositionFeature { - private val qr by lazy { - DecompositionFactory_DSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) } - } - - override val q: Matrix by lazy { - qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature) - } - - override val r: Matrix by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) } - } - - CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature { - override val l: Matrix by lazy { - val cholesky = - DecompositionFactory_DSCC.cholesky().apply { decompose(origin.copy()) } - - (cholesky.getT(null) as DMatrix).wrapMatrix().withFeature(LFeature) - } - } - - LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object : - LUDecompositionFeature, DeterminantFeature, InverseMatrixFeature { - private val lu by lazy { - DecompositionFactory_DSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) } - } - - override val l: Matrix by lazy { - lu.getLower(null).wrapMatrix().withFeature(LFeature) - } - - override val u: Matrix by lazy { - lu.getUpper(null).wrapMatrix().withFeature(UFeature) - } - - override val inverse: Matrix by lazy { - var a = origin - val inverse = DMatrixRMaj(1, 1) - val solver = LinearSolverFactory_DSCC.lu(FillReducing.NONE) - if (solver.modifiesA()) a = a.copy() - val i = CommonOps_DDRM.identity(a.numRows) - solver.solve(i, inverse) - inverse.wrapMatrix() - } - - override val determinant: Double by lazy { elementAlgebra.number(lu.computeDeterminant().real) } - } - - else -> null - }?.let{ - type.cast(it) - } - } - - /** - * Solves for *x* in the following equation: *x = [a] -1 · [b]*. - * - * @param a the base matrix. - * @param b n by p matrix. - * @return the solution for *x* that is n by p. - */ - public fun solve(a: Matrix, b: Matrix): EjmlDoubleMatrix { - val res = DMatrixSparseCSC(1, 1) - CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res) - return res.wrapMatrix() - } - - /** - * Solves for *x* in the following equation: *x = [a] -1 · [b]*. - * - * @param a the base matrix. - * @param b n by p vector. - * @return the solution for *x* that is n by p. - */ - public fun solve(a: Matrix, b: Point): EjmlDoubleVector { - val res = DMatrixSparseCSC(1, 1) - CommonOps_DSCC.solve(DMatrixSparseCSC(a.toEjml().origin), DMatrixSparseCSC(b.toEjml().origin), res) - return EjmlDoubleVector(res) - } -} - -/** - * [EjmlLinearSpace] implementation based on [CommonOps_FSCC], [DecompositionFactory_FSCC] operations and - * [FMatrixSparseCSC] matrices. - */ -public object EjmlLinearSpaceFSCC : EjmlLinearSpace() { - /** - * The [FloatField] reference. - */ - override val elementAlgebra: FloatField get() = FloatField - - @Suppress("UNCHECKED_CAST") - override fun Matrix.toEjml(): EjmlFloatMatrix = when { - this is EjmlFloatMatrix<*> && origin is FMatrixSparseCSC -> this as EjmlFloatMatrix - else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) } - } - - @Suppress("UNCHECKED_CAST") - override fun Point.toEjml(): EjmlFloatVector = when { - this is EjmlFloatVector<*> && origin is FMatrixSparseCSC -> this as EjmlFloatVector - else -> EjmlFloatVector(FMatrixSparseCSC(size, 1).also { - (0 until it.numRows).forEach { row -> it[row, 0] = get(row) } - }) - } - - override fun buildMatrix( - rows: Int, - columns: Int, - initializer: FloatField.(i: Int, j: Int) -> Float, - ): EjmlFloatMatrix = FMatrixSparseCSC(rows, columns).also { - (0 until rows).forEach { row -> - (0 until columns).forEach { col -> it[row, col] = elementAlgebra.initializer(row, col) } - } - }.wrapMatrix() - - override fun buildVector( - size: Int, - initializer: FloatField.(Int) -> Float, - ): EjmlFloatVector = EjmlFloatVector(FMatrixSparseCSC(size, 1).also { - (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) } - }) - - private fun T.wrapMatrix() = EjmlFloatMatrix(this) - private fun T.wrapVector() = EjmlFloatVector(this) - - override fun Matrix.unaryMinus(): Matrix = this * elementAlgebra { -one } - - override fun Matrix.dot(other: Matrix): EjmlFloatMatrix { - val out = FMatrixSparseCSC(1, 1) - CommonOps_FSCC.mult(toEjml().origin, other.toEjml().origin, out) - return out.wrapMatrix() - } - - override fun Matrix.dot(vector: Point): EjmlFloatVector { - val out = FMatrixSparseCSC(1, 1) - CommonOps_FSCC.mult(toEjml().origin, vector.toEjml().origin, out) - return out.wrapVector() - } - - override operator fun Matrix.minus(other: Matrix): EjmlFloatMatrix { - val out = FMatrixSparseCSC(1, 1) - - CommonOps_FSCC.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra { -one }, - other.toEjml().origin, - out, - null, - null, - ) - - return out.wrapMatrix() - } - - override operator fun Matrix.times(value: Float): EjmlFloatMatrix { - val res = FMatrixSparseCSC(1, 1) - CommonOps_FSCC.scale(value, toEjml().origin, res) - return res.wrapMatrix() - } - - override fun Point.unaryMinus(): EjmlFloatVector { - val res = FMatrixSparseCSC(1, 1) - CommonOps_FSCC.changeSign(toEjml().origin, res) - return res.wrapVector() - } - - override fun Matrix.plus(other: Matrix): EjmlFloatMatrix { - val out = FMatrixSparseCSC(1, 1) - - CommonOps_FSCC.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra.one, - other.toEjml().origin, - out, - null, - null, - ) - - return out.wrapMatrix() - } - - override fun Point.plus(other: Point): EjmlFloatVector { - val out = FMatrixSparseCSC(1, 1) - - CommonOps_FSCC.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra.one, - other.toEjml().origin, - out, - null, - null, - ) - - return out.wrapVector() - } - - override fun Point.minus(other: Point): EjmlFloatVector { - val out = FMatrixSparseCSC(1, 1) - - CommonOps_FSCC.add( - elementAlgebra.one, - toEjml().origin, - elementAlgebra { -one }, - other.toEjml().origin, - out, - null, - null, - ) - - return out.wrapVector() - } - - override fun Float.times(m: Matrix): EjmlFloatMatrix = m * this - - override fun Point.times(value: Float): EjmlFloatVector { - val res = FMatrixSparseCSC(1, 1) - CommonOps_FSCC.scale(value, toEjml().origin, res) - return res.wrapVector() - } - - override fun Float.times(v: Point): EjmlFloatVector = v * this - - @UnstableKMathAPI - override fun computeFeature(structure: Matrix, type: KClass): F? { - structure.getFeature(type)?.let { return it } - val origin = structure.toEjml().origin - - return when (type) { - QRDecompositionFeature::class -> object : QRDecompositionFeature { - private val qr by lazy { - DecompositionFactory_FSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) } - } - - override val q: Matrix by lazy { - qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature) - } - - override val r: Matrix by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) } - } - - CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature { - override val l: Matrix by lazy { - val cholesky = - DecompositionFactory_FSCC.cholesky().apply { decompose(origin.copy()) } - - (cholesky.getT(null) as FMatrix).wrapMatrix().withFeature(LFeature) - } - } - - LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object : - LUDecompositionFeature, DeterminantFeature, InverseMatrixFeature { - private val lu by lazy { - DecompositionFactory_FSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) } - } - - override val l: Matrix by lazy { - lu.getLower(null).wrapMatrix().withFeature(LFeature) - } - - override val u: Matrix by lazy { - lu.getUpper(null).wrapMatrix().withFeature(UFeature) - } - - override val inverse: Matrix by lazy { - var a = origin - val inverse = FMatrixRMaj(1, 1) - val solver = LinearSolverFactory_FSCC.lu(FillReducing.NONE) - if (solver.modifiesA()) a = a.copy() - val i = CommonOps_FDRM.identity(a.numRows) - solver.solve(i, inverse) - inverse.wrapMatrix() - } - - override val determinant: Float by lazy { elementAlgebra.number(lu.computeDeterminant().real) } - } - - else -> null - }?.let{ - type.cast(it) - } - } - - /** - * Solves for *x* in the following equation: *x = [a] -1 · [b]*. - * - * @param a the base matrix. - * @param b n by p matrix. - * @return the solution for *x* that is n by p. - */ - public fun solve(a: Matrix, b: Matrix): EjmlFloatMatrix { - val res = FMatrixSparseCSC(1, 1) - CommonOps_FSCC.solve(FMatrixSparseCSC(a.toEjml().origin), FMatrixSparseCSC(b.toEjml().origin), res) - return res.wrapMatrix() - } - - /** - * Solves for *x* in the following equation: *x = [a] -1 · [b]*. - * - * @param a the base matrix. - * @param b n by p vector. - * @return the solution for *x* that is n by p. - */ - public fun solve(a: Matrix, b: Point): EjmlFloatVector { - val res = FMatrixSparseCSC(1, 1) - CommonOps_FSCC.solve(FMatrixSparseCSC(a.toEjml().origin), FMatrixSparseCSC(b.toEjml().origin), res) - return EjmlFloatVector(res) - } -} - diff --git a/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt b/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt index e89810e0d..04cc69708 100644 --- a/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt +++ b/kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt @@ -61,9 +61,9 @@ internal class EjmlMatrixTest { fun features() { val m = randomMatrix val w = EjmlDoubleMatrix(m) - val det: DeterminantFeature = EjmlLinearSpaceDDRM.computeFeature(w) ?: fail() + val det: Determinant = EjmlLinearSpaceDDRM.attributeFor(w) ?: fail() assertEquals(CommonOps_DDRM.det(m), det.determinant) - val lup: LupDecompositionFeature = EjmlLinearSpaceDDRM.computeFeature(w) ?: fail() + val lup: LupDecompositionAttribute = EjmlLinearSpaceDDRM.attributeFor(w) ?: fail() val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows, m.numCols) .also { it.decompose(m.copy()) } diff --git a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt index 411a35188..08dc75789 100644 --- a/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt +++ b/kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt @@ -16,7 +16,6 @@ import kotlin.math.pow public typealias DoubleVector = Point -@Suppress("FunctionName") public fun DoubleVector(vararg doubles: Double): DoubleVector = doubles.asBuffer() /** diff --git a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt index af84f47f2..b8a11c472 100644 --- a/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt +++ b/kmath-functions/src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt @@ -7,11 +7,14 @@ package space.kscience.kmath.functions +import space.kscience.kmath.UnstableKMathAPI import space.kscience.kmath.operations.Ring import space.kscience.kmath.operations.ScaleOperations import space.kscience.kmath.operations.invoke import kotlin.math.max import kotlin.math.min +import kotlin.reflect.KType +import kotlin.reflect.typeOf /** @@ -48,7 +51,7 @@ public data class Polynomial( * * @usesMathJax */ - public val coefficients: List + public val coefficients: List, ) { override fun toString(): String = "Polynomial$coefficients" } @@ -62,16 +65,17 @@ public data class Polynomial( * @param ring underlying ring of constants of type [A]. */ public open class PolynomialSpace( - /** - * Underlying ring of constants. Its operations on constants are used by local operations on constants and polynomials. - */ public val ring: A, ) : Ring>, ScaleOperations> where A : Ring, A : ScaleOperations { + @UnstableKMathAPI + override val elementType: KType get() = typeOf>() + /** * Instance of zero constant (zero of the underlying ring). */ public val constantZero: C get() = ring.zero + /** * Instance of unit constant (unit of the underlying ring). */ @@ -95,6 +99,7 @@ public open class PolynomialSpace( ) } } + /** * Returns difference between the constant represented as a polynomial and the polynomial. */ @@ -115,6 +120,7 @@ public open class PolynomialSpace( ) } } + /** * Returns product of the constant represented as a polynomial and the polynomial. */ @@ -147,6 +153,7 @@ public open class PolynomialSpace( ) } } + /** * Returns difference between the constant represented as a polynomial and the polynomial. */ @@ -165,6 +172,7 @@ public open class PolynomialSpace( ) } } + /** * Returns product of the constant represented as a polynomial and the polynomial. */ @@ -183,6 +191,7 @@ public open class PolynomialSpace( * Converts the constant [value] to polynomial. */ public fun number(value: C): Polynomial = Polynomial(listOf(value)) + /** * Converts the constant to polynomial. */ @@ -194,6 +203,7 @@ public open class PolynomialSpace( public override operator fun Polynomial.unaryMinus(): Polynomial = ring { Polynomial(coefficients.map { -it }) } + /** * Returns sum of the polynomials. */ @@ -210,6 +220,7 @@ public open class PolynomialSpace( } ) } + /** * Returns difference of the polynomials. */ @@ -226,6 +237,7 @@ public open class PolynomialSpace( } ) } + /** * Returns product of the polynomials. */ @@ -245,6 +257,7 @@ public open class PolynomialSpace( * Instance of zero polynomial (zero of the polynomial ring). */ override val zero: Polynomial = Polynomial(emptyList()) + /** * Instance of unit polynomial (unit of the polynomial ring). */ diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/NormalDistribution.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/NormalDistribution.kt index ae814254b..4cd1b89af 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/NormalDistribution.kt +++ b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/NormalDistribution.kt @@ -5,13 +5,12 @@ package space.kscience.kmath.distributions +import space.kscience.kmath.chains.BlockingDoubleChain import space.kscience.kmath.chains.Chain import space.kscience.kmath.operations.DoubleField.pow import space.kscience.kmath.random.RandomGenerator -import space.kscience.kmath.samplers.GaussianSampler +import space.kscience.kmath.samplers.* import space.kscience.kmath.samplers.InternalErf -import space.kscience.kmath.samplers.NormalizedGaussianSampler -import space.kscience.kmath.samplers.ZigguratNormalizedGaussianSampler import kotlin.math.* /** @@ -24,7 +23,7 @@ public class NormalDistribution(public val sampler: GaussianSampler) : Distribut return exp(-0.5 * x1 * x1 - (ln(sampler.standardDeviation) + 0.5 * ln(2 * PI))) } - override fun sample(generator: RandomGenerator): Chain = sampler.sample(generator) + override fun sample(generator: RandomGenerator): BlockingDoubleChain = sampler.sample(generator) override fun cumulative(arg: Double): Double { val dev = arg - sampler.mean