From 6da51b7794f9ccf6d24c29b85b5081bba1cedf65 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sun, 9 Jul 2023 15:51:50 +0300
Subject: [PATCH] [WIP] Features to Attributes refactoring

---
 CHANGELOG.md                                  |    4 +-
 .../space/kscience/attributes/Attribute.kt    |   21 +-
 .../space/kscience/attributes/Attributes.kt   |   48 +-
 .../space/kscience/attributes/SafeType.kt     |   26 +
 .../space/kscience/attributes/annotations.kt  |   17 +
 buildSrc/settings.gradle.kts                  |    3 +-
 .../kmath/ejml/codegen/ejmlCodegen.kt         |    2 +
 gradle.properties                             |    3 +-
 gradle/wrapper/gradle-wrapper.properties      |    2 +-
 .../kscience/kmath/commons/linear/CMMatrix.kt |   31 +-
 .../kscience/kmath/expressions/DSAlgebra.kt   |    2 +
 .../kmath/linear/BufferedLinearSpace.kt       |    3 +
 .../kmath/linear/DoubleLinearSpace.kt         |    4 +-
 .../kscience/kmath/linear/LinearSpace.kt      |   40 +-
 .../kscience/kmath/linear/LupDecomposition.kt |  188 +--
 .../kscience/kmath/linear/MatrixBuilder.kt    |    7 +-
 .../kscience/kmath/linear/MatrixFeatures.kt   |  158 +--
 .../kscience/kmath/linear/MatrixWrapper.kt    |   47 +-
 .../kscience/kmath/linear/VirtualMatrix.kt    |    9 +-
 .../space/kscience/kmath/nd/AlgebraND.kt      |    4 +-
 .../space/kscience/kmath/nd/BufferND.kt       |    4 +-
 .../space/kscience/kmath/nd/Structure2D.kt    |    2 +-
 .../space/kscience/kmath/nd/StructureND.kt    |   22 +-
 .../kscience/kmath/operations/Algebra.kt      |   10 +-
 .../kmath/operations/BufferAlgebra.kt         |    4 +
 .../space/kscience/kmath/structures/Buffer.kt |   22 +-
 .../kmath/structures/BufferAccessor2D.kt      |    9 +-
 .../kmath/structures/MutableBuffer.kt         |   25 +-
 .../kscience/kmath/ejml/EjmlLinearSpace.kt    |    4 +-
 .../space/kscience/kmath/ejml/_generated.kt   | 1003 -----------------
 .../kscience/kmath/ejml/EjmlMatrixTest.kt     |    4 +-
 .../space/kscience/kmath/real/DoubleVector.kt |    1 -
 .../kscience/kmath/functions/Polynomial.kt    |   21 +-
 .../kmath/distributions/NormalDistribution.kt |    7 +-
 34 files changed, 457 insertions(+), 1300 deletions(-)
 create mode 100644 attributes-kt/src/commonMain/kotlin/space/kscience/attributes/SafeType.kt
 create mode 100644 attributes-kt/src/commonMain/kotlin/space/kscience/attributes/annotations.kt
 delete mode 100644 kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt

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<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> {
-    override fun equals(other: Any?): Boolean = (other as? PolymorphicAttribute<*>)?.type == this.type
+public abstract class PolymorphicAttribute<T>(public val type: SafeType<T>) : Attribute<T> {
+    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<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?,
-): Attributes = Attributes(
-    if (attrValue == null) {
-        content - attribute
-    } else {
-        content + (attribute to attrValue)
-    }
-)
+    attrValue: T,
+): Attributes = Attributes(content + (attribute to attrValue))
+
+public fun <A : Attribute<Unit>> 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 <T, A : SetAttribute<T>> Attributes.withoutAttributeElement(
     )
 }
 
+/**
+ * Create [Attributes] with a single key
+ */
 public fun <T : Any, A : Attribute<T>> 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<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>
