[WIP] Features to Attributes refactoring

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
 

attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attribute.kt

@@ -9,20 +9,31 @@ import kotlin.reflect.KType
 
 public interface Attribute<T>
 
+/**
+ * An attribute that could be either present or absent
+ */
+public interface FlagAttribute : Attribute<Unit>
 
+/**
+ * An attribute with a default value
+ */
 public interface AttributeWithDefault<T> : Attribute<T> {
     public val default: T
 }
 
+/**
+ * Attribute containing a set of values
+ */
 public interface SetAttribute<V> : Attribute<Set<V>>
 
 /**
  * An attribute that has a type parameter for value
+ * @param type parameter-type
  */
-public abstract class PolymorphicAttribute<T>(public val type: KType) : Attribute<T> {
+public abstract class PolymorphicAttribute<T>(public val type: SafeType<T>) : Attribute<T> {
-    override fun equals(other: Any?): Boolean = (other as? PolymorphicAttribute<*>)?.type == this.type
+    override fun equals(other: Any?): Boolean = other != null &&
+            (this::class == other::class) &&
+            (other as? PolymorphicAttribute<*>)?.type == this.type
 
-    override fun hashCode(): Int {
+    override fun hashCode(): Int = this::class.hashCode() + type.hashCode()
-        return type.hashCode()
-    }
 }

attributes-kt/src/commonMain/kotlin/space/kscience/attributes/Attributes.kt

@@ -24,18 +24,45 @@ public value class Attributes internal constructor(public val content: Map<out A
 
 public fun Attributes.isEmpty(): Boolean = content.isEmpty()
 
+/**
+ * Get attribute value or default
+ */
 public fun <T> Attributes.getOrDefault(attribute: AttributeWithDefault<T>): T = get(attribute) ?: attribute.default
 
-public fun <T, A : Attribute<T>> Attributes.withAttribute(
+/**
+ * Check if there is an attribute that matches given key by type and adheres to [predicate].
+ */
+@Suppress("UNCHECKED_CAST")
+public inline fun <T, reified A : Attribute<T>> 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 <T, reified A : Attribute<T>> 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 <reified A : FlagAttribute> Attributes.has(): Boolean =
+    content.keys.any { it is A }
+
+/**
+ * Create [Attributes] with an added or replaced attribute key.
+ */
+public fun <T : Any, A : Attribute<T>> Attributes.withAttribute(
     attribute: A,
-    attrValue: T?,
+    attrValue: T,
-): Attributes = Attributes(
+): Attributes = Attributes(content + (attribute to attrValue))
-    if (attrValue == null) {
+
-        content - attribute
+public fun <A : Attribute<Unit>> Attributes.withAttribute(attribute: A): Attributes =
-    } else {
+    withAttribute(attribute, Unit)
-        content + (attribute to attrValue)
+
-    }
+/**
-)
+ * 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 <T, A : SetAttribute<T>> Attributes.withoutAttributeElement(
     )
 }
 
+/**
+ * Create [Attributes] with a single key
+ */
 public fun <T : Any, A : Attribute<T>> Attributes(
     attribute: A,
     attrValue: T,

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<T> @PublishedApi internal constructor(public val kType: KType)
+
+public inline fun <reified T> safeTypeOf(): SafeType<T> = SafeType(typeOf<T>())
+
+/**
+ * Derive Kotlin [KClass] from this type and fail if the type is not a class (should not happen)
+ */
+@Suppress("UNCHECKED_CAST")
+@UnstableAttributesAPI
+public val <T> SafeType<T>.kClass: KClass<T & Any> get() = kType.classifier as KClass<T & Any>

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

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")

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}>

gradle.properties

@@ -12,4 +12,3 @@ org.gradle.parallel=true
 org.gradle.workers.max=4
 org.gradle.configureondemand=true
 org.gradle.jvmargs=-Xmx4096m
-

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

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<Double> {
     override val rowNum: Int get() = origin.rowDimension
@@ -38,6 +40,8 @@ public fun RealVector.toPoint(): CMVector = CMVector(this)
 public object CMLinearSpace : LinearSpace<Double, DoubleField> {
     override val elementAlgebra: DoubleField get() = DoubleField
 
+    override val elementType: KType = typeOf<Double>()
+
     override fun buildMatrix(
         rows: Int,
         columns: Int,
@@ -99,7 +103,7 @@ public object CMLinearSpace : LinearSpace<Double, DoubleField> {
         v * this
 
     @UnstableKMathAPI
-    override fun <F : StructureFeature> computeFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
+    override fun <F : StructureAttribute> computeFeature(structure: Matrix<Double>, type: KClass<out F>): 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<Double, DoubleField> {
         return when (type) {
             IsDiagonal::class -> if (origin is DiagonalMatrix) IsDiagonal else null
 
-            DeterminantFeature::class, LupDecompositionFeature::class -> object :
+            Determinant::class, LupDecompositionAttribute::class -> object :
-                DeterminantFeature<Double>,
+                Determinant<Double>,
-                LupDecompositionFeature<Double> {
+                LupDecompositionAttribute<Double> {
                 private val lup by lazy { LUDecomposition(origin) }
                 override val determinant: Double by lazy { lup.determinant }
-                override val l: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(lup.l).withFeature(LFeature) }
+                override val l: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(lup.l).withAttribute(LowerTriangular) }
-                override val u: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(lup.u).withFeature(UFeature) }
+                override val u: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(lup.u).withAttribute(UpperTriangular) }
                 override val p: Matrix<Double> by lazy { CMMatrix(lup.p) }
             }
 
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
+            CholeskyDecompositionAttribute::class -> object : CholeskyDecompositionAttribute<Double> {
                 override val l: Matrix<Double> by lazy<Matrix<Double>> {
                     val cholesky = CholeskyDecomposition(origin)
-                    CMMatrix(cholesky.l).withFeature(LFeature)
+                    CMMatrix(cholesky.l).withAttribute(LowerTriangular)
                 }
             }
 
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
+            QRDecompositionAttribute::class -> object : QRDecompositionAttribute<Double> {
                 private val qr by lazy { QRDecomposition(origin) }
-                override val q: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(qr.q).withFeature(OrthogonalFeature) }
+                override val q: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(qr.q).withAttribute(OrthogonalAttribute) }
-                override val r: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(qr.r).withFeature(UFeature) }
+                override val r: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(qr.r).withAttribute(UpperTriangular) }
             }
 
-            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Double> {
+            SingularValueDecompositionAttribute::class -> object : SingularValueDecompositionAttribute<Double> {
                 private val sv by lazy { SingularValueDecomposition(origin) }
                 override val u: Matrix<Double> by lazy { CMMatrix(sv.u) }
                 override val s: Matrix<Double> by lazy { CMMatrix(sv.s) }
                 override val v: Matrix<Double> by lazy { CMMatrix(sv.v) }
                 override val singularValues: Point<Double> by lazy { DoubleBuffer(sv.singularValues) }
             }
+
             else -> null
         }?.let(type::cast)
     }

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.

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<T, out A : Ring<T>>(

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<Double, DoubleField> {
 
@@ -20,7 +22,7 @@ public object DoubleLinearSpace : LinearSpace<Double, DoubleField> {
     override fun buildMatrix(
         rows: Int,
         columns: Int,
-        initializer: DoubleField.(i: Int, j: Int) -> Double
+        initializer: DoubleField.(i: Int, j: Int) -> Double,
     ): Matrix<Double> = DoubleFieldOpsND.structureND(ShapeND(rows, columns)) { (i, j) ->
         DoubleField.initializer(i, j)
     }.as2D()

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.*
-import space.kscience.kmath.nd.Structure2D
-import space.kscience.kmath.nd.StructureFeature
-import space.kscience.kmath.nd.as1D
 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<T> = MutableStructure2D<T>
  */
 public typealias Point<T> = Buffer<T>
 
+/**
+ * A marker interface for algebras that operate on matrices
+ * @param T type of matrix element
+ */
+public interface MatrixOperations<T>
+
 /**
  * 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<T, out A : Ring<T>> {
+public interface LinearSpace<T, out A : Ring<T>> : MatrixOperations<T> {
     public val elementAlgebra: A
 
     /**
@@ -167,16 +172,16 @@ public interface LinearSpace<T, out A : Ring<T>> {
     public operator fun T.times(v: Point<T>): Point<T> = 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 <F : StructureFeature> computeFeature(structure: Matrix<T>, type: KClass<out F>): F? =
+    public fun <T, A : StructureAttribute<T>> attributeFor(structure: StructureND<*>, attribute: A): T? =
-        structure.getFeature(type)
+        structure.attributes[attribute]
 
     public companion object {
 
@@ -184,23 +189,12 @@ public interface LinearSpace<T, out A : Ring<T>> {
          * A structured matrix with custom buffer
          */
         public fun <T : Any, A : Ring<T>> buffered(
-            algebra: A
+            algebra: A,
         ): LinearSpace<T, A> = 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 <T : Any, reified F : StructureFeature> LinearSpace<T, *>.computeFeature(structure: Matrix<T>): F? =
-    computeFeature(structure, F::class)
-
 
 public inline operator fun <LS : LinearSpace<*, *>, R> LS.invoke(block: LS.() -> R): R = run(block)
 

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.*
-import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.MutableBuffer
-import space.kscience.kmath.structures.MutableBufferFactory
 
 /**
- * Common implementation of [LupDecompositionFeature].
+ * Matrices with this feature support LU factorization with partial pivoting: *[p] · a = [l] · [u]* where
- */
+ * *a* is the owning matrix.
-public class LupDecomposition<T : Any>(
-    public val context: LinearSpace<T, *>,
-    public val elementContext: Field<T>,
-    public val lu: Matrix<T>,
-    public val pivot: IntArray,
-    private val even: Boolean,
-) : LupDecompositionFeature<T>, DeterminantFeature<T> {
-    /**
-     * Returns the matrix L of the decomposition.
  *
-     * L is a lower-triangular matrix with [Ring.one] in diagonal
+ * @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
  */
-    override val l: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j ->
+public data class LupDecomposition<T>(
-        when {
+    public val l: Matrix<T>,
-            j < i -> lu[i, j]
+    public val u: Matrix<T>,
-            j == i -> elementContext.one
+    public val p: Matrix<T>,
-            else -> elementContext.zero
+)
-        }
-    }.withFeature(LFeature)
 
 
-    /**
+public class LupDecompositionAttribute<T>(type: SafeType<LupDecomposition<T>>) :
-     * Returns the matrix U of the decomposition.
+    PolymorphicAttribute<LupDecomposition<T>>(type),
-     *
+    MatrixAttribute<LupDecomposition<T>>
-     * U is an upper-triangular matrix including the diagonal
-     */
-    override val u: Matrix<T> = VirtualMatrix(lu.shape[0], lu.shape[1]) { i, j ->
-        if (j >= i) lu[i, j] else elementContext.zero
-    }.withFeature(UFeature)
 
-    /**
+public val <T> MatrixOperations<T>.LUP: LupDecompositionAttribute<T>
-     * Returns the P rows permutation matrix.
+    get() = LupDecompositionAttribute(safeTypeOf())
-     *
-     * 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<T> = 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] } }
-    }
 
-}
+///**
+// * Common implementation of [LupDecomposition].
+// */
+//private class LupDecompositionImpl<T : Any>(
+//    public val elementContext: Field<T>,
+//    public val lu: Matrix<T>,
+//    public val pivot: IntBuffer,
+//    private val even: Boolean,
+//) : LupDecomposition<T> {
+//    /**
+//     * Returns the matrix L of the decomposition.
+//     *
+//     * L is a lower-triangular matrix with [Ring.one] in diagonal
+//     */
+//    override val l: Matrix<T> = 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<T> = 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<T> = 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 <T : Comparable<T>> LinearSpace<T, Ring<T>>.abs(value: T): T =
@@ -73,7 +98,6 @@ internal fun <T : Comparable<T>> LinearSpace<T, Ring<T>>.abs(value: T): T =
  * Create a lup decomposition of generic matrix.
  */
 public fun <T : Comparable<T>> LinearSpace<T, Field<T>>.lup(
-    factory: MutableBufferFactory<T>,
     matrix: Matrix<T>,
     checkSingular: (T) -> Boolean,
 ): LupDecomposition<T> {
@@ -82,15 +106,15 @@ public fun <T : Comparable<T>> LinearSpace<T, Field<T>>.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 <T : Comparable<T>> LinearSpace<T, Field<T>>.lup(
                 for (row in col + 1 until m) lu[row, col] /= luDiag
             }
 
-            return LupDecomposition(this@lup, elementAlgebra, lu.collect(), pivot, even)
+            val l: MatrixWrapper<T> = 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 <reified T : Comparable<T>> LinearSpace<T, Field<T>>.lup(
-    matrix: Matrix<T>,
-    noinline checkSingular: (T) -> Boolean,
-): LupDecomposition<T> = lup(MutableBuffer.Companion::auto, matrix, checkSingular)
 
 public fun LinearSpace<Double, DoubleField>.lup(
     matrix: Matrix<Double>,
     singularityThreshold: Double = 1e-11,
-): LupDecomposition<Double> =
+): LupDecomposition<Double> = lup(matrix) { it < singularityThreshold }
-    lup(::DoubleBuffer, matrix) { it < singularityThreshold }
 
-internal fun <T : Any> LupDecomposition<T>.solve(
+internal fun <T : Any, A : Field<T>> LinearSpace<T, A>.solve(
-    factory: MutableBufferFactory<T>,
+    lup: LupDecomposition<T>,
     matrix: Matrix<T>,
 ): Matrix<T> {
-    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 {
+    BufferAccessor2D(matrix.rowNum, matrix.colNum, elementAlgebra.bufferFactory).run {
-        elementContext {
+        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 <T : Any> LupDecomposition<T>.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 <T : Any> LupDecomposition<T>.solve(
  */
 @OptIn(UnstableKMathAPI::class)
 public fun <T : Comparable<T>, F : Field<T>> LinearSpace<T, F>.lupSolver(
-    bufferFactory: MutableBufferFactory<T>,
     singularityCheck: (T) -> Boolean,
 ): LinearSolver<T> = object : LinearSolver<T> {
     override fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
         // Use existing decomposition if it is provided by matrix
-        val decomposition = computeFeature(a) ?: lup(bufferFactory, a, singularityCheck)
+        val decomposition = attributeFor(a, LUP) ?: lup(a, singularityCheck)
-        return decomposition.solve(bufferFactory, b)
+        return solve(decomposition, b)
     }
 
     override fun inverse(matrix: Matrix<T>): Matrix<T> = solve(matrix, one(matrix.rowNum, matrix.colNum))
 }
 
 public fun LinearSpace<Double, DoubleField>.lupSolver(singularityThreshold: Double = 1e-11): LinearSolver<Double> =
-    lupSolver(::DoubleBuffer) { it < singularityThreshold }
+    lupSolver { it < singularityThreshold }

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 <T : Any> LinearSpace<T, Ring<T>>.column(
 
 public fun <T : Any> LinearSpace<T, Ring<T>>.column(vararg values: T): Matrix<T> = column(values.size, values::get)
 
-public object SymmetricMatrixFeature : MatrixFeature
+public object Symmetric : MatrixAttribute<Unit>, 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 <T : Any, A : Ring<T>> MatrixBuilder<T, A>.symmetric(
             } else {
                 cached
             }
-        }.withFeature(SymmetricMatrixFeature)
+        }.withAttribute(Symmetric)
     }
 }

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<T> : Attribute<T>
+public interface MatrixAttribute<T> : StructureAttribute<T>
 
 /**
  * Matrices with this feature are considered to have only diagonal non-zero elements.
  */
-public interface IsDiagonal : MatrixFeature<Unit> {
+public interface IsDiagonal : MatrixAttribute<Unit>, 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<sup>&minus;1</sup>* where *a* is the owning matrix.
+ * Matrices with this feature can be inverted.
  *
  * @param T the type of matrices' items.
  */
-public class Inverted<T> private constructor() : MatrixFeature<Matrix<T>> {
+public class Inverted<T>(type: SafeType<Matrix<T>>) :
-    internal val instance: Inverted<Nothing> = Inverted()
+    PolymorphicAttribute<Matrix<T>>(type),
-}
+    MatrixAttribute<Matrix<T>>
 
-@Suppress("UNCHECKED_CAST")
+public val <T> MatrixOperations<T>.Inverted: Inverted<T> get() = Inverted(safeTypeOf())
-public val <T> LinearSpace<T, *>.Inverted: Inverted<T> get() = Inverted.instance as Inverted<T>
 
 /**
  * Matrices with this feature can compute their determinant.
  *
  * @param T the type of matrices' items.
  */
-public class DeterminantFeature<T : Any> : MatrixFeature<T>
+public class Determinant<T>(type: SafeType<T>) :
+    PolymorphicAttribute<T>(type),
+    MatrixAttribute<T>
 
-/**
+public inline val <reified T> MatrixOperations<T>.Determinant: Determinant<T> get() = Determinant(safeTypeOf())
- * Produces a [DeterminantFeature] where the [DeterminantFeature.determinant] is [determinant].
- *
- * @param determinant the value of determinant.
- * @return a new [DeterminantFeature].
- */
-@Suppress("FunctionName")
-public fun <T : Any> DeterminantFeature(determinant: T): DeterminantFeature<T> = object : DeterminantFeature<T> {
-    override val determinant: T = determinant
-}
 
 /**
  * Matrices with this feature are lower triangular ones.
  */
-public object LFeature : MatrixFeature<Unit>
+public object LowerTriangular : MatrixAttribute<Unit>, FlagAttribute
 
 /**
  * Matrices with this feature are upper triangular ones.
  */
-public object UFeature : MatrixFeature<Unit>
+public object UpperTriangular : MatrixAttribute<Unit>, FlagAttribute
+
+/**
+ * Matrices with this feature support LU factorization: *a = [l] &middot; [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<T>(val l: Matrix<T>, val u: Matrix<T>)
 
 /**
  * Matrices with this feature support LU factorization: *a = [l] &middot; [u]* where *a* is the owning matrix.
  *
  * @param T the type of matrices' items.
  */
-public interface LUDecompositionFeature<out T : Any> : MatrixFeature {
+public class LuDecompositionAttribute<T>(type: SafeType<LUDecomposition<T>>) :
-    /**
+    PolymorphicAttribute<LUDecomposition<T>>(type),
-     * The lower triangular matrix in this decomposition. It may have [LFeature].
+    MatrixAttribute<LUDecomposition<T>>
-     */
-    public val l: Matrix<T>
 
-    /**
+public val <T> MatrixOperations<T>.LU: LuDecompositionAttribute<T> get() = LuDecompositionAttribute(safeTypeOf())
-     * The upper triangular matrix in this decomposition. It may have [UFeature].
-     */
-    public val u: Matrix<T>
-}
 
-/**
- * Matrices with this feature support LU factorization with partial pivoting: *[p] &middot; a = [l] &middot; [u]* where
- * *a* is the owning matrix.
- *
- * @param T the type of matrices' items.
- */
-public interface LupDecompositionFeature<out T : Any> : MatrixFeature {
-    /**
-     * The lower triangular matrix in this decomposition. It may have [LFeature].
-     */
-    public val l: Matrix<T>
-
-    /**
-     * The upper triangular matrix in this decomposition. It may have [UFeature].
-     */
-    public val u: Matrix<T>
-
-    /**
-     * The permutation matrix in this decomposition.
-     */
-    public val p: Matrix<T>
-}
 
 /**
  * Matrices with this feature are orthogonal ones: *a &middot; a<sup>T</sup> = u* where *a* is the owning matrix, *u*
  * is the unit matrix ([IsUnit]).
  */
-public object OrthogonalFeature : MatrixFeature
+public object OrthogonalAttribute : MatrixAttribute<Unit>, FlagAttribute
 
+
+public interface QRDecomposition<out T> {
     /**
- * Matrices with this feature support QR factorization: *a = [q] &middot; [r]* where *a* is the owning matrix.
+     * The orthogonal matrix in this decomposition. It may have [OrthogonalAttribute].
- *
- * @param T the type of matrices' items.
- */
-public interface QRDecompositionFeature<out T : Any> : MatrixFeature {
-    /**
-     * The orthogonal matrix in this decomposition. It may have [OrthogonalFeature].
      */
     public val q: Matrix<T>
 
     /**
-     * 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<T>
 }
 
+/**
+ * Matrices with this feature support QR factorization: *a = [QR.q] &middot; [QR.r]* where *a* is the owning matrix.
+ *
+ * @param T the type of matrices' items.
+ */
+public class QRDecompositionAttribute<T>(type: SafeType<QRDecomposition<T>>) :
+    PolymorphicAttribute<QRDecomposition<T>>(type),
+    MatrixAttribute<QRDecomposition<T>>
+
+public val <T> MatrixOperations<T>.QR: QRDecompositionAttribute<T>
+    get() = QRDecompositionAttribute(safeTypeOf())
+
+public interface CholeskyDecomposition<T> {
+    /**
+     * The triangular matrix in this decomposition. It may have either [UpperTriangular] or [LowerTriangular].
+     */
+    public val l: Matrix<T>
+}
+
 /**
  * Matrices with this feature support Cholesky factorization: *a = [l] &middot; [l]<sup>H</sup>* where *a* is the
  * owning matrix.
  *
  * @param T the type of matrices' items.
  */
-public interface CholeskyDecompositionFeature<out T : Any> : MatrixFeature {
+public class CholeskyDecompositionAttribute<T>(type: SafeType<CholeskyDecomposition<T>>) :
-    /**
+    PolymorphicAttribute<CholeskyDecomposition<T>>(type),
-     * The triangular matrix in this decomposition. It may have either [UFeature] or [LFeature].
+    MatrixAttribute<CholeskyDecomposition<T>>
-     */
-    public val l: Matrix<T>
-}
 
+public val <T> MatrixOperations<T>.Cholesky: CholeskyDecompositionAttribute<T>
+    get() = CholeskyDecompositionAttribute(safeTypeOf())
+
+public interface SingularValueDecomposition<T> {
     /**
- * Matrices with this feature support SVD: *a = [u] &middot; [s] &middot; [v]<sup>H</sup>* where *a* is the owning
+     * The matrix in this decomposition. It is unitary, and it consists of left singular vectors.
- * matrix.
- *
- * @param T the type of matrices' items.
- */
-public interface SingularValueDecompositionFeature<out T : Any> : MatrixFeature {
-    /**
-     * The matrix in this decomposition. It is unitary, and it consists from left singular vectors.
      */
     public val u: Matrix<T>
 
@@ -164,14 +151,27 @@ public interface SingularValueDecompositionFeature<out T : Any> : MatrixFeature
     public val s: Matrix<T>
 
     /**
-     * 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<T>
 
     /**
-     * The buffer of singular values of this SVD.
+     * The buffer of singular values for this SVD.
      */
     public val singularValues: Point<T>
 }
 
+/**
+ * Matrices with this feature support SVD: *a = [u] &middot; [s] &middot; [v]<sup>H</sup>* where *a* is the owning
+ * matrix.
+ *
+ * @param T the type of matrices' items.
+ */
+public class SingularValueDecompositionAttribute<T>(type: SafeType<SingularValueDecomposition<T>>) :
+    PolymorphicAttribute<SingularValueDecomposition<T>>(type),
+    MatrixAttribute<SingularValueDecomposition<T>>
+
+public val <T> MatrixOperations<T>.SVD: SingularValueDecompositionAttribute<T>
+    get() = SingularValueDecompositionAttribute(safeTypeOf())
+
 //TODO add sparse matrix feature

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<out T : Any> internal constructor(
     public val origin: Matrix<T>,
-    public val attributes: Attributes,
+    override val attributes: Attributes,
 ) : Matrix<T> by origin {
 
     override fun toString(): String = "MatrixWrapper(matrix=$origin, features=$attributes)"
@@ -34,23 +35,31 @@ public val <T : Any> Matrix<T>.origin: Matrix<T>
 /**
  * Add a single feature to a [Matrix]
  */
-public fun <T : Any> Matrix<T>.withFeature(newFeature: MatrixFeature): MatrixWrapper<T> = if (this is MatrixWrapper) {
+public fun <T : Any, A : Attribute<T>> Matrix<T>.withAttribute(
-    MatrixWrapper(origin, attributes.with(newFeature))
+    attribute: A,
+    attrValue: T,
+): MatrixWrapper<T> = 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 fun <T : Any, A : Attribute<Unit>> Matrix<T>.withAttribute(
-public operator fun <T : Any> Matrix<T>.plus(newFeature: MatrixFeature): MatrixWrapper<T> = withFeature(newFeature)
+    attribute: A,
+): MatrixWrapper<T> = 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 <T : Any> Matrix<T>.withFeatures(newFeatures: Iterable<MatrixFeature>): MatrixWrapper<T> =
+public fun <T : Any> Matrix<T>.modifyAttributes(modifier: (Attributes) -> Attributes): MatrixWrapper<T> =
     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 <T : Any> Matrix<T>.withFeatures(newFeatures: Iterable<MatrixFeature>
 public fun <T : Any> LinearSpace<T, Ring<T>>.one(
     rows: Int,
     columns: Int,
-): Matrix<T> = VirtualMatrix(rows, columns) { i, j ->
+): MatrixWrapper<T> = VirtualMatrix(rows, columns) { i, j ->
     if (i == j) elementAlgebra.one else elementAlgebra.zero
-}.withFeature(IsUnit)
+}.withAttribute(IsUnit)
 
 
 /**
@@ -70,16 +79,16 @@ public fun <T : Any> LinearSpace<T, Ring<T>>.one(
 public fun <T : Any> LinearSpace<T, Ring<T>>.zero(
     rows: Int,
     columns: Int,
-): Matrix<T> = VirtualMatrix(rows, columns) { _, _ ->
+): MatrixWrapper<T> = VirtualMatrix(rows, columns) { _, _ ->
     elementAlgebra.zero
-}.withFeature(IsZero)
+}.withAttribute(IsZero)
 
-public class TransposedFeature<out T : Any>(public val original: Matrix<T>) : MatrixFeature
+public class TransposedAttribute<out T : Any>(public val original: Matrix<T>) : MatrixAttribute
 
 /**
  * Create a virtual transposed matrix without copying anything. `A.transpose().transpose() === A`
  */
 @Suppress("UNCHECKED_CAST")
 @OptIn(UnstableKMathAPI::class)
-public fun <T : Any> Matrix<T>.transpose(): Matrix<T> = getFeature(TransposedFeature::class)?.original as? Matrix<T>
+public fun <T : Any> Matrix<T>.transpose(): Matrix<T> = getFeature(TransposedAttribute::class)?.original as? Matrix<T>
-    ?: VirtualMatrix(colNum, rowNum) { i, j -> get(j, i) }.withFeature(TransposedFeature(this))
+    ?: VirtualMatrix(colNum, rowNum) { i, j -> get(j, i) }.withAttribute(TransposedAttribute(this))

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<out T : Any>(
     override val rowNum: Int,
     override val colNum: Int,
+    override val attributes: Attributes = Attributes.EMPTY,
     public val generator: (i: Int, j: Int) -> T,
 ) : Matrix<T> {
 
@@ -24,5 +26,8 @@ public class VirtualMatrix<out T : Any>(
     override operator fun get(i: Int, j: Int): T = generator(i, j)
 }
 
-public fun <T : Any> MatrixBuilder<T, *>.virtual(generator: (i: Int, j: Int) -> T): VirtualMatrix<T> =
+public fun <T : Any> MatrixBuilder<T, *>.virtual(
-    VirtualMatrix(rows, columns, generator)
+    attributes: Attributes = Attributes.EMPTY,
+    generator: (i: Int, j: Int) -> T,
+): VirtualMatrix<T> =
+    VirtualMatrix(rows, columns, attributes, generator)

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt

@@ -79,7 +79,7 @@ public interface AlgebraND<T, out C : Algebra<T>> : Algebra<StructureND<T>> {
      * @return a feature object or `null` if it isn't present.
      */
     @UnstableKMathAPI
-    public fun <F : StructureFeature> getFeature(structure: StructureND<T>, type: KClass<out F>): F? =
+    public fun <F : StructureAttribute> getFeature(structure: StructureND<T>, type: KClass<out F>): F? =
         structure.getFeature(type)
 
     public companion object
@@ -93,7 +93,7 @@ public interface AlgebraND<T, out C : Algebra<T>> : Algebra<StructureND<T>> {
  * @return a feature object or `null` if it isn't present.
  */
 @UnstableKMathAPI
-public inline fun <T : Any, reified F : StructureFeature> AlgebraND<T, *>.getFeature(structure: StructureND<T>): F? =
+public inline fun <T : Any, reified F : StructureAttribute> AlgebraND<T, *>.getFeature(structure: StructureND<T>): F? =
     getFeature(structure, F::class)
 
 /**

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt

@@ -34,9 +34,9 @@ public open class BufferND<out T>(
 /**
  * Create a generic [BufferND] using provided [initializer]
  */
-public fun <T> BufferND(
+public inline fun <reified T> BufferND(
     shape: ShapeND,
-    bufferFactory: BufferFactory<T> = BufferFactory.boxing(),
+    bufferFactory: BufferFactory<T> = BufferFactory.auto(),
     initializer: (IntArray) -> T,
 ): BufferND<T> {
     val strides = Strides(shape)

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt

@@ -110,7 +110,7 @@ private value class Structure2DWrapper<out T>(val structure: StructureND<T>) : S
     @PerformancePitfall
     override operator fun get(i: Int, j: Int): T = structure[i, j]
 
-    override fun <F : StructureFeature> getFeature(type: KClass<out F>): F? = structure.getFeature(type)
+    override fun <F : StructureAttribute> getFeature(type: KClass<out F>): F? = structure.getFeature(type)
 
     @PerformancePitfall
     override fun elements(): Sequence<Pair<IntArray, T>> = structure.elements()

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<T> : Attribute<T>
+public interface StructureAttribute<T> : Attribute<T>
 
 /**
  * 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<out T> : AttributeContainer, WithShape {
          */
         public fun <T> buffered(
             strides: Strides,
-            bufferFactory: BufferFactory<T> = BufferFactory.boxing(),
             initializer: (IntArray) -> T,
-        ): BufferND<T> = BufferND(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
+        ): BufferND<T> = BufferND(strides, Buffer.boxing(strides.linearSize) { i -> initializer(strides.index(i)) })
+
+
+        public fun <T> buffered(
+            shape: ShapeND,
+            initializer: (IntArray) -> T,
+        ): BufferND<T> = buffered(ColumnStrides(shape), initializer)
 
         /**
          * Inline create NDStructure with non-boxing buffer implementation if it is possible
@@ -135,16 +140,11 @@ public interface StructureND<out T> : AttributeContainer, WithShape {
         ): BufferND<T> = BufferND(strides, Buffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) })
 
         public inline fun <T : Any> auto(
-            type: KClass<T>,
+            type: SafeType<T>,
             strides: Strides,
             crossinline initializer: (IntArray) -> T,
         ): BufferND<T> = BufferND(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) })
 
-        public fun <T> buffered(
-            shape: ShapeND,
-            bufferFactory: BufferFactory<T> = BufferFactory.boxing(),
-            initializer: (IntArray) -> T,
-        ): BufferND<T> = buffered(ColumnStrides(shape), bufferFactory, initializer)
 
         public inline fun <reified T : Any> auto(
             shape: ShapeND,
@@ -159,7 +159,7 @@ public interface StructureND<out T> : AttributeContainer, WithShape {
             auto(ColumnStrides(ShapeND(shape)), initializer)
 
         public inline fun <T : Any> auto(
-            type: KClass<T>,
+            type: SafeType<T>,
             vararg shape: Int,
             crossinline initializer: (IntArray) -> T,
         ): BufferND<T> = auto(type, ColumnStrides(ShapeND(shape)), initializer)

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<T> {
-
     /**
      * Provide a factory for buffers, associated with this [Algebra]
      */
@@ -67,12 +67,12 @@ public interface Algebra<T> {
      *
      * @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<T> {
         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:
      *

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<T, out A : Algebra<T>> : Algebra<Buffer<T>> {
     public val elementAlgebra: A
+
+    public val elementType: KType
+
     public val elementBufferFactory: MutableBufferFactory<T> get() = elementAlgebra.bufferFactory
 
     public fun buffer(size: Int, vararg elements: T): Buffer<T> {

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<out T> : WithSize {
          * The [size] is specified, and each element is calculated by calling the specified [initializer] function.
          */
         @Suppress("UNCHECKED_CAST")
-        public inline fun <T : Any> auto(type: KClass<T>, size: Int, initializer: (Int) -> T): Buffer<T> =
+        public inline fun <T> auto(type: SafeType<T>, size: Int, initializer: (Int) -> T): Buffer<T> =
-            when (type) {
+            when (type.kType) {
-                Double::class -> MutableBuffer.double(size) { initializer(it) as Double } as Buffer<T>
+                typeOf<Double>() -> MutableBuffer.double(size) { initializer(it) as Double } as Buffer<T>
-                Short::class -> MutableBuffer.short(size) { initializer(it) as Short } as Buffer<T>
+                typeOf<Short>() -> MutableBuffer.short(size) { initializer(it) as Short } as Buffer<T>
-                Int::class -> MutableBuffer.int(size) { initializer(it) as Int } as Buffer<T>
+                typeOf<Int>() -> MutableBuffer.int(size) { initializer(it) as Int } as Buffer<T>
-                Long::class -> MutableBuffer.long(size) { initializer(it) as Long } as Buffer<T>
+                typeOf<Long>() -> MutableBuffer.long(size) { initializer(it) as Long } as Buffer<T>
-                Float::class -> MutableBuffer.float(size) { initializer(it) as Float } as Buffer<T>
+                typeOf<Float>() -> MutableBuffer.float(size) { initializer(it) as Float } as Buffer<T>
                 else -> boxing(size, initializer)
             }
 
@@ -115,8 +117,8 @@ public interface Buffer<out T> : WithSize {
          *
          * The [size] is specified, and each element is calculated by calling the specified [initializer] function.
          */
-        public inline fun <reified T : Any> auto(size: Int, initializer: (Int) -> T): Buffer<T> =
+        public inline fun <reified T> auto(size: Int, initializer: (Int) -> T): Buffer<T> =
-            auto(T::class, size, initializer)
+            auto(safeTypeOf<T>(), size, initializer)
     }
 }
 

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/BufferAccessor2D.kt

@@ -27,12 +27,9 @@ internal class BufferAccessor2D<T>(
     fun create(mat: Structure2D<T>): MutableBuffer<T> = create { i, j -> mat[i, j] }
 
     //TODO optimize wrapper
-    fun MutableBuffer<T>.collect(): Structure2D<T> = StructureND.buffered(
+    fun MutableBuffer<T>.toStructure2D(): Structure2D<T> = StructureND.buffered(
-        ColumnStrides(ShapeND(rowNum, colNum)),
+        ColumnStrides(ShapeND(rowNum, colNum))
-        factory
+    ) { (i, j) -> get(i, j) }.as2D()
-    ) { (i, j) ->
-        get(i, j)
-    }.as2D()
 
     inner class Row(val buffer: MutableBuffer<T>, val rowIndex: Int) : MutableBuffer<T> {
         override val size: Int get() = colNum

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<T> : Buffer<T> {
          * The [size] is specified, and each element is calculated by calling the specified [initializer] function.
          */
         @Suppress("UNCHECKED_CAST")
-        public inline fun <T : Any> auto(type: KClass<out T>, size: Int, initializer: (Int) -> T): MutableBuffer<T> =
+        public inline fun <T> auto(type: SafeType<T>, size: Int, initializer: (Int) -> T): MutableBuffer<T> =
-            when (type) {
+            when (type.kType) {
-                Double::class -> double(size) { initializer(it) as Double } as MutableBuffer<T>
+                typeOf<Double>() -> double(size) { initializer(it) as Double } as MutableBuffer<T>
-                Short::class -> short(size) { initializer(it) as Short } as MutableBuffer<T>
+                typeOf<Short>() -> short(size) { initializer(it) as Short } as MutableBuffer<T>
-                Int::class -> int(size) { initializer(it) as Int } as MutableBuffer<T>
+                typeOf<Int>() -> int(size) { initializer(it) as Int } as MutableBuffer<T>
-                Float::class -> float(size) { initializer(it) as Float } as MutableBuffer<T>
+                typeOf<Float>() -> float(size) { initializer(it) as Float } as MutableBuffer<T>
-                Long::class -> long(size) { initializer(it) as Long } as MutableBuffer<T>
+                typeOf<Long>() -> long(size) { initializer(it) as Long } as MutableBuffer<T>
                 else -> boxing(size, initializer)
             }
 
@@ -90,9 +92,8 @@ public interface MutableBuffer<T> : Buffer<T> {
          *
          * The [size] is specified, and each element is calculated by calling the specified [initializer] function.
          */
-        @Suppress("UNCHECKED_CAST")
+        public inline fun <reified T> auto(size: Int, initializer: (Int) -> T): MutableBuffer<T> =
-        public inline fun <reified T : Any> auto(size: Int, initializer: (Int) -> T): MutableBuffer<T> =
+            auto(safeTypeOf<T>(), size, initializer)
-            auto(T::class, size, initializer)
     }
 }
 

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<T : Any, out A : Ring<T>, out M : org.ejml
     @Suppress("UNCHECKED_CAST")
     @UnstableKMathAPI
     public fun EjmlMatrix<T, *>.inverse(): Structure2D<Double> =
-        computeFeature(this, InverseMatrixFeature::class)?.inverse as Structure2D<Double>
+        attributeFor(this, InverseMatrixFeature::class)?.inverse as Structure2D<Double>
 }

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<Double> = EjmlLinearSpaceDDRM.computeFeature(w) ?: fail()
+        val det: Determinant<Double> = EjmlLinearSpaceDDRM.attributeFor(w) ?: fail()
         assertEquals(CommonOps_DDRM.det(m), det.determinant)
-        val lup: LupDecompositionFeature<Double> = EjmlLinearSpaceDDRM.computeFeature(w) ?: fail()
+        val lup: LupDecompositionAttribute<Double> = EjmlLinearSpaceDDRM.attributeFor(w) ?: fail()
 
         val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows, m.numCols)
             .also { it.decompose(m.copy()) }

kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt

@@ -16,7 +16,6 @@ import kotlin.math.pow
 
 public typealias DoubleVector = Point<Double>
 
-@Suppress("FunctionName")
 public fun DoubleVector(vararg doubles: Double): DoubleVector = doubles.asBuffer()
 
 /**

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<out C>(
      *
      * @usesMathJax
      */
-    public val coefficients: List<C>
+    public val coefficients: List<C>,
 ) {
     override fun toString(): String = "Polynomial$coefficients"
 }
@@ -62,16 +65,17 @@ public data class Polynomial<out C>(
  * @param ring underlying ring of constants of type [A].
  */
 public open class PolynomialSpace<C, A>(
-    /**
-     * Underlying ring of constants. Its operations on constants are used by local operations on constants and polynomials.
-     */
     public val ring: A,
 ) : Ring<Polynomial<C>>, ScaleOperations<Polynomial<C>> where A : Ring<C>, A : ScaleOperations<C> {
 
+    @UnstableKMathAPI
+    override val elementType: KType get() = typeOf<Polynomial<C>>()
+
     /**
      * 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<C, A>(
                 )
             }
         }
+
     /**
      * Returns difference between the constant represented as a polynomial and the polynomial.
      */
@@ -115,6 +120,7 @@ public open class PolynomialSpace<C, A>(
                 )
             }
         }
+
     /**
      * Returns product of the constant represented as a polynomial and the polynomial.
      */
@@ -147,6 +153,7 @@ public open class PolynomialSpace<C, A>(
                 )
             }
         }
+
     /**
      * Returns difference between the constant represented as a polynomial and the polynomial.
      */
@@ -165,6 +172,7 @@ public open class PolynomialSpace<C, A>(
                 )
             }
         }
+
     /**
      * Returns product of the constant represented as a polynomial and the polynomial.
      */
@@ -183,6 +191,7 @@ public open class PolynomialSpace<C, A>(
      * Converts the constant [value] to polynomial.
      */
     public fun number(value: C): Polynomial<C> = Polynomial(listOf(value))
+
     /**
      * Converts the constant to polynomial.
      */
@@ -194,6 +203,7 @@ public open class PolynomialSpace<C, A>(
     public override operator fun Polynomial<C>.unaryMinus(): Polynomial<C> = ring {
         Polynomial(coefficients.map { -it })
     }
+
     /**
      * Returns sum of the polynomials.
      */
@@ -210,6 +220,7 @@ public open class PolynomialSpace<C, A>(
             }
         )
     }
+
     /**
      * Returns difference of the polynomials.
      */
@@ -226,6 +237,7 @@ public open class PolynomialSpace<C, A>(
             }
         )
     }
+
     /**
      * Returns product of the polynomials.
      */
@@ -245,6 +257,7 @@ public open class PolynomialSpace<C, A>(
      * Instance of zero polynomial (zero of the polynomial ring).
      */
     override val zero: Polynomial<C> = Polynomial(emptyList())
+
     /**
      * Instance of unit polynomial (unit of the polynomial ring).
      */

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<Double> = sampler.sample(generator)
+    override fun sample(generator: RandomGenerator): BlockingDoubleChain = sampler.sample(generator)
 
     override fun cumulative(arg: Double): Double {
         val dev = arg - sampler.mean