\ 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<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 :
-                DeterminantFeature<Double>,
-                LupDecompositionFeature<Double> {
+            Determinant::class, LupDecompositionAttribute::class -> object :
+                Determinant<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 u: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(lup.u).withFeature(UFeature) }
+                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).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 r: Matrix<Double> by lazy<Matrix<Double>> { CMMatrix(qr.r).withFeature(UFeature) }
+                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).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)
     }
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<T, out A : Ring<T>>(
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<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()
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<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? =
-        structure.getFeature(type)
+    public fun <T, A : StructureAttribute<T>> attributeFor(structure: StructureND<*>, attribute: A): T? =
+        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)
 
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] &middot; a = [l] &middot; [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<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
-     */
-    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(LFeature)
+public data class LupDecomposition<T>(
+    public val l: Matrix<T>,
+    public val u: Matrix<T>,
+    public val p: Matrix<T>,
+)
 
 
-    /**
-     * 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(UFeature)
+public class LupDecompositionAttribute<T>(type: SafeType<LupDecomposition<T>>) :
+    PolymorphicAttribute<LupDecomposition<T>>(type),
+    MatrixAttribute<LupDecomposition<T>>
 
-    /**
-     * 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
-    }
+public val <T> MatrixOperations<T>.LUP: LupDecompositionAttribute<T>
+    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<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> =
-    lup(::DoubleBuffer, matrix) { it < singularityThreshold }
+): LupDecomposition<Double> = lup(matrix) { it < singularityThreshold }
 
-internal fun <T : Any> LupDecomposition<T>.solve(
-    factory: MutableBufferFactory<T>,
+internal fun <T : Any, A : Field<T>> LinearSpace<T, A>.solve(
+    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 {
-        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 <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)
-        return decomposition.solve(bufferFactory, b)
+        val decomposition = attributeFor(a, LUP) ?: lup(a, singularityCheck)
+        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 }
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 <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)
     }
 }
\ 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<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>> {
-    internal val instance: Inverted<Nothing> = Inverted()
-}
+public class Inverted<T>(type: SafeType<Matrix<T>>) :
+    PolymorphicAttribute<Matrix<T>>(type),
+    MatrixAttribute<Matrix<T>>
 
-@Suppress("UNCHECKED_CAST")
-public val <T> LinearSpace<T, *>.Inverted: Inverted<T> get() = Inverted.instance as Inverted<T>
+public val <T> MatrixOperations<T>.Inverted: Inverted<T> get() = Inverted(safeTypeOf())
 
 /**
  * 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>
 
-/**
- * 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
-}
+public inline val <reified T> MatrixOperations<T>.Determinant: Determinant<T> get() = Determinant(safeTypeOf())
 
 /**
  * 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 {
-    /**
-     * The lower triangular matrix in this decomposition. It may have [LFeature].
-     */
-    public val l: Matrix<T>
+public class LuDecompositionAttribute<T>(type: SafeType<LUDecomposition<T>>) :
+    PolymorphicAttribute<LUDecomposition<T>>(type),
+    MatrixAttribute<LUDecomposition<T>>
 
-    /**
-     * The upper triangular matrix in this decomposition. It may have [UFeature].
-     */
-    public val u: Matrix<T>
-}
+public val <T> MatrixOperations<T>.LU: LuDecompositionAttribute<T> get() = LuDecompositionAttribute(safeTypeOf())
 
-/**
- * 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
 
-/**
- * Matrices with this feature support QR factorization: *a = [q] &middot; [r]* where *a* is the owning matrix.
- *
- * @param T the type of matrices' items.
- */
-public interface QRDecompositionFeature<out T : Any> : MatrixFeature {
+
+public interface QRDecomposition<out T> {
     /**
-     * The orthogonal matrix in this decomposition. It may have [OrthogonalFeature].
+     * The orthogonal matrix in this decomposition. It may have [OrthogonalAttribute].
      */
     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 {
-    /**
-     * The triangular matrix in this decomposition. It may have either [UFeature] or [LFeature].
-     */
-    public val l: Matrix<T>
-}
+public class CholeskyDecompositionAttribute<T>(type: SafeType<CholeskyDecomposition<T>>) :
+    PolymorphicAttribute<CholeskyDecomposition<T>>(type),
+    MatrixAttribute<CholeskyDecomposition<T>>
 
-/**
- * Matrices with this feature support SVD: *a = [u] &middot; [s] &middot; [v]<sup>H</sup>* where *a* is the owning
- * matrix.
- *
- * @param T the type of matrices' items.
- */
-public interface SingularValueDecompositionFeature<out T : Any> : MatrixFeature {
+public val <T> MatrixOperations<T>.Cholesky: CholeskyDecompositionAttribute<T>
+    get() = CholeskyDecompositionAttribute(safeTypeOf())
+
+public interface SingularValueDecomposition<T> {
     /**
-     * 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<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
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<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) {
-    MatrixWrapper(origin, attributes.with(newFeature))
+public fun <T : Any, A : Attribute<T>> Matrix<T>.withAttribute(
+    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 operator fun <T : Any> Matrix<T>.plus(newFeature: MatrixFeature): MatrixWrapper<T> = withFeature(newFeature)
+public fun <T : Any, A : Attribute<Unit>> Matrix<T>.withAttribute(
+    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>
-    ?: VirtualMatrix(colNum, rowNum) { i, j -> get(j, i) }.withFeature(TransposedFeature(this))
+public fun <T : Any> Matrix<T>.transpose(): Matrix<T> = getFeature(TransposedAttribute::class)?.original as? Matrix<T>
+    ?: 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<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> =
-    VirtualMatrix(rows, columns, generator)
+public fun <T : Any> MatrixBuilder<T, *>.virtual(
+    attributes: Attributes = Attributes.EMPTY,
+    generator: (i: Int, j: Int) -> T,
+): VirtualMatrix<T> =
+    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<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)
 
 /**
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<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)
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<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()
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<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)
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<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:
      *
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<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> {
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<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> =
-            when (type) {
-                Double::class -> MutableBuffer.double(size) { initializer(it) as Double } as Buffer<T>
-                Short::class -> MutableBuffer.short(size) { initializer(it) as Short } as Buffer<T>
-                Int::class -> MutableBuffer.int(size) { initializer(it) as Int } as Buffer<T>
-                Long::class -> MutableBuffer.long(size) { initializer(it) as Long } as Buffer<T>
-                Float::class -> MutableBuffer.float(size) { initializer(it) as Float } as Buffer<T>
+        public inline fun <T> auto(type: SafeType<T>, size: Int, initializer: (Int) -> T): Buffer<T> =
+            when (type.kType) {
+                typeOf<Double>() -> MutableBuffer.double(size) { initializer(it) as Double } as Buffer<T>
+                typeOf<Short>() -> MutableBuffer.short(size) { initializer(it) as Short } as Buffer<T>
+                typeOf<Int>() -> MutableBuffer.int(size) { initializer(it) as Int } as Buffer<T>
+                typeOf<Long>() -> MutableBuffer.long(size) { initializer(it) as Long } 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> =
-            auto(T::class, size, initializer)
+        public inline fun <reified T> auto(size: Int, initializer: (Int) -> T): Buffer<T> =
+            auto(safeTypeOf<T>(), 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<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(
-        ColumnStrides(ShapeND(rowNum, colNum)),
-        factory
-    ) { (i, j) ->
-        get(i, j)
-    }.as2D()
+    fun MutableBuffer<T>.toStructure2D(): Structure2D<T> = StructureND.buffered(
+        ColumnStrides(ShapeND(rowNum, colNum))
+    ) { (i, j) -> get(i, j) }.as2D()
 
     inner class Row(val buffer: MutableBuffer<T>, val rowIndex: Int) : MutableBuffer<T> {
         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<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> =
-            when (type) {
-                Double::class -> double(size) { initializer(it) as Double } as MutableBuffer<T>
-                Short::class -> short(size) { initializer(it) as Short } as MutableBuffer<T>
-                Int::class -> int(size) { initializer(it) as Int } as MutableBuffer<T>
-                Float::class -> float(size) { initializer(it) as Float } as MutableBuffer<T>
-                Long::class -> long(size) { initializer(it) as Long } as MutableBuffer<T>
+        public inline fun <T> auto(type: SafeType<T>, size: Int, initializer: (Int) -> T): MutableBuffer<T> =
+            when (type.kType) {
+                typeOf<Double>() -> double(size) { initializer(it) as Double } as MutableBuffer<T>
+                typeOf<Short>() -> short(size) { initializer(it) as Short } as MutableBuffer<T>
+                typeOf<Int>() -> int(size) { initializer(it) as Int } as MutableBuffer<T>
+                typeOf<Float>() -> float(size) { initializer(it) as Float } as MutableBuffer<T>
+                typeOf<Long>() -> long(size) { initializer(it) as Long } as MutableBuffer<T>
                 else -> boxing(size, initializer)
             }
 
@@ -90,11 +92,10 @@ 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 : Any> auto(size: Int, initializer: (Int) -> T): MutableBuffer<T> =
-            auto(T::class, size, initializer)
+        public inline fun <reified T> auto(size: Int, initializer: (Int) -> T): MutableBuffer<T> =
+            auto(safeTypeOf<T>(), size, initializer)
     }
 }
 
 
-public sealed interface PrimitiveBuffer<T>: MutableBuffer<T>
\ No newline at end of file
+public sealed interface PrimitiveBuffer<T> : MutableBuffer<T>
\ 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<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>
 }
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<out M : DMatrix>(override val origin: M) : EjmlVector<Double, M>(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<out M : FMatrix>(override val origin: M) : EjmlVector<Float, M>(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<out M : DMatrix>(override val origin: M) : EjmlMatrix<Double, M>(origin) {
-    override operator fun get(i: Int, j: Int): Double = origin[i, j]
-}
-
-/**
- * [EjmlMatrix] specialization for [Float].
- */
-public class EjmlFloatMatrix<out M : FMatrix>(override val origin: M) : EjmlMatrix<Float, M>(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<Double, DoubleField, DMatrixRMaj>() {
-    /**
-     * The [DoubleField] reference.
-     */
-    override val elementAlgebra: DoubleField get() = DoubleField
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixRMaj> = when {
-        this is EjmlDoubleMatrix<*> && origin is DMatrixRMaj -> this as EjmlDoubleMatrix<DMatrixRMaj>
-        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixRMaj> = when {
-        this is EjmlDoubleVector<*> && origin is DMatrixRMaj -> this as EjmlDoubleVector<DMatrixRMaj>
-        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> = 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<DMatrixRMaj> = EjmlDoubleVector(DMatrixRMaj(size, 1).also {
-        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
-    })
-
-    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
-    private fun <T : DMatrix> T.wrapVector() = EjmlDoubleVector(this)
-
-    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
-
-    override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.mult(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.mult(toEjml().origin, vector.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        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<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixRMaj> {
-        val res = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.scale(value, toEjml().origin, res)
-        return res.wrapMatrix()
-    }
-
-    override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixRMaj> {
-        val res = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.changeSign(toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-        
-        CommonOps_DDRM.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin, 
-            out,
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override fun Point<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        val out = DMatrixRMaj(1, 1)
-
-        CommonOps_DDRM.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin,
-            out,
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Point<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        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<Double>): EjmlDoubleMatrix<DMatrixRMaj> = m * this
-
-    override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixRMaj> {
-        val res = DMatrixRMaj(1, 1)
-        CommonOps_DDRM.scale(value, toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixRMaj> = v * this
-
-    @UnstableKMathAPI
-    override fun <F : StructureFeature> computeFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
-        structure.getFeature(type)?.let { return it }
-        val origin = structure.toEjml().origin
-
-        return when (type) {
-            InverseMatrixFeature::class -> object : InverseMatrixFeature<Double> {
-                override val inverse: Matrix<Double> by lazy {
-                    val res = origin.copy()
-                    CommonOps_DDRM.invert(res)
-                    res.wrapMatrix()
-                }
-            }
-
-            DeterminantFeature::class -> object : DeterminantFeature<Double> {
-                override val determinant: Double by lazy { CommonOps_DDRM.det(origin) }
-            }
-
-            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Double> {
-                private val svd by lazy {
-                    DecompositionFactory_DDRM.svd(origin.numRows, origin.numCols, true, true, false)
-                        .apply { decompose(origin.copy()) }
-                }
-
-                override val u: Matrix<Double> by lazy { svd.getU(null, false).wrapMatrix() }
-                override val s: Matrix<Double> by lazy { svd.getW(null).wrapMatrix() }
-                override val v: Matrix<Double> by lazy { svd.getV(null, false).wrapMatrix() }
-                override val singularValues: Point<Double> by lazy { DoubleBuffer(svd.singularValues) }
-            }
-
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
-                private val qr by lazy {
-                    DecompositionFactory_DDRM.qr().apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<Double> by lazy {
-                    qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature)
-                }
-
-                override val r: Matrix<Double> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
-                override val l: Matrix<Double> by lazy {
-                    val cholesky =
-                        DecompositionFactory_DDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
-
-                    cholesky.getT(null).wrapMatrix().withFeature(LFeature)
-                }
-            }
-
-            LupDecompositionFeature::class -> object : LupDecompositionFeature<Double> {
-                private val lup by lazy {
-                    DecompositionFactory_DDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<Double> by lazy {
-                    lup.getLower(null).wrapMatrix().withFeature(LFeature)
-                }
-
-                override val u: Matrix<Double> by lazy {
-                    lup.getUpper(null).wrapMatrix().withFeature(UFeature)
-                }
-
-                override val p: Matrix<Double> by lazy { lup.getRowPivot(null).wrapMatrix() }
-            }
-
-            else -> null
-        }?.let{
-            type.cast(it)
-        }
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [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<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixRMaj> {
-        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] <sup>-1</sup> &middot; [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<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixRMaj> {
-        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<Float, FloatField, FMatrixRMaj>() {
-    /**
-     * The [FloatField] reference.
-     */
-    override val elementAlgebra: FloatField get() = FloatField
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Matrix<Float>.toEjml(): EjmlFloatMatrix<FMatrixRMaj> = when {
-        this is EjmlFloatMatrix<*> && origin is FMatrixRMaj -> this as EjmlFloatMatrix<FMatrixRMaj>
-        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Point<Float>.toEjml(): EjmlFloatVector<FMatrixRMaj> = when {
-        this is EjmlFloatVector<*> && origin is FMatrixRMaj -> this as EjmlFloatVector<FMatrixRMaj>
-        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> = 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<FMatrixRMaj> = EjmlFloatVector(FMatrixRMaj(size, 1).also {
-        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
-    })
-
-    private fun <T : FMatrix> T.wrapMatrix() = EjmlFloatMatrix(this)
-    private fun <T : FMatrix> T.wrapVector() = EjmlFloatVector(this)
-
-    override fun Matrix<Float>.unaryMinus(): Matrix<Float> = this * elementAlgebra { -one }
-
-    override fun Matrix<Float>.dot(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
-        val out = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.mult(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    override fun Matrix<Float>.dot(vector: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
-        val out = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.mult(toEjml().origin, vector.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    override operator fun Matrix<Float>.minus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
-        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<Float>.times(value: Float): EjmlFloatMatrix<FMatrixRMaj> {
-        val res = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.scale(value, toEjml().origin, res)
-        return res.wrapMatrix()
-    }
-
-    override fun Point<Float>.unaryMinus(): EjmlFloatVector<FMatrixRMaj> {
-        val res = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.changeSign(toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Matrix<Float>.plus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
-        val out = FMatrixRMaj(1, 1)
-        
-        CommonOps_FDRM.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin, 
-            out,
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override fun Point<Float>.plus(other: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
-        val out = FMatrixRMaj(1, 1)
-
-        CommonOps_FDRM.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin,
-            out,
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Point<Float>.minus(other: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
-        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<Float>): EjmlFloatMatrix<FMatrixRMaj> = m * this
-
-    override fun Point<Float>.times(value: Float): EjmlFloatVector<FMatrixRMaj> {
-        val res = FMatrixRMaj(1, 1)
-        CommonOps_FDRM.scale(value, toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixRMaj> = v * this
-
-    @UnstableKMathAPI
-    override fun <F : StructureFeature> computeFeature(structure: Matrix<Float>, type: KClass<out F>): F? {
-        structure.getFeature(type)?.let { return it }
-        val origin = structure.toEjml().origin
-
-        return when (type) {
-            InverseMatrixFeature::class -> object : InverseMatrixFeature<Float> {
-                override val inverse: Matrix<Float> by lazy {
-                    val res = origin.copy()
-                    CommonOps_FDRM.invert(res)
-                    res.wrapMatrix()
-                }
-            }
-
-            DeterminantFeature::class -> object : DeterminantFeature<Float> {
-                override val determinant: Float by lazy { CommonOps_FDRM.det(origin) }
-            }
-
-            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<Float> {
-                private val svd by lazy {
-                    DecompositionFactory_FDRM.svd(origin.numRows, origin.numCols, true, true, false)
-                        .apply { decompose(origin.copy()) }
-                }
-
-                override val u: Matrix<Float> by lazy { svd.getU(null, false).wrapMatrix() }
-                override val s: Matrix<Float> by lazy { svd.getW(null).wrapMatrix() }
-                override val v: Matrix<Float> by lazy { svd.getV(null, false).wrapMatrix() }
-                override val singularValues: Point<Float> by lazy { FloatBuffer(svd.singularValues) }
-            }
-
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<Float> {
-                private val qr by lazy {
-                    DecompositionFactory_FDRM.qr().apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<Float> by lazy {
-                    qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature)
-                }
-
-                override val r: Matrix<Float> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Float> {
-                override val l: Matrix<Float> by lazy {
-                    val cholesky =
-                        DecompositionFactory_FDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
-
-                    cholesky.getT(null).wrapMatrix().withFeature(LFeature)
-                }
-            }
-
-            LupDecompositionFeature::class -> object : LupDecompositionFeature<Float> {
-                private val lup by lazy {
-                    DecompositionFactory_FDRM.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<Float> by lazy {
-                    lup.getLower(null).wrapMatrix().withFeature(LFeature)
-                }
-
-                override val u: Matrix<Float> by lazy {
-                    lup.getUpper(null).wrapMatrix().withFeature(UFeature)
-                }
-
-                override val p: Matrix<Float> by lazy { lup.getRowPivot(null).wrapMatrix() }
-            }
-
-            else -> null
-        }?.let{
-            type.cast(it)
-        }
-    }
-
-    /**
-     * Solves for *x* in the following equation: *x = [a] <sup>-1</sup> &middot; [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<Float>, b: Matrix<Float>): EjmlFloatMatrix<FMatrixRMaj> {
-        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] <sup>-1</sup> &middot; [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<Float>, b: Point<Float>): EjmlFloatVector<FMatrixRMaj> {
-        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<Double, DoubleField, DMatrixSparseCSC>() {
-    /**
-     * The [DoubleField] reference.
-     */
-    override val elementAlgebra: DoubleField get() = DoubleField
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Matrix<Double>.toEjml(): EjmlDoubleMatrix<DMatrixSparseCSC> = when {
-        this is EjmlDoubleMatrix<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleMatrix<DMatrixSparseCSC>
-        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Point<Double>.toEjml(): EjmlDoubleVector<DMatrixSparseCSC> = when {
-        this is EjmlDoubleVector<*> && origin is DMatrixSparseCSC -> this as EjmlDoubleVector<DMatrixSparseCSC>
-        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> = 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<DMatrixSparseCSC> = EjmlDoubleVector(DMatrixSparseCSC(size, 1).also {
-        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
-    })
-
-    private fun <T : DMatrix> T.wrapMatrix() = EjmlDoubleMatrix(this)
-    private fun <T : DMatrix> T.wrapVector() = EjmlDoubleVector(this)
-
-    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = this * elementAlgebra { -one }
-
-    override fun Matrix<Double>.dot(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
-        val out = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.mult(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    override fun Matrix<Double>.dot(vector: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
-        val out = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.mult(toEjml().origin, vector.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    override operator fun Matrix<Double>.minus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
-        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<Double>.times(value: Double): EjmlDoubleMatrix<DMatrixSparseCSC> {
-        val res = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.scale(value, toEjml().origin, res)
-        return res.wrapMatrix()
-    }
-
-    override fun Point<Double>.unaryMinus(): EjmlDoubleVector<DMatrixSparseCSC> {
-        val res = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.changeSign(toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Matrix<Double>.plus(other: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
-        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<Double>.plus(other: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
-        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<Double>.minus(other: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
-        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<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> = m * this
-
-    override fun Point<Double>.times(value: Double): EjmlDoubleVector<DMatrixSparseCSC> {
-        val res = DMatrixSparseCSC(1, 1)
-        CommonOps_DSCC.scale(value, toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Double.times(v: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> = v * this
-
-    @UnstableKMathAPI
-    override fun <F : StructureFeature> computeFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
-        structure.getFeature(type)?.let { return it }
-        val origin = structure.toEjml().origin
-
-        return when (type) {
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<Double> {
-                private val qr by lazy {
-                    DecompositionFactory_DSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<Double> by lazy {
-                    qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature)
-                }
-
-                override val r: Matrix<Double> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
-                override val l: Matrix<Double> 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<Double>, DeterminantFeature<Double>, InverseMatrixFeature<Double> {
-                private val lu by lazy {
-                    DecompositionFactory_DSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<Double> by lazy {
-                    lu.getLower(null).wrapMatrix().withFeature(LFeature)
-                }
-
-                override val u: Matrix<Double> by lazy {
-                    lu.getUpper(null).wrapMatrix().withFeature(UFeature)
-                }
-
-                override val inverse: Matrix<Double> 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] <sup>-1</sup> &middot; [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<Double>, b: Matrix<Double>): EjmlDoubleMatrix<DMatrixSparseCSC> {
-        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] <sup>-1</sup> &middot; [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<Double>, b: Point<Double>): EjmlDoubleVector<DMatrixSparseCSC> {
-        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<Float, FloatField, FMatrixSparseCSC>() {
-    /**
-     * The [FloatField] reference.
-     */
-    override val elementAlgebra: FloatField get() = FloatField
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Matrix<Float>.toEjml(): EjmlFloatMatrix<FMatrixSparseCSC> = when {
-        this is EjmlFloatMatrix<*> && origin is FMatrixSparseCSC -> this as EjmlFloatMatrix<FMatrixSparseCSC>
-        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Point<Float>.toEjml(): EjmlFloatVector<FMatrixSparseCSC> = when {
-        this is EjmlFloatVector<*> && origin is FMatrixSparseCSC -> this as EjmlFloatVector<FMatrixSparseCSC>
-        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> = 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<FMatrixSparseCSC> = EjmlFloatVector(FMatrixSparseCSC(size, 1).also {
-        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
-    })
-
-    private fun <T : FMatrix> T.wrapMatrix() = EjmlFloatMatrix(this)
-    private fun <T : FMatrix> T.wrapVector() = EjmlFloatVector(this)
-
-    override fun Matrix<Float>.unaryMinus(): Matrix<Float> = this * elementAlgebra { -one }
-
-    override fun Matrix<Float>.dot(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
-        val out = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.mult(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    override fun Matrix<Float>.dot(vector: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
-        val out = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.mult(toEjml().origin, vector.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    override operator fun Matrix<Float>.minus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
-        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<Float>.times(value: Float): EjmlFloatMatrix<FMatrixSparseCSC> {
-        val res = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.scale(value, toEjml().origin, res)
-        return res.wrapMatrix()
-    }
-
-    override fun Point<Float>.unaryMinus(): EjmlFloatVector<FMatrixSparseCSC> {
-        val res = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.changeSign(toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Matrix<Float>.plus(other: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
-        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<Float>.plus(other: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
-        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<Float>.minus(other: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
-        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<Float>): EjmlFloatMatrix<FMatrixSparseCSC> = m * this
-
-    override fun Point<Float>.times(value: Float): EjmlFloatVector<FMatrixSparseCSC> {
-        val res = FMatrixSparseCSC(1, 1)
-        CommonOps_FSCC.scale(value, toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Float.times(v: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> = v * this
-
-    @UnstableKMathAPI
-    override fun <F : StructureFeature> computeFeature(structure: Matrix<Float>, type: KClass<out F>): F? {
-        structure.getFeature(type)?.let { return it }
-        val origin = structure.toEjml().origin
-
-        return when (type) {
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<Float> {
-                private val qr by lazy {
-                    DecompositionFactory_FSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<Float> by lazy {
-                    qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature)
-                }
-
-                override val r: Matrix<Float> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Float> {
-                override val l: Matrix<Float> 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<Float>, DeterminantFeature<Float>, InverseMatrixFeature<Float> {
-                private val lu by lazy {
-                    DecompositionFactory_FSCC.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<Float> by lazy {
-                    lu.getLower(null).wrapMatrix().withFeature(LFeature)
-                }
-
-                override val u: Matrix<Float> by lazy {
-                    lu.getUpper(null).wrapMatrix().withFeature(UFeature)
-                }
-
-                override val inverse: Matrix<Float> 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] <sup>-1</sup> &middot; [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<Float>, b: Matrix<Float>): EjmlFloatMatrix<FMatrixSparseCSC> {
-        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] <sup>-1</sup> &middot; [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<Float>, b: Point<Float>): EjmlFloatVector<FMatrixSparseCSC> {
-        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<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()) }
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<Double>
 
-@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<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).
      */
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<Double> = sampler.sample(generator)
+    override fun sample(generator: RandomGenerator): BlockingDoubleChain = sampler.sample(generator)
 
     override fun cumulative(arg: Double): Double {
         val dev = arg - sampler.mean