Prepare for 0.4.0 release

# Conflicts:
#	kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/_generated.kt

.gitignore

@@ -5,6 +5,7 @@ out/
 .idea/
 .vscode/
 .fleet/
+.kotlin/
 
 
 # Avoid ignoring Gradle wrapper jar file (.jar files are usually ignored)

CHANGELOG.md

@@ -3,14 +3,33 @@
 ## Unreleased
 
 ### Added
+
+### Changed
+
+### Deprecated
+
+### Removed
+
+### Fixed
+
+### Security
+
+## 0.4.0-dev-3 - 2024-02-18
+
+### Added
+
 - Reification. Explicit `SafeType` for algebras and buffers.
 - Integer division algebras.
 - Float32 geometries.
 - New Attributes-kt module that could be used as stand-alone. It declares. type-safe attributes containers.
 - Explicit `mutableStructureND` builders for mutable structures.
-- `Buffer.asList()` zero-copy transformation. 
+- `Buffer.asList()` zero-copy transformation.
+- Wasm support.
+- Parallel implementation of `LinearSpace` for Float64
+- Parallel buffer factories
 
 ### Changed
+
 - Default naming for algebra and buffers now uses IntXX/FloatXX notation instead of Java types.
 - Remove unnecessary inlines in basic algebras.
 - QuaternionField -> QuaternionAlgebra and does not implement `Field` anymore since it is non-commutative
@@ -21,18 +40,23 @@
 
 ### Deprecated
 
+- ND4J engine
+
 ### Removed
+
 - `asPolynomial` function due to scope pollution
+- Codegend for ejml (450 lines of codegen for 1000 lines of code is too much)
 
 ### Fixed
+
 - Median statistics
 - Complex power of negative real numbers
-
-### Security
+- Add proper mutability for MutableBufferND rows and columns
 
 ## 0.3.1 - 2023-04-09
 
 ### Added
+
 - Wasm support for `memory`, `core`, `complex` and `functions` modules.
 - Generic builders for `BufferND` and `MutableBufferND`
 - `NamedMatrix` - matrix with symbol-based indexing
@@ -42,6 +66,7 @@
 - Algebra now has an obligatory `bufferFactory` (#477).
 
 ### Changed
+
 - Removed marker `Vector` type for geometry
 - Geometry uses type-safe angles
 - Tensor operations switched to prefix notation
@@ -54,12 +79,14 @@
 - Multik went MPP
 
 ### Removed
+
 - Trajectory moved to https://github.com/SciProgCentre/maps-kt
 - Polynomials moved to https://github.com/SciProgCentre/kmath-polynomial
 
 ## 0.3.0
 
 ### Added
+
 - `ScaleOperations` interface
 - `Field` extends `ScaleOperations`
 - Basic integration API
@@ -84,6 +111,7 @@
 - Compilation to TeX for MST: #254
 
 ### Changed
+
 - Annotations moved to `space.kscience.kmath`
 - Exponential operations merged with hyperbolic functions
 - Space is replaced by Group. Space is reserved for vector spaces.
@@ -117,9 +145,11 @@
 - `UnivariateFunction` -> `Function1D`, `MultivariateFunction` -> `FunctionND`
 
 ### Deprecated
+
 - Specialized `DoubleBufferAlgebra`
 
 ### Removed
+
 - Nearest in Domain. To be implemented in geometry package.
 - Number multiplication and division in main Algebra chain
 - `contentEquals` from Buffer. It moved to the companion.
@@ -130,12 +160,14 @@
 - Algebra elements are completely removed. Use algebra contexts instead.
 
 ### Fixed
+
 - Ring inherits RingOperations, not GroupOperations
 - Univariate histogram filling
 
 ## 0.2.0
 
 ### Added
+
 - `fun` annotation for SAM interfaces in library
 - Explicit `public` visibility for all public APIs
 - Better trigonometric and hyperbolic functions for `AutoDiffField` (https://github.com/mipt-npm/kmath/pull/140)
@@ -155,6 +187,7 @@
 - Basic Quaternion vector support in `kmath-complex`.
 
 ### Changed
+
 - Package changed from `scientifik` to `space.kscience`
 - Gradle version: 6.6 -> 6.8.2
 - Minor exceptions refactor (throwing `IllegalArgumentException` by argument checks instead of `IllegalStateException`)
@@ -179,6 +212,7 @@
 - Add `out` projection to `Buffer` generic
 
 ### Removed
+
 - `kmath-koma` module because it doesn't support Kotlin 1.4.
 - Support of `legacy` JS backend (we will support only IR)
 - `toGrid` method.
@@ -187,11 +221,13 @@
 - StructureND identity and equals
 
 ### Fixed
+
 - `symbol` method in `MstExtendedField` (https://github.com/mipt-npm/kmath/pull/140)
 
 ## 0.1.4
 
 ### Added
+
 - Functional Expressions API
 - Mathematical Syntax Tree, its interpreter and API
 - String to MST parser (https://github.com/mipt-npm/kmath/pull/120)
@@ -209,6 +245,7 @@
 - Norm support for `Complex`
 
 ### Changed
+
 - `readAsMemory` now has `throws IOException` in JVM signature.
 - Several functions taking functional types were made `inline`.
 - Several functions taking functional types now have `callsInPlace` contracts.
@@ -220,6 +257,7 @@
 - Moved probability distributions to commons-rng and to `kmath-prob`
 
 ### Fixed
+
 - Missing copy method in Memory implementation on JS (https://github.com/mipt-npm/kmath/pull/106)
 - D3.dim value in `kmath-dimensions`
 - Multiplication in integer rings in `kmath-core` (https://github.com/mipt-npm/kmath/pull/101)

README.md

@@ -2,7 +2,7 @@
 [![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382)
 ![Gradle build](https://github.com/SciProgCentre/kmath/workflows/Gradle%20build/badge.svg)
 [![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22)
-[![Space](https://img.shields.io/badge/dynamic/xml?color=orange&label=Space&query=//metadata/versioning/latest&url=https%3A%2F%2Fmaven.pkg.jetbrains.space%2Fmipt-npm%2Fp%2Fsci%2Fmaven%2Fspace%2Fkscience%2Fkmath-core%2Fmaven-metadata.xml)](https://maven.pkg.jetbrains.space/spc/p/sci/maven/space/kscience/)
+[![Space](https://img.shields.io/badge/dynamic/xml?color=orange&label=Space&query=//metadata/versioning/latest&url=https%3A%2F%2Fmaven.pkg.jetbrains.space%2Fmipt-npm%2Fp%2Fsci%2Fmaven%2Fspace%2Fkscience%2Fkmath-core%2Fmaven-metadata.xml)](https://maven.pkg.jetbrains.space/mipt-npm/p/sci/maven/space/kscience/)
 
 # KMath
 
@@ -11,18 +11,21 @@ analog to Python's NumPy library. Later we found that kotlin is much more flexib
 architecture designs. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. The `numpy`-like
 experience could be achieved with [kmath-for-real](/kmath-for-real) extension module.
 
-[Documentation site (**WIP**)](https://SciProgCentre.github.io/kmath/)
+[Documentation site](https://SciProgCentre.github.io/kmath/)
 
 ## Publications and talks
 
 * [A conceptual article about context-oriented design](https://proandroiddev.com/an-introduction-context-oriented-programming-in-kotlin-2e79d316b0a2)
 * [Another article about context-oriented design](https://proandroiddev.com/diving-deeper-into-context-oriented-programming-in-kotlin-3ecb4ec38814)
 * [ACAT 2019 conference paper](https://aip.scitation.org/doi/abs/10.1063/1.5130103)
+* [A talk at KotlinConf 2019 about using kotlin for science](https://youtu.be/LI_5TZ7tnOE?si=4LknX41gl_YeUbIe)
+* [A talk on architecture at Joker-2021 (in Russian)](https://youtu.be/1bZ2doHiRRM?si=9w953ro9yu98X_KJ)
+* [The same talk in English](https://youtu.be/yP5DIc2fVwQ?si=louZzQ1dcXV6gP10)
+* [A seminar on tensor API](https://youtu.be/0H99wUs0xTM?si=6c__04jrByFQtVpo)
 
 # Goal
 
-* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM, JS and Native)
-  .
+* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM, JS, Native and Wasm).
 * Provide basic multiplatform implementations for those abstractions (without significant performance optimization).
 * Provide bindings and wrappers with those abstractions for popular optimized platform libraries.
 
@@ -53,18 +56,20 @@ module definitions below. The module stability could have the following levels:
 ## Modules
 
 
+### [attributes-kt](attributes-kt)
+> An API and basic implementation for arranging objects in a continuous memory block.
+>
+> **Maturity**: DEVELOPMENT
+
 ### [benchmarks](benchmarks)
-> 
 >
 > **Maturity**: EXPERIMENTAL
 
 ### [examples](examples)
-> 
 >
 > **Maturity**: EXPERIMENTAL
 
 ### [kmath-ast](kmath-ast)
-> 
 >
 > **Maturity**: EXPERIMENTAL
 >
@@ -76,7 +81,7 @@ module definitions below. The module stability could have the following levels:
 
 
 ### [kmath-commons](kmath-commons)
-> 
+> Commons math binding for kmath
 >
 > **Maturity**: EXPERIMENTAL
 
@@ -108,17 +113,15 @@ performance calculations to code generation.
 
 
 ### [kmath-coroutines](kmath-coroutines)
-> 
 >
 > **Maturity**: EXPERIMENTAL
 
 ### [kmath-dimensions](kmath-dimensions)
-> 
+> A proof of concept module for adding type-safe dimensions to structures
 >
 > **Maturity**: PROTOTYPE
 
 ### [kmath-ejml](kmath-ejml)
-> 
 >
 > **Maturity**: PROTOTYPE
 >
@@ -142,7 +145,7 @@ One can still use generic algebras though.
 
 
 ### [kmath-functions](kmath-functions)
-> 
+> Functions, integration and interpolation
 >
 > **Maturity**: EXPERIMENTAL
 >
@@ -155,31 +158,28 @@ One can still use generic algebras though.
 
 
 ### [kmath-geometry](kmath-geometry)
-> 
 >
 > **Maturity**: PROTOTYPE
 
 ### [kmath-histograms](kmath-histograms)
-> 
 >
 > **Maturity**: PROTOTYPE
 
 ### [kmath-jafama](kmath-jafama)
-> 
+> Jafama integration module
 >
-> **Maturity**: PROTOTYPE
+> **Maturity**: DEPRECATED
 >
 > **Features:**
 > - [jafama-double](kmath-jafama/src/main/kotlin/space/kscience/kmath/jafama/) : Double ExtendedField implementations based on Jafama
 
 
 ### [kmath-jupyter](kmath-jupyter)
-> 
 >
 > **Maturity**: PROTOTYPE
 
 ### [kmath-kotlingrad](kmath-kotlingrad)
-> 
+> Kotlin∇ integration module
 >
 > **Maturity**: EXPERIMENTAL
 >
@@ -194,14 +194,14 @@ One can still use generic algebras though.
 > **Maturity**: DEVELOPMENT
 
 ### [kmath-multik](kmath-multik)
-> 
+> JetBrains Multik connector
 >
 > **Maturity**: PROTOTYPE
 
 ### [kmath-nd4j](kmath-nd4j)
-> 
+> ND4J NDStructure implementation and according NDAlgebra classes
 >
-> **Maturity**: EXPERIMENTAL
+> **Maturity**: DEPRECATED
 >
 > **Features:**
 > - [nd4jarraystructure](kmath-nd4j/#) : NDStructure wrapper for INDArray
@@ -210,27 +210,24 @@ One can still use generic algebras though.
 
 
 ### [kmath-optimization](kmath-optimization)
-> 
 >
 > **Maturity**: EXPERIMENTAL
 
 ### [kmath-stat](kmath-stat)
-> 
 >
 > **Maturity**: EXPERIMENTAL
 
 ### [kmath-symja](kmath-symja)
-> 
+> Symja integration module
 >
 > **Maturity**: PROTOTYPE
 
 ### [kmath-tensorflow](kmath-tensorflow)
-> 
+> Google tensorflow connector
 >
 > **Maturity**: PROTOTYPE
 
 ### [kmath-tensors](kmath-tensors)
-> 
 >
 > **Maturity**: PROTOTYPE
 >
@@ -241,12 +238,11 @@ One can still use generic algebras though.
 
 
 ### [kmath-viktor](kmath-viktor)
-> 
+> Binding for https://github.com/JetBrains-Research/viktor
 >
 > **Maturity**: DEVELOPMENT
 
 ### [test-utils](test-utils)
-> 
 >
 > **Maturity**: EXPERIMENTAL
 
@@ -256,22 +252,21 @@ One can still use generic algebras though.
 KMath is developed as a multi-platform library, which means that most of the interfaces are declared in the
 [common source sets](/kmath-core/src/commonMain) and implemented there wherever it is possible. In some cases, features
 are delegated to platform-specific implementations even if they could be provided in the common module for performance
-reasons. Currently, the Kotlin/JVM is the primary platform, however Kotlin/Native and Kotlin/JS contributions and
+reasons. Currently, Kotlin/JVM is the primary platform, however, Kotlin/Native and Kotlin/JS contributions and
 feedback are also welcome.
 
 ## Performance
 
-Calculation performance is one of major goals of KMath in the future, but in some cases it is impossible to achieve both
+Calculation of performance is one of the major goals of KMath in the future, but in some cases it is impossible to achieve both
 performance and flexibility.
 
-We expect to focus on creating convenient universal API first and then work on increasing performance for specific
+We expect to focus on creating a convenient universal API first and then work on increasing performance for specific
 cases. We expect the worst KMath benchmarks will perform better than native Python, but worse than optimized
 native/SciPy (mostly due to boxing operations on primitive numbers). The best performance of optimized parts could be better than SciPy.
 
 ## Requirements
 
-KMath currently relies on JDK 11 for compilation and execution of Kotlin-JVM part. We recommend to use GraalVM-CE 11 for
-execution to get better performance.
+KMath currently relies on JDK 11 for compilation and execution of Kotlin-JVM part. We recommend using GraalVM-CE or Oracle GraalVM for execution to get better performance.
 
 ### Repositories
 
@@ -291,10 +286,7 @@ dependencies {
 }
 ```
 
-Gradle `6.0+` is required for multiplatform artifacts.
-
 ## Contributing
 
-The project requires a lot of additional work. The most important thing we need is a feedback about what features are
-required the most. Feel free to create feature requests. We are also welcome to code contributions, especially in issues
-marked with [waiting for a hero](https://github.com/SciProgCentre/kmath/labels/waiting%20for%20a%20hero) label.
+The project requires a lot of additional work. The most important thing we need is feedback about what features are
+required the most. Feel free to create feature requests. We are also welcome to code contributions, especially in issues marked with [good first issue](hhttps://github.com/SciProgCentre/kmath/issues?q=is%3Aissue+is%3Aopen+label%3A%22good+first+issue%22) label.

attributes-kt/README.md

@@ -0,0 +1,4 @@
+# Module attributes-kt
+
+
+

attributes-kt/api/attributes-kt.api

@@ -0,0 +1,97 @@
+public abstract interface class space/kscience/attributes/Attribute {
+}
+
+public abstract interface class space/kscience/attributes/AttributeContainer {
+	public abstract fun getAttributes ()Lspace/kscience/attributes/Attributes;
+}
+
+public abstract interface class space/kscience/attributes/AttributeScope {
+}
+
+public abstract interface class space/kscience/attributes/AttributeWithDefault : space/kscience/attributes/Attribute {
+	public abstract fun getDefault ()Ljava/lang/Object;
+}
+
+public abstract interface class space/kscience/attributes/Attributes {
+	public static final field Companion Lspace/kscience/attributes/Attributes$Companion;
+	public fun get (Lspace/kscience/attributes/Attribute;)Ljava/lang/Object;
+	public abstract fun getContent ()Ljava/util/Map;
+	public fun getKeys ()Ljava/util/Set;
+}
+
+public final class space/kscience/attributes/Attributes$Companion {
+	public final fun getEMPTY ()Lspace/kscience/attributes/Attributes;
+}
+
+public final class space/kscience/attributes/AttributesBuilder : space/kscience/attributes/Attributes {
+	public fun <init> ()V
+	public final fun add (Lspace/kscience/attributes/SetAttribute;Ljava/lang/Object;)V
+	public final fun build ()Lspace/kscience/attributes/Attributes;
+	public final fun from (Lspace/kscience/attributes/Attributes;)V
+	public fun getContent ()Ljava/util/Map;
+	public final fun invoke (Lspace/kscience/attributes/Attribute;Ljava/lang/Object;)V
+	public final fun remove (Lspace/kscience/attributes/SetAttribute;Ljava/lang/Object;)V
+	public final fun set (Lspace/kscience/attributes/Attribute;Ljava/lang/Object;)V
+}
+
+public final class space/kscience/attributes/AttributesBuilderKt {
+	public static final fun Attributes (Lkotlin/jvm/functions/Function1;)Lspace/kscience/attributes/Attributes;
+}
+
+public final class space/kscience/attributes/AttributesKt {
+	public static final fun Attributes (Lspace/kscience/attributes/Attribute;)Lspace/kscience/attributes/Attributes;
+	public static final fun Attributes (Lspace/kscience/attributes/Attribute;Ljava/lang/Object;)Lspace/kscience/attributes/Attributes;
+	public static final fun getOrDefault (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/AttributeWithDefault;)Ljava/lang/Object;
+	public static final fun isEmpty (Lspace/kscience/attributes/Attributes;)Z
+	public static final fun modify (Lspace/kscience/attributes/Attributes;Lkotlin/jvm/functions/Function1;)Lspace/kscience/attributes/Attributes;
+	public static final fun plus (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/Attributes;)Lspace/kscience/attributes/Attributes;
+	public static final fun withAttribute (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/Attribute;)Lspace/kscience/attributes/Attributes;
+	public static final fun withAttribute (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/Attribute;Ljava/lang/Object;)Lspace/kscience/attributes/Attributes;
+	public static final fun withAttributeElement (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/SetAttribute;Ljava/lang/Object;)Lspace/kscience/attributes/Attributes;
+	public static final fun withoutAttribute (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/Attribute;)Lspace/kscience/attributes/Attributes;
+	public static final fun withoutAttributeElement (Lspace/kscience/attributes/Attributes;Lspace/kscience/attributes/SetAttribute;Ljava/lang/Object;)Lspace/kscience/attributes/Attributes;
+}
+
+public abstract interface class space/kscience/attributes/FlagAttribute : space/kscience/attributes/Attribute {
+}
+
+public abstract class space/kscience/attributes/PolymorphicAttribute : space/kscience/attributes/Attribute {
+	public synthetic fun <init> (Lkotlin/reflect/KType;Lkotlin/jvm/internal/DefaultConstructorMarker;)V
+	public fun equals (Ljava/lang/Object;)Z
+	public final fun getType-V0oMfBY ()Lkotlin/reflect/KType;
+	public fun hashCode ()I
+}
+
+public final class space/kscience/attributes/PolymorphicAttributeKt {
+	public static final fun get (Lspace/kscience/attributes/Attributes;Lkotlin/jvm/functions/Function0;)Ljava/lang/Object;
+	public static final fun set (Lspace/kscience/attributes/AttributesBuilder;Lkotlin/jvm/functions/Function0;Ljava/lang/Object;)V
+}
+
+public final class space/kscience/attributes/SafeType {
+	public static final synthetic fun box-impl (Lkotlin/reflect/KType;)Lspace/kscience/attributes/SafeType;
+	public static fun constructor-impl (Lkotlin/reflect/KType;)Lkotlin/reflect/KType;
+	public fun equals (Ljava/lang/Object;)Z
+	public static fun equals-impl (Lkotlin/reflect/KType;Ljava/lang/Object;)Z
+	public static final fun equals-impl0 (Lkotlin/reflect/KType;Lkotlin/reflect/KType;)Z
+	public final fun getKType ()Lkotlin/reflect/KType;
+	public fun hashCode ()I
+	public static fun hashCode-impl (Lkotlin/reflect/KType;)I
+	public fun toString ()Ljava/lang/String;
+	public static fun toString-impl (Lkotlin/reflect/KType;)Ljava/lang/String;
+	public final synthetic fun unbox-impl ()Lkotlin/reflect/KType;
+}
+
+public final class space/kscience/attributes/SafeTypeKt {
+	public static final fun getKClass-X0YbwmU (Lkotlin/reflect/KType;)Lkotlin/reflect/KClass;
+}
+
+public abstract interface class space/kscience/attributes/SetAttribute : space/kscience/attributes/Attribute {
+}
+
+public abstract interface annotation class space/kscience/attributes/UnstableAttributesAPI : java/lang/annotation/Annotation {
+}
+
+public abstract interface class space/kscience/attributes/WithType {
+	public abstract fun getType-V0oMfBY ()Lkotlin/reflect/KType;
+}
+

attributes-kt/build.gradle.kts

@@ -2,6 +2,8 @@ plugins {
     id("space.kscience.gradle.mpp")
 }
 
+version = "0.1.0"
+
 kscience {
     jvm()
     js()

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

@@ -5,6 +5,9 @@
 
 package space.kscience.attributes
 
+/**
+ * A marker interface for an attribute. Attributes are used as keys to access contents of type [T] in the container.
+ */
 public interface Attribute<T>
 
 /**

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

@@ -13,7 +13,8 @@ public interface AttributeContainer {
 }
 
 /**
- * A scope, where attribute keys could be resolved
+ * A scope, where attribute keys could be resolved.
+ * [O] is used only to resolve types in compile-time.
  */
 public interface AttributeScope<O>
 

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

@@ -11,10 +11,19 @@ import kotlin.jvm.JvmInline
  * A set of attributes. The implementation must guarantee that [content] keys correspond to its value types.
  */
 public interface Attributes {
+    /**
+     * Raw content for this [Attributes]
+     */
     public val content: Map<out Attribute<*>, Any?>
 
+    /**
+     * Attribute keys contained in this [Attributes]
+     */
     public val keys: Set<Attribute<*>> get() = content.keys
 
+    /**
+     * Provide an attribute value. Return null if attribute is not present or if its value is null.
+     */
     @Suppress("UNCHECKED_CAST")
     public operator fun <T> get(attribute: Attribute<T>): T? = content[attribute] as? T
 

benchmarks/src/jsMain/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt

@@ -99,7 +99,7 @@ class ExpressionsInterpretersBenchmark {
         private val estree = node.estreeCompileToExpression(Float64Field)
 
         private val raw = Expression<Double> { args ->
-            val x = args[x]!!
+            val x = args.getValue(x)
             x * 2.0 + 2.0 / x - 16.0 / sin(x)
         }
     }

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/DotBenchmark.kt

@@ -11,12 +11,11 @@ import kotlinx.benchmark.Scope
 import kotlinx.benchmark.State
 import space.kscience.kmath.commons.linear.CMLinearSpace
 import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
+import space.kscience.kmath.linear.Float64ParallelLinearSpace
 import space.kscience.kmath.linear.invoke
 import space.kscience.kmath.linear.linearSpace
 import space.kscience.kmath.operations.Float64Field
-import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.tensorflow.produceWithTF
-import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
 import space.kscience.kmath.tensors.core.tensorAlgebra
 import kotlin.random.Random
 
@@ -72,12 +71,12 @@ internal class DotBenchmark {
     }
 
     @Benchmark
-    fun tensorDot(blackhole: Blackhole) = with(Float64Field.tensorAlgebra) {
+    fun multikDot(blackhole: Blackhole) = with(multikAlgebra) {
         blackhole.consume(matrix1 dot matrix2)
     }
 
     @Benchmark
-    fun multikDot(blackhole: Blackhole) = with(multikAlgebra) {
+    fun tensorDot(blackhole: Blackhole) = with(Float64Field.tensorAlgebra) {
         blackhole.consume(matrix1 dot matrix2)
     }
 
@@ -87,12 +86,8 @@ internal class DotBenchmark {
     }
 
     @Benchmark
-    fun doubleDot(blackhole: Blackhole) = with(Float64Field.linearSpace) {
+    fun parallelDot(blackhole: Blackhole) = with(Float64ParallelLinearSpace) {
         blackhole.consume(matrix1 dot matrix2)
     }
 
-    @Benchmark
-    fun doubleTensorDot(blackhole: Blackhole) = DoubleTensorAlgebra.invoke {
-        blackhole.consume(matrix1 dot matrix2)
-    }
 }

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/MatrixInverseBenchmark.kt

@@ -15,6 +15,7 @@ import space.kscience.kmath.ejml.EjmlLinearSpaceDDRM
 import space.kscience.kmath.linear.invoke
 import space.kscience.kmath.linear.linearSpace
 import space.kscience.kmath.linear.lupSolver
+import space.kscience.kmath.linear.parallel
 import space.kscience.kmath.operations.algebra
 import kotlin.random.Random
 
@@ -38,16 +39,19 @@ internal class MatrixInverseBenchmark {
     }
 
     @Benchmark
-    fun cmLUPInversion(blackhole: Blackhole) {
-        CMLinearSpace {
-            blackhole.consume(lupSolver().inverse(matrix))
-        }
+    fun kmathParallelLupInversion(blackhole: Blackhole) {
+        blackhole.consume(Double.algebra.linearSpace.parallel.lupSolver().inverse(matrix))
     }
 
     @Benchmark
-    fun ejmlInverse(blackhole: Blackhole) {
-        EjmlLinearSpaceDDRM {
-            blackhole.consume(matrix.toEjml().inverted())
-        }
+    fun cmLUPInversion(blackhole: Blackhole) = CMLinearSpace {
+        blackhole.consume(lupSolver().inverse(matrix))
     }
+
+
+    @Benchmark
+    fun ejmlInverse(blackhole: Blackhole) = EjmlLinearSpaceDDRM {
+        blackhole.consume(matrix.toEjml().inverted())
+    }
+
 }

build.gradle.kts

@@ -3,7 +3,7 @@ import space.kscience.gradle.useSPCTeam
 
 plugins {
     id("space.kscience.gradle.project")
-    id("org.jetbrains.kotlinx.kover") version "0.6.0"
+    id("org.jetbrains.kotlinx.kover") version "0.7.6"
 }
 
 allprojects {
@@ -14,7 +14,7 @@ allprojects {
     }
 
     group = "space.kscience"
-    version = "0.4.0-dev-3"
+    version = "0.4.0"
 }
 
 subprojects {
@@ -34,7 +34,7 @@ subprojects {
                     localDirectory.set(kotlinDir)
 
                     remoteUrl.set(
-                        java.net.URL("https://github.com/SciProgCentre/kmath/tree/master/${this@subprojects.name}/$kotlinDirPath")
+                        uri("https://github.com/SciProgCentre/kmath/tree/master/${this@subprojects.name}/$kotlinDirPath").toURL()
                     )
                 }
 

buildSrc/build.gradle.kts

@@ -17,7 +17,7 @@ val benchmarksVersion = spclibs.versions.kotlinx.benchmark.get()
 dependencies {
     api("space.kscience:gradle-tools:$toolsVersion")
     //plugins form benchmarks
-    api("org.jetbrains.kotlinx:kotlinx-benchmark-plugin:0.4.9")
+    api("org.jetbrains.kotlinx:kotlinx-benchmark-plugin:$benchmarksVersion")
     //api("org.jetbrains.kotlin:kotlin-allopen:$kotlinVersion")
     //to be used inside build-script only
     //implementation(spclibs.kotlinx.serialization.json)

buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt

@@ -1,443 +0,0 @@
-/*
- * Copyright 2018-2024 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.
- */
-
-@file:Suppress("KDocUnresolvedReference")
-
-package space.kscience.kmath.ejml.codegen
-
-import org.intellij.lang.annotations.Language
-import java.io.File
-
-private fun Appendable.appendEjmlVector(type: String, ejmlMatrixType: String) {
-    @Language("kotlin") val text = """/**
- * [EjmlVector] specialization for [$type].
- */
-public class Ejml${type}Vector<out M : $ejmlMatrixType>(override val origin: M) : EjmlVector<$type, M>(origin) {
-    init {
-        require(origin.numRows == 1) { "The origin matrix must have only one row to form a vector" }
-    }
-
-    override val type: SafeType<${type}> get() = safeTypeOf()
-
-    override operator fun get(index: Int): $type = origin[0, index]
-}"""
-    appendLine(text)
-    appendLine()
-}
-
-private fun Appendable.appendEjmlMatrix(type: String, ejmlMatrixType: String) {
-    val text = """/**
- * [EjmlMatrix] specialization for [$type].
- */
-public class Ejml${type}Matrix<out M : $ejmlMatrixType>(override val origin: M) : EjmlMatrix<$type, M>(origin) {
-    override val type: SafeType<${type}> get() = safeTypeOf()
-
-    override operator fun get(i: Int, j: Int): $type = origin[i, j]
-}"""
-    appendLine(text)
-    appendLine()
-}
-
-private fun Appendable.appendEjmlLinearSpace(
-    type: String,
-    kmathAlgebra: String,
-    ejmlMatrixParentTypeMatrix: String,
-    ejmlMatrixType: String,
-    ejmlMatrixDenseType: String,
-    ops: String,
-    denseOps: String,
-    isDense: Boolean,
-) {
-    @Language("kotlin") val text = """
-
-/**
- * [EjmlLinearSpace] implementation based on [CommonOps_$ops], [DecompositionFactory_${ops}] operations and
- * [${ejmlMatrixType}] matrices.
- */
-public object EjmlLinearSpace${ops} : EjmlLinearSpace<${type}, ${kmathAlgebra}, $ejmlMatrixType>() {
-    /**
-     * The [${kmathAlgebra}] reference.
-     */
-    override val elementAlgebra: $kmathAlgebra get() = $kmathAlgebra
-
-    override val type: SafeType<${type}> get() = safeTypeOf()
-
-    @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}>
-        else -> buildMatrix(rowNum, colNum) { i, j -> get(i, j) }
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    override fun Point<${type}>.toEjml(): Ejml${type}Vector<${ejmlMatrixType}> = when {
-        this is Ejml${type}Vector<*> && origin is $ejmlMatrixType -> this as Ejml${type}Vector<${ejmlMatrixType}>
-        else -> Ejml${type}Vector(${ejmlMatrixType}(size, 1).also {
-            (0 until it.numRows).forEach { row -> it[row, 0] = get(row) }
-        })
-    }
-
-    override fun buildMatrix(
-        rows: Int,
-        columns: Int,
-        initializer: ${kmathAlgebra}.(i: Int, j: Int) -> ${type},
-    ): Ejml${type}Matrix<${ejmlMatrixType}> = ${ejmlMatrixType}(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: ${kmathAlgebra}.(Int) -> ${type},
-    ): Ejml${type}Vector<${ejmlMatrixType}> = Ejml${type}Vector(${ejmlMatrixType}(size, 1).also {
-        (0 until it.numRows).forEach { row -> it[row, 0] = elementAlgebra.initializer(row) }
-    })
-
-    private fun <T : ${ejmlMatrixParentTypeMatrix}> T.wrapMatrix() = Ejml${type}Matrix(this)
-    private fun <T : ${ejmlMatrixParentTypeMatrix}> T.wrapVector() = Ejml${type}Vector(this)
-
-    override fun Matrix<${type}>.unaryMinus(): Matrix<${type}> = this * elementAlgebra { -one }
-
-    override fun Matrix<${type}>.dot(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
-        val out = ${ejmlMatrixType}(1, 1)
-        CommonOps_${ops}.mult(toEjml().origin, other.toEjml().origin, out)
-        return out.wrapMatrix()
-    }
-
-    override fun Matrix<${type}>.dot(vector: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
-        val out = ${ejmlMatrixType}(1, 1)
-        CommonOps_${ops}.mult(toEjml().origin, vector.toEjml().origin, out)
-        return out.wrapVector()
-    }
-
-    override operator fun Matrix<${type}>.minus(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
-        val out = ${ejmlMatrixType}(1, 1)
-
-        CommonOps_${ops}.add(
-            elementAlgebra.one, 
-            toEjml().origin, 
-            elementAlgebra { -one },
-            other.toEjml().origin, 
-            out,${
-        if (isDense) "" else
-            """
-            null, 
-            null,"""
-    }
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override operator fun Matrix<${type}>.times(value: ${type}): Ejml${type}Matrix<${ejmlMatrixType}> {
-        val res = ${ejmlMatrixType}(1, 1)
-        CommonOps_${ops}.scale(value, toEjml().origin, res)
-        return res.wrapMatrix()
-    }
-
-    override fun Point<${type}>.unaryMinus(): Ejml${type}Vector<${ejmlMatrixType}> {
-        val res = ${ejmlMatrixType}(1, 1)
-        CommonOps_${ops}.changeSign(toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun Matrix<${type}>.plus(other: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
-        val out = ${ejmlMatrixType}(1, 1)
-        
-        CommonOps_${ops}.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin, 
-            out,${
-        if (isDense) "" else
-            """
-            null, 
-            null,"""
-    }
-        )
-
-        return out.wrapMatrix()
-    }
-
-    override fun Point<${type}>.plus(other: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
-        val out = ${ejmlMatrixType}(1, 1)
-
-        CommonOps_${ops}.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra.one,
-            other.toEjml().origin,
-            out,${
-        if (isDense) "" else
-            """
-            null, 
-            null,"""
-    }
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun Point<${type}>.minus(other: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
-        val out = ${ejmlMatrixType}(1, 1)
-
-        CommonOps_${ops}.add(
-            elementAlgebra.one,
-            toEjml().origin,
-            elementAlgebra { -one },
-            other.toEjml().origin,
-            out,${
-        if (isDense) "" else
-            """
-            null, 
-            null,"""
-    }
-        )
-
-        return out.wrapVector()
-    }
-
-    override fun ${type}.times(m: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> = m * this
-
-    override fun Point<${type}>.times(value: ${type}): Ejml${type}Vector<${ejmlMatrixType}> {
-        val res = ${ejmlMatrixType}(1, 1)
-        CommonOps_${ops}.scale(value, toEjml().origin, res)
-        return res.wrapVector()
-    }
-
-    override fun ${type}.times(v: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> = v * this
-
-    @UnstableKMathAPI
-    override fun <F : StructureFeature> computeFeature(structure: Matrix<${type}>, type: KClass<out F>): F? {
-        structure.getFeature(type)?.let { return it }
-        val origin = structure.toEjml().origin
-
-        return when (type) {
-        ${
-        if (isDense)
-            """    InverseMatrixFeature::class -> object : InverseMatrixFeature<${type}> {
-                override val inverse: Matrix<${type}> by lazy {
-                    val res = origin.copy()
-                    CommonOps_${ops}.invert(res)
-                    res.wrapMatrix()
-                }
-            }
-
-            DeterminantFeature::class -> object : DeterminantFeature<${type}> {
-                override val determinant: $type by lazy { CommonOps_${ops}.det(origin) }
-            }
-
-            SingularValueDecompositionFeature::class -> object : SingularValueDecompositionFeature<${type}> {
-                private val svd by lazy {
-                    DecompositionFactory_${ops}.svd(origin.numRows, origin.numCols, true, true, false)
-                        .apply { decompose(origin.copy()) }
-                }
-
-                override val u: Matrix<${type}> by lazy { svd.getU(null, false).wrapMatrix() }
-                override val s: Matrix<${type}> by lazy { svd.getW(null).wrapMatrix() }
-                override val v: Matrix<${type}> by lazy { svd.getV(null, false).wrapMatrix() }
-                override val singularValues: Point<${type}> by lazy { ${type}Buffer(svd.singularValues) }
-            }
-
-            QRDecompositionFeature::class -> object : QRDecompositionFeature<${type}> {
-                private val qr by lazy {
-                    DecompositionFactory_${ops}.qr().apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<${type}> by lazy {
-                    qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature)
-                }
-
-                override val r: Matrix<${type}> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<${type}> {
-                override val l: Matrix<${type}> by lazy {
-                    val cholesky =
-                        DecompositionFactory_${ops}.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
-
-                    cholesky.getT(null).wrapMatrix().withFeature(LFeature)
-                }
-            }
-
-            LupDecompositionFeature::class -> object : LupDecompositionFeature<${type}> {
-                private val lup by lazy {
-                    DecompositionFactory_${ops}.lu(origin.numRows, origin.numCols).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<${type}> by lazy {
-                    lup.getLower(null).wrapMatrix().withFeature(LFeature)
-                }
-
-                override val u: Matrix<${type}> by lazy {
-                    lup.getUpper(null).wrapMatrix().withFeature(UFeature)
-                }
-
-                override val p: Matrix<${type}> by lazy { lup.getRowPivot(null).wrapMatrix() }
-            }""" else """    QRDecompositionFeature::class -> object : QRDecompositionFeature<$type> {
-                private val qr by lazy {
-                    DecompositionFactory_${ops}.qr(FillReducing.NONE).apply { decompose(origin.copy()) }
-                }
-
-                override val q: Matrix<${type}> by lazy {
-                    qr.getQ(null, false).wrapMatrix().withFeature(OrthogonalFeature)
-                }
-
-                override val r: Matrix<${type}> by lazy { qr.getR(null, false).wrapMatrix().withFeature(UFeature) }
-            }
-
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<${type}> {
-                override val l: Matrix<${type}> by lazy {
-                    val cholesky =
-                        DecompositionFactory_${ops}.cholesky().apply { decompose(origin.copy()) }
-
-                    (cholesky.getT(null) as ${ejmlMatrixParentTypeMatrix}).wrapMatrix().withFeature(LFeature)
-                }
-            }
-
-            LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
-                LUDecompositionFeature<${type}>, DeterminantFeature<${type}>, InverseMatrixFeature<${type}> {
-                private val lu by lazy {
-                    DecompositionFactory_${ops}.lu(FillReducing.NONE).apply { decompose(origin.copy()) }
-                }
-
-                override val l: Matrix<${type}> by lazy {
-                    lu.getLower(null).wrapMatrix().withFeature(LFeature)
-                }
-
-                override val u: Matrix<${type}> by lazy {
-                    lu.getUpper(null).wrapMatrix().withFeature(UFeature)
-                }
-
-                override val inverse: Matrix<${type}> by lazy {
-                    var a = origin
-                    val inverse = ${ejmlMatrixDenseType}(1, 1)
-                    val solver = LinearSolverFactory_${ops}.lu(FillReducing.NONE)
-                    if (solver.modifiesA()) a = a.copy()
-                    val i = CommonOps_${denseOps}.identity(a.numRows)
-                    solver.solve(i, inverse)
-                    inverse.wrapMatrix()
-                }
-
-                override val determinant: $type 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<${type}>, b: Matrix<${type}>): Ejml${type}Matrix<${ejmlMatrixType}> {
-        val res = ${ejmlMatrixType}(1, 1)
-        CommonOps_${ops}.solve(${ejmlMatrixType}(a.toEjml().origin), ${ejmlMatrixType}(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<${type}>, b: Point<${type}>): Ejml${type}Vector<${ejmlMatrixType}> {
-        val res = ${ejmlMatrixType}(1, 1)
-        CommonOps_${ops}.solve(${ejmlMatrixType}(a.toEjml().origin), ${ejmlMatrixType}(b.toEjml().origin), res)
-        return Ejml${type}Vector(res)
-    }
-}"""
-    appendLine(text)
-    appendLine()
-}
-
-
-/**
- * Generates routine EJML classes.
- */
-fun ejmlCodegen(outputFile: String): Unit = File(outputFile).run {
-    parentFile.mkdirs()
-
-    writer().use {
-        it.appendLine("/*")
-        it.appendLine(" * Copyright 2018-2024 KMath contributors.")
-        it.appendLine(" * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.")
-        it.appendLine(" */")
-        it.appendLine()
-        it.appendLine("/* This file is generated with buildSrc/src/main/kotlin/space/kscience/kmath/ejml/codegen/ejmlCodegen.kt */")
-        it.appendLine()
-        it.appendLine("package space.kscience.kmath.ejml")
-        it.appendLine()
-        it.appendLine("""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.attributes.SafeType
-import space.kscience.attributes.safeTypeOf
-import space.kscience.kmath.linear.*
-import space.kscience.kmath.linear.Matrix
-import space.kscience.kmath.UnstableKMathAPI
-import space.kscience.kmath.nd.StructureFeature
-import space.kscience.kmath.structures.Float64
-import space.kscience.kmath.structures.Float32
-import space.kscience.kmath.operations.Float64Field
-import space.kscience.kmath.operations.Float32Field
-import space.kscience.kmath.operations.DoubleField
-import space.kscience.kmath.operations.FloatField
-import space.kscience.kmath.operations.invoke
-import space.kscience.kmath.structures.Float64Buffer
-import space.kscience.kmath.structures.Float32Buffer
-import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.FloatBuffer
-import kotlin.reflect.KClass
-import kotlin.reflect.cast""")
-        it.appendLine()
-        it.appendEjmlVector("Double", "DMatrix")
-        it.appendEjmlVector("Float", "FMatrix")
-        it.appendEjmlMatrix("Double", "DMatrix")
-        it.appendEjmlMatrix("Float", "FMatrix")
-        it.appendEjmlLinearSpace("Double", "Float64Field", "DMatrix", "DMatrixRMaj", "DMatrixRMaj", "DDRM", "DDRM", true)
-        it.appendEjmlLinearSpace("Float", "Float32Field", "FMatrix", "FMatrixRMaj", "FMatrixRMaj", "FDRM", "FDRM", true)
-
-        it.appendEjmlLinearSpace(
-            type = "Double",
-            kmathAlgebra = "Float64Field",
-            ejmlMatrixParentTypeMatrix = "DMatrix",
-            ejmlMatrixType = "DMatrixSparseCSC",
-            ejmlMatrixDenseType = "DMatrixRMaj",
-            ops = "DSCC",
-            denseOps = "DDRM",
-            isDense = false,
-        )
-
-        it.appendEjmlLinearSpace(
-            type = "Float",
-            kmathAlgebra = "Float32Field",
-            ejmlMatrixParentTypeMatrix = "FMatrix",
-            ejmlMatrixType = "FMatrixSparseCSC",
-            ejmlMatrixDenseType = "FMatrixRMaj",
-            ops = "FSCC",
-            denseOps = "FDRM",
-            isDense = false,
-        )
-    }
-}

docs/templates/ARTIFACT-TEMPLATE.md

@@ -3,25 +3,10 @@
 The Maven coordinates of this project are `${group}:${name}:${version}`.
 
 **Gradle:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-    // development and snapshot versions
-    maven { url 'https://maven.pkg.jetbrains.space/spc/p/sci/dev' }
-}
-
-dependencies {
-    implementation '${group}:${name}:${version}'
-}
-```
-**Gradle Kotlin DSL:**
 ```kotlin
 repositories {
     maven("https://repo.kotlin.link")
     mavenCentral()
-    // development and snapshot versions
-    maven("https://maven.pkg.jetbrains.space/spc/p/sci/dev")
 }
 
 dependencies {

docs/templates/README-TEMPLATE.md

@@ -11,18 +11,21 @@ analog to Python's NumPy library. Later we found that kotlin is much more flexib
 architecture designs. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. The `numpy`-like
 experience could be achieved with [kmath-for-real](/kmath-for-real) extension module.
 
-[Documentation site (**WIP**)](https://SciProgCentre.github.io/kmath/)
+[Documentation site](https://SciProgCentre.github.io/kmath/)
 
 ## Publications and talks
 
 * [A conceptual article about context-oriented design](https://proandroiddev.com/an-introduction-context-oriented-programming-in-kotlin-2e79d316b0a2)
 * [Another article about context-oriented design](https://proandroiddev.com/diving-deeper-into-context-oriented-programming-in-kotlin-3ecb4ec38814)
 * [ACAT 2019 conference paper](https://aip.scitation.org/doi/abs/10.1063/1.5130103)
+* [A talk at KotlinConf 2019 about using kotlin for science](https://youtu.be/LI_5TZ7tnOE?si=4LknX41gl_YeUbIe)
+* [A talk on architecture at Joker-2021 (in Russian)](https://youtu.be/1bZ2doHiRRM?si=9w953ro9yu98X_KJ)
+* [The same talk in English](https://youtu.be/yP5DIc2fVwQ?si=louZzQ1dcXV6gP10)
+* [A seminar on tensor API](https://youtu.be/0H99wUs0xTM?si=6c__04jrByFQtVpo)
 
 # Goal
 
-* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM, JS and Native)
-  .
+* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM, JS, Native and Wasm).
 * Provide basic multiplatform implementations for those abstractions (without significant performance optimization).
 * Provide bindings and wrappers with those abstractions for popular optimized platform libraries.
 
@@ -73,7 +76,7 @@ native/SciPy (mostly due to boxing operations on primitive numbers). The best pe
 
 ## Requirements
 
-KMath currently relies on JDK 11 for compilation and execution of Kotlin-JVM part. We recommend using GraalVM-CE 11/17 for execution to get better performance.
+KMath currently relies on JDK 11 for compilation and execution of Kotlin-JVM part. We recommend using GraalVM-CE or Oracle GraalVM for execution to get better performance.
 
 ### Repositories
 
@@ -95,7 +98,5 @@ dependencies {
 
 ## Contributing
 
-The project requires a lot of additional work. The most important thing we need is a feedback about what features are
-required the most. Feel free to create feature requests. We are also welcome to code contributions, especially in issues
-marked with
-[good first issue](hhttps://github.com/SciProgCentre/kmath/issues?q=is%3Aissue+is%3Aopen+label%3A%22good+first+issue%22) label.
+The project requires a lot of additional work. The most important thing we need is feedback about what features are
+required the most. Feel free to create feature requests. We are also welcome to code contributions, especially in issues marked with [good first issue](hhttps://github.com/SciProgCentre/kmath/issues?q=is%3Aissue+is%3Aopen+label%3A%22good+first+issue%22) label.

examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt

@@ -13,10 +13,7 @@ import space.kscience.kmath.expressions.autodiff
 import space.kscience.kmath.expressions.symbol
 import space.kscience.kmath.operations.asIterable
 import space.kscience.kmath.operations.toList
-import space.kscience.kmath.optimization.FunctionOptimizationTarget
-import space.kscience.kmath.optimization.optimizeWith
-import space.kscience.kmath.optimization.result
-import space.kscience.kmath.optimization.resultValue
+import space.kscience.kmath.optimization.*
 import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.real.DoubleVector
 import space.kscience.kmath.real.map
@@ -80,8 +77,9 @@ suspend fun main() {
     val result = chi2.optimizeWith(
         CMOptimizer,
         mapOf(a to 1.5, b to 0.9, c to 1.0),
-        FunctionOptimizationTarget.MINIMIZE
-    )
+    ){
+        FunctionOptimizationTarget(OptimizationDirection.MINIMIZE)
+    }
 
     //display a page with plot and numerical results
     val page = Plotly.page {

examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt

@@ -7,6 +7,7 @@ package space.kscience.kmath.fit
 
 import kotlinx.html.br
 import kotlinx.html.h3
+import space.kscience.attributes.Attributes
 import space.kscience.kmath.data.XYErrorColumnarData
 import space.kscience.kmath.distributions.NormalDistribution
 import space.kscience.kmath.expressions.Symbol
@@ -64,7 +65,7 @@ suspend fun main() {
         QowOptimizer,
         Double.autodiff,
         mapOf(a to 0.9, b to 1.2, c to 2.0, e to 1.0, d to 1.0, e to 0.0),
-        OptimizationParameters(a, b, c, d)
+        attributes = Attributes(OptimizationParameters, listOf(a, b, c, d))
     ) { arg ->
         //bind variables to autodiff context
         val a by binding

examples/src/main/kotlin/space/kscience/kmath/functions/interpolate.kt

@@ -8,7 +8,6 @@ package space.kscience.kmath.functions
 import space.kscience.kmath.interpolation.SplineInterpolator
 import space.kscience.kmath.interpolation.interpolatePolynomials
 import space.kscience.kmath.operations.Float64Field
-import space.kscience.kmath.structures.Float64Buffer
 import space.kscience.plotly.Plotly
 import space.kscience.plotly.UnstablePlotlyAPI
 import space.kscience.plotly.makeFile
@@ -24,9 +23,7 @@ fun main() {
         x to sin(x)
     }
 
-    val polynomial: PiecewisePolynomial<Double> = SplineInterpolator(
-        Float64Field, ::Float64Buffer
-    ).interpolatePolynomials(data)
+    val polynomial: PiecewisePolynomial<Double> = SplineInterpolator(Float64Field).interpolatePolynomials(data)
 
     val function = polynomial.asFunction(Float64Field, 0.0)
 

examples/src/main/kotlin/space/kscience/kmath/linear/dotPerformance.kt

@@ -5,29 +5,24 @@
 
 package space.kscience.kmath.linear
 
-import space.kscience.kmath.operations.algebra
 import kotlin.random.Random
-import kotlin.time.ExperimentalTime
 import kotlin.time.measureTime
 
-@OptIn(ExperimentalTime::class)
-fun main() {
+fun main() = with(Float64ParallelLinearSpace) {
     val random = Random(12224)
     val dim = 1000
 
     //creating invertible matrix
-    val matrix1 = Double.algebra.linearSpace.buildMatrix(dim, dim) { i, j ->
+    val matrix1 = buildMatrix(dim, dim) { i, j ->
         if (i <= j) random.nextDouble() else 0.0
     }
-    val matrix2 = Double.algebra.linearSpace.buildMatrix(dim, dim) { i, j ->
+    val matrix2 = buildMatrix(dim, dim) { i, j ->
         if (i <= j) random.nextDouble() else 0.0
     }
 
     val time = measureTime {
-        with(Double.algebra.linearSpace) {
-            repeat(10) {
-                matrix1 dot matrix2
-            }
+        repeat(30) {
+            matrix1 dot matrix2
         }
     }
 

examples/src/main/kotlin/space/kscience/kmath/linear/lupPerformance.kt

@@ -0,0 +1,27 @@
+/*
+ * Copyright 2018-2024 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.kmath.linear
+
+import kotlin.random.Random
+import kotlin.time.measureTime
+
+fun main(): Unit = with(Float64LinearSpace) {
+    val random = Random(1224)
+    val dim = 500
+
+    //creating invertible matrix
+    val u = buildMatrix(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+    val l = buildMatrix(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 }
+    val matrix = l dot u
+
+    val time = measureTime {
+        repeat(20) {
+            lupSolver().inverse(matrix)
+        }
+    }
+
+    println(time)
+}

examples/src/main/kotlin/space/kscience/kmath/series/seriesBuilder.kt

@@ -7,7 +7,7 @@ package space.kscience.kmath.series
 
 
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.Float64Buffer
 import space.kscience.kmath.structures.asBuffer
 import space.kscience.kmath.structures.toDoubleArray
 import space.kscience.plotly.*
@@ -21,7 +21,7 @@ fun main(): Unit = with(Double.seriesAlgebra()) {
 
     val arrayOfRandoms = DoubleArray(20) { random.nextDouble() }
 
-    val series1: DoubleBuffer = arrayOfRandoms.asBuffer()
+    val series1: Float64Buffer = arrayOfRandoms.asBuffer()
     val series2: Series<Double> = series1.moveBy(3)
 
     val res = series2 - series1
@@ -42,7 +42,7 @@ fun main(): Unit = with(Double.seriesAlgebra()) {
     Plotly.plot {
         series("series1", series1)
         series("series2", series2)
-        series("dif", res){
+        series("dif", res) {
             mode = ScatterMode.lines
             line.color("magenta")
         }

examples/src/main/kotlin/space/kscience/kmath/structures/ComplexND.kt

@@ -8,7 +8,7 @@
 package space.kscience.kmath.structures
 
 import space.kscience.kmath.complex.*
-import space.kscience.kmath.linear.transpose
+import space.kscience.kmath.linear.transposed
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.nd.as2D
 import space.kscience.kmath.nd.ndAlgebra
@@ -60,7 +60,7 @@ fun complexExample() {
             val sum = matrix + x + 1.0
 
             //Represent the sum as 2d-structure and transpose
-            sum.as2D().transpose()
+            sum.as2D().transposed()
         }
     }
 }

examples/src/main/kotlin/space/kscience/kmath/structures/StreamDoubleFieldND.kt

@@ -7,6 +7,7 @@ package space.kscience.kmath.structures
 
 import space.kscience.kmath.PerformancePitfall
 import space.kscience.kmath.nd.*
+import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.ExtendedField
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.NumbersAddOps

examples/src/main/kotlin/space/kscience/kmath/structures/StructureWriteBenchmark.kt

@@ -6,20 +6,18 @@
 package space.kscience.kmath.structures
 
 import space.kscience.kmath.nd.BufferND
-import space.kscience.kmath.nd.ShapeND
-import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.operations.mapToBuffer
 import kotlin.system.measureTimeMillis
 
 private inline fun <T, reified R : Any> BufferND<T>.mapToBufferND(
-    bufferFactory: BufferFactory<R> = BufferFactory.auto(),
+    bufferFactory: BufferFactory<R> = BufferFactory(),
     crossinline block: (T) -> R,
 ): BufferND<R> = BufferND(indices, buffer.mapToBuffer(bufferFactory, block))
 
 @Suppress("UNUSED_VARIABLE")
 fun main() {
     val n = 6000
-    val structure = StructureND.buffered(ShapeND(n, n), Buffer.Companion::auto) { 1.0 }
+    val structure = BufferND(n, n) { 1.0 }
     structure.mapToBufferND { it + 1 } // warm-up
     val time1 = measureTimeMillis { val res = structure.mapToBufferND { it + 1 } }
     println("Structure mapping finished in $time1 millis")

examples/src/main/kotlin/space/kscience/kmath/structures/buffers.kt

@@ -13,7 +13,7 @@ import space.kscience.kmath.operations.withSize
 inline fun <reified R : Any> MutableBuffer.Companion.same(
     n: Int,
     value: R
-): MutableBuffer<R> = auto(n) { value }
+): MutableBuffer<R> = MutableBuffer(n) { value }
 
 
 fun main() {

gradle/wrapper/gradle-wrapper.properties

@@ -1,5 +1,5 @@
 distributionBase=GRADLE_USER_HOME
 distributionPath=wrapper/dists
-distributionUrl=https\://services.gradle.org/distributions/gradle-8.4-bin.zip
+distributionUrl=https\://services.gradle.org/distributions/gradle-8.6-bin.zip
 zipStoreBase=GRADLE_USER_HOME
 zipStorePath=wrapper/dists

kmath-ast/README.md

@@ -10,19 +10,8 @@ Extensions to MST API: transformations, dynamic compilation and visualization.
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-ast:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-ast:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-ast:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -31,7 +20,7 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-ast:0.4.0-dev-1")
+    implementation("space.kscience:kmath-ast:0.4.0-dev-3")
 }
 ```
 

kmath-commons/README.md

@@ -6,19 +6,8 @@ Commons math binding for kmath
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-commons:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-commons:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-commons:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -27,6 +16,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-commons:0.4.0-dev-1")
+    implementation("space.kscience:kmath-commons:0.4.0-dev-3")
 }
 ```

kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt

@@ -7,16 +7,20 @@ package space.kscience.kmath.commons.linear
 
 import org.apache.commons.math3.linear.*
 import org.apache.commons.math3.linear.LUDecomposition
+import org.apache.commons.math3.linear.SingularValueDecomposition
 import space.kscience.attributes.SafeType
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.linear.*
+import space.kscience.kmath.linear.CholeskyDecomposition
+import space.kscience.kmath.linear.QRDecomposition
 import space.kscience.kmath.nd.Structure2D
 import space.kscience.kmath.nd.StructureAttribute
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.Float64Buffer
-import kotlin.reflect.cast
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.structures.IntBuffer
+import space.kscience.kmath.structures.asBuffer
 
 public class CMMatrix(public val origin: RealMatrix) : Matrix<Double> {
     override val type: SafeType<Double> get() = DoubleField.type
@@ -109,45 +113,44 @@ public object CMLinearSpace : LinearSpace<Double, Float64Field> {
 
         val origin = structure.toCM().origin
 
-        return when (attribute) {
-            IsDiagonal -> if (origin is DiagonalMatrix) IsDiagonal else null
+        val raw: Any? = when (attribute) {
+            IsDiagonal -> if (origin is DiagonalMatrix) Unit else null
             Determinant -> LUDecomposition(origin).determinant
-            LUP -> GenericLupDecomposition {
-                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).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) }
+
+            LUP -> object : LupDecomposition<Float64> {
+                val lup by lazy { LUDecomposition(origin) }
+                override val pivot: IntBuffer get() = lup.pivot.asBuffer()
+                override val l: Matrix<Float64> get() = lup.l.wrap()
+                override val u: Matrix<Float64> get() = lup.u.wrap()
             }
 
-            CholeskyDecompositionAttribute -> object : CholeskyDecompositionAttribute<Double> {
-                override val l: Matrix<Double> by lazy<Matrix<Double>> {
-                    val cholesky = CholeskyDecomposition(origin)
-                    CMMatrix(cholesky.l).withAttribute(LowerTriangular)
-                }
+            Cholesky -> object : CholeskyDecomposition<Float64> {
+                val cmCholesky by lazy { org.apache.commons.math3.linear.CholeskyDecomposition(origin) }
+                override val l: Matrix<Double> get() = cmCholesky.l.wrap()
             }
 
-            QRDecompositionAttribute -> object : QRDecompositionAttribute<Double> {
-                private val qr by lazy { QRDecomposition(origin) }
-                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) }
+            QR -> object : QRDecomposition<Float64> {
+                val cmQr by lazy { org.apache.commons.math3.linear.QRDecomposition(origin) }
+                override val q: Matrix<Float64> get() = cmQr.q.wrap().withAttribute(OrthogonalAttribute)
+                override val r: Matrix<Float64> get() = cmQr.r.wrap().withAttribute(UpperTriangular)
             }
 
-            SVDAttribute -> object : SVDAttribute<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 { Float64Buffer(sv.singularValues) }
+            SVD -> object : space.kscience.kmath.linear.SingularValueDecomposition<Float64> {
+                val cmSvd by lazy { SingularValueDecomposition(origin) }
+
+                override val u: Matrix<Float64> get() = cmSvd.u.wrap()
+                override val s: Matrix<Float64> get() = cmSvd.s.wrap()
+                override val v: Matrix<Float64> get() = cmSvd.v.wrap()
+                override val singularValues: Point<Float64> get() = cmSvd.singularValues.asBuffer()
+
             }
 
             else -> null
-        }?.let(type::cast)
+        }
+        @Suppress("UNCHECKED_CAST")
+        return raw as V?
     }
+
 }
 
 public operator fun CMMatrix.plus(other: CMMatrix): CMMatrix = CMMatrix(origin.add(other.origin))

kmath-commons/src/jvmMain/kotlin/space/kscience/kmath/commons/optimization/CMOptimizer.kt

@@ -3,6 +3,7 @@
  * 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)
+
 package space.kscience.kmath.commons.optimization
 
 import org.apache.commons.math3.optim.*
@@ -20,7 +21,6 @@ import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.expressions.SymbolIndexer
 import space.kscience.kmath.expressions.derivative
 import space.kscience.kmath.expressions.withSymbols
-import space.kscience.kmath.misc.log
 import space.kscience.kmath.optimization.*
 import kotlin.collections.set
 import kotlin.reflect.KClass
@@ -28,7 +28,7 @@ import kotlin.reflect.KClass
 public operator fun PointValuePair.component1(): DoubleArray = point
 public operator fun PointValuePair.component2(): Double = value
 
-public object CMOptimizerEngine: OptimizationAttribute<() -> MultivariateOptimizer>
+public object CMOptimizerEngine : OptimizationAttribute<() -> MultivariateOptimizer>
 
 /**
  * Specify a Commons-maths optimization engine
@@ -37,7 +37,7 @@ public fun AttributesBuilder<FunctionOptimization<Double>>.cmEngine(optimizerBui
     set(CMOptimizerEngine, optimizerBuilder)
 }
 
-public object CMOptimizerData: SetAttribute<SymbolIndexer.() -> OptimizationData>
+public object CMOptimizerData : SetAttribute<SymbolIndexer.() -> OptimizationData>
 
 /**
  * Specify Commons-maths optimization data.
@@ -118,21 +118,24 @@ public object CMOptimizer : Optimizer<Double, FunctionOptimization<Double>> {
 
             val logger = problem.attributes[OptimizationLog]
 
-            for (feature in problem.attributes) {
-                when (feature) {
-                    is CMOptimizerData -> feature.data.forEach { dataBuilder ->
-                        addOptimizationData(dataBuilder())
-                    }
-                    is FunctionOptimizationTarget -> when (feature) {
-                        FunctionOptimizationTarget.MAXIMIZE -> addOptimizationData(GoalType.MAXIMIZE)
-                        FunctionOptimizationTarget.MINIMIZE -> addOptimizationData(GoalType.MINIMIZE)
-                    }
-                    else -> logger?.log { "The feature $feature is unused in optimization" }
+            problem.attributes[CMOptimizerData]?.let { builders: Set<SymbolIndexer.() -> OptimizationData> ->
+                builders.forEach { dataBuilder ->
+                    addOptimizationData(dataBuilder())
+                }
+            }
+
+            problem.attributes[FunctionOptimizationTarget]?.let { direction: OptimizationDirection ->
+                when (direction) {
+                    OptimizationDirection.MAXIMIZE -> addOptimizationData(GoalType.MAXIMIZE)
+                    OptimizationDirection.MINIMIZE -> addOptimizationData(GoalType.MINIMIZE)
                 }
             }
 
             val (point, value) = cmOptimizer.optimize(*optimizationData.values.toTypedArray())
-            return problem.withAttributes(OptimizationResult(point.toMap()), OptimizationValue(value))
+            return problem.withAttributes {
+                result(point.toMap())
+                value(value)
+            }
         }
     }
 }

kmath-commons/src/jvmTest/kotlin/space/kscience/kmath/commons/integration/IntegrationTest.kt

@@ -7,6 +7,8 @@ package space.kscience.kmath.commons.integration
 
 import org.junit.jupiter.api.Test
 import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.integration.IntegrandAbsoluteAccuracy
+import space.kscience.kmath.integration.IntegrandRelativeAccuracy
 import space.kscience.kmath.integration.integrate
 import space.kscience.kmath.integration.value
 import space.kscience.kmath.operations.Float64Field.sin
@@ -27,8 +29,8 @@ internal class IntegrationTest {
     @Test
     fun customSimpson() {
         val res = CMIntegrator.simpson().integrate(0.0..PI, {
-            targetRelativeAccuracy = 1e-4
-            targetAbsoluteAccuracy = 1e-4
+            IntegrandRelativeAccuracy(1e-4)
+            IntegrandAbsoluteAccuracy(1e-4)
         }, function).value
         assertTrue { abs(res - 2) < 1e-3 }
         assertTrue { abs(res - 2) > 1e-12 }

kmath-commons/src/jvmTest/kotlin/space/kscience/kmath/commons/optimization/OptimizeTest.kt

@@ -73,8 +73,9 @@ internal class OptimizeTest {
         val result: FunctionOptimization<Double> = chi2.optimizeWith(
             CMOptimizer,
             mapOf(a to 1.5, b to 0.9, c to 1.0),
-            FunctionOptimizationTarget.MINIMIZE
-        )
+        ){
+            FunctionOptimizationTarget(OptimizationDirection.MINIMIZE)
+        }
         println(result)
         println("Chi2/dof = ${result.resultValue / (x.size - 3)}")
     }

kmath-complex/README.md

@@ -8,19 +8,8 @@ Complex and hypercomplex number systems in KMath.
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-complex:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-complex:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-complex:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -29,6 +18,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-complex:0.4.0-dev-1")
+    implementation("space.kscience:kmath-complex:0.4.0-dev-3")
 }
 ```

kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/Complex.kt

@@ -144,12 +144,19 @@ public object ComplexField :
                 zero
             }
         }
-
     } else {
         exp(pow * ln(arg))
     }
 
-    public fun power(arg: Complex, pow: Complex): Complex = exp(pow * ln(arg))
+    public fun power(arg: Complex, pow: Complex): Complex = if(arg == zero || arg == (-0.0).toComplex()){
+        if(pow == zero){
+            one
+        } else {
+            zero
+        }
+    } else {
+        exp(pow * ln(arg))
+    }
 
     public fun Complex.pow(power: Complex): Complex = power(this, power)
 

kmath-complex/src/commonTest/kotlin/space/kscience/kmath/complex/ComplexFieldTest.kt

@@ -71,6 +71,11 @@ internal class ComplexFieldTest {
             (i * 8).let { it.im.toInt() to it.re.toInt() },
             (Complex(2, 2) pow 2).let { it.im.toInt() to it.re.toInt() })
 
+        assertEquals(1.0, Complex(0.0).pow(Complex(0.0)).re, 0.01)
+        assertEquals(1.0, Complex(-0.0).pow(Complex(0.0)).re, 0.01)
+        assertEquals(0.0, Complex(0.0).pow(Complex(2.0)).re, 0.01)
+        assertEquals(0.0, Complex(-0.0).pow(Complex(2.0)).re, 0.01)
+        assertEquals(1.0, Complex(-1.0).pow(Complex(2.0)).re, 0.01)
         assertEquals(2.0, power(Complex(-4.0, 0.0), 0.5 + 0 * i).im, 0.01)
         assertEquals(2.0, power(Complex(-4.0, 0.0), 0.5).im, 0.01)
     }

kmath-core/README.md

@@ -11,23 +11,13 @@ objects to the expression by providing a context. Expressions can be used for a
 performance calculations to code generation.
  - [domains](src/commonMain/kotlin/space/kscience/kmath/domains) : Domains
  - [autodiff](src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt) : Automatic differentiation
+ - [linear.parallel](#) : Parallel implementation for `LinearAlgebra`
 
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-core:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-core:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-core:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -36,6 +26,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-core:0.4.0-dev-1")
+    implementation("space.kscience:kmath-core:0.4.0-dev-3")
 }
 ```

kmath-core/build.gradle.kts

@@ -71,4 +71,12 @@ readme {
         id = "autodiff",
         ref = "src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt"
     ) { "Automatic differentiation" }
+
+    feature(
+        id="Parallel linear algebra"
+    ){
+        """
+            Parallel implementation for `LinearAlgebra`
+        """.trimIndent()
+    }
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/Expression.kt

@@ -7,6 +7,7 @@ package space.kscience.kmath.expressions
 
 import space.kscience.attributes.SafeType
 import space.kscience.attributes.WithType
+import space.kscience.attributes.safeTypeOf
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.operations.Algebra
 import space.kscience.kmath.operations.DoubleField
@@ -30,12 +31,18 @@ public interface Expression<T> : WithType<T> {
     public operator fun invoke(arguments: Map<Symbol, T>): T
 }
 
+/**
+ * Create an expression from a functional block.
+ */
 public fun <T> Expression(type: SafeType<T>, block: (Map<Symbol, T>) -> T): Expression<T> = object : Expression<T> {
     override fun invoke(arguments: Map<Symbol, T>): T = block(arguments)
 
     override val type: SafeType<T> = type
 }
 
+public inline fun <reified T> Expression(noinline block: (Map<Symbol, T>) -> T): Expression<T> =
+    Expression(safeTypeOf<T>(), block)
+
 /**
  * Specialization of [Expression] for [Double] allowing better performance because of using array.
  */

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/Float64LinearSpace.kt

@@ -54,11 +54,12 @@ public object Float64LinearSpace : LinearSpace<Double, Float64Field> {
         require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
         val rows = this@dot.rows.map { it.linearize() }
         val columns = other.columns.map { it.linearize() }
+        val indices = 0 until this.colNum
         return buildMatrix(rowNum, other.colNum) { i, j ->
             val r = rows[i]
             val c = columns[j]
             var res = 0.0
-            for (l in r.indices) {
+            for (l in indices) {
                 res += r[l] * c[l]
             }
             res
@@ -69,10 +70,11 @@ public object Float64LinearSpace : LinearSpace<Double, Float64Field> {
     override fun Matrix<Double>.dot(vector: Point<Double>): Float64Buffer {
         require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
         val rows = this@dot.rows.map { it.linearize() }
+        val indices = 0 until this.colNum
         return Float64Buffer(rowNum) { i ->
             val r = rows[i]
             var res = 0.0
-            for (j in r.indices) {
+            for (j in indices) {
                 res += r[j] * vector[j]
             }
             res

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LupDecomposition.kt

@@ -4,25 +4,33 @@
  */
 
 @file:Suppress("UnusedReceiverParameter")
+@file:OptIn(PerformancePitfall::class)
 
 package space.kscience.kmath.linear
 
 import space.kscience.attributes.Attributes
 import space.kscience.attributes.PolymorphicAttribute
 import space.kscience.attributes.safeTypeOf
+import space.kscience.kmath.PerformancePitfall
 import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.*
 import space.kscience.kmath.structures.*
 
 public interface LupDecomposition<T> {
-    public val linearSpace: LinearSpace<T, Field<T>>
-    public val elementAlgebra: Field<T> get() = linearSpace.elementAlgebra
-
     public val pivot: IntBuffer
     public val l: Matrix<T>
     public val u: Matrix<T>
 }
 
+/**
+ * Create a pivot matrix from pivot vector using provided [LinearSpace]
+ */
+public fun <T> LupDecomposition<T>.pivotMatrix(linearSpace: LinearSpace<T, Ring<T>>): Matrix<T> =
+    VirtualMatrix(linearSpace.type, l.rowNum, l.colNum) { row, column ->
+        if (column == pivot[row]) linearSpace.elementAlgebra.one else linearSpace.elementAlgebra.zero
+    }
+
 /**
  * Matrices with this feature support LU factorization with partial pivoting: *[p] &middot; a = [l] &middot; [u]* where
  * *a* is the owning matrix.
@@ -31,12 +39,13 @@ public interface LupDecomposition<T> {
  * @param lu combined L and U matrix
  */
 public class GenericLupDecomposition<T>(
-    override val linearSpace: LinearSpace<T, Field<T>>,
+    public val elementAlgebra: Field<T>,
     private val lu: Matrix<T>,
     override val pivot: IntBuffer,
     private val even: Boolean,
 ) : LupDecomposition<T> {
 
+
     override val l: Matrix<T>
         get() = VirtualMatrix(lu.type, lu.rowNum, lu.colNum, attributes = Attributes(LowerTriangular)) { i, j ->
             when {
@@ -51,11 +60,6 @@ public class GenericLupDecomposition<T>(
             if (j >= i) lu[i, j] else elementAlgebra.zero
         }
 
-    public val pivotMatrix: Matrix<T>
-        get() = VirtualMatrix(linearSpace.type, l.rowNum, l.colNum) { row, column ->
-            if (column == pivot[row]) elementAlgebra.one else elementAlgebra.zero
-        }
-
     public val determinant: T by lazy {
         elementAlgebra { (0 until l.shape[0]).fold(if (even) one else -one) { value, i -> value * lu[i, i] } }
     }
@@ -76,136 +80,145 @@ 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(
+public fun <T : Comparable<T>> Field<T>.lup(
     matrix: Matrix<T>,
     checkSingular: (T) -> Boolean,
-): LupDecomposition<T> = elementAlgebra {
+): GenericLupDecomposition<T>  {
     require(matrix.rowNum == matrix.colNum) { "LU decomposition supports only square matrices" }
     val m = matrix.colNum
     val pivot = IntArray(matrix.rowNum)
 
-    //TODO just waits for multi-receivers
-    with(BufferAccessor2D(matrix.rowNum, matrix.colNum, elementAlgebra.bufferFactory)) {
+    val strides = RowStrides(ShapeND(matrix.rowNum, matrix.colNum))
 
-        val lu = create(matrix)
+    val lu = MutableBufferND(
+        strides,
+        bufferAlgebra.buffer(strides.linearSize) { offset ->
+            matrix[strides.index(offset)]
+        }
+    )
 
-        // Initialize the permutation array and parity
-        for (row in 0 until m) pivot[row] = row
-        var even = true
 
-        // Initialize the permutation array and parity
-        for (row in 0 until m) pivot[row] = row
+    // Initialize the permutation array and parity
+    for (row in 0 until m) pivot[row] = row
+    var even = true
 
-        // Loop over columns
-        for (col in 0 until m) {
-            // upper
-            for (row in 0 until col) {
-                val luRow = lu.row(row)
-                var sum = luRow[col]
-                for (i in 0 until row) sum -= luRow[i] * lu[i, col]
-                luRow[col] = sum
+    // Initialize the permutation array and parity
+    for (row in 0 until m) pivot[row] = row
+
+    // Loop over columns
+    for (col in 0 until m) {
+        // upper
+        for (row in 0 until col) {
+            var sum = lu[row, col]
+            for (i in 0 until row){
+                sum -= lu[row, i] * lu[i, col]
             }
-
-            // lower
-            var max = col // permutation row
-            var largest = -one
-
-            for (row in col until m) {
-                val luRow = lu.row(row)
-                var sum = luRow[col]
-                for (i in 0 until col) sum -= luRow[i] * lu[i, col]
-                luRow[col] = sum
-
-                // maintain the best permutation choice
-                if (abs(sum) > largest) {
-                    largest = abs(sum)
-                    max = row
-                }
-            }
-
-            // Singularity check
-            check(!checkSingular(abs(lu[max, col]))) { "The matrix is singular" }
-
-            // Pivot if necessary
-            if (max != col) {
-                val luMax = lu.row(max)
-                val luCol = lu.row(col)
-
-                for (i in 0 until m) {
-                    val tmp = luMax[i]
-                    luMax[i] = luCol[i]
-                    luCol[i] = tmp
-                }
-
-                val temp = pivot[max]
-                pivot[max] = pivot[col]
-                pivot[col] = temp
-                even = !even
-            }
-
-            // Divide the lower elements by the "winning" diagonal elt.
-            val luDiag = lu[col, col]
-            for (row in col + 1 until m) lu[row, col] /= luDiag
+            lu[row, col] = sum
         }
 
-        return GenericLupDecomposition(this@lup, lu.toStructure2D(), pivot.asBuffer(), even)
+        // lower
+        var max = col // permutation row
+        var largest = -one
+
+        for (row in col until m) {
+            var sum = lu[row, col]
+            for (i in 0 until col){
+                sum -= lu[row, i] * lu[i, col]
+            }
+            lu[row, col] = sum
+
+            // maintain the best permutation choice
+            if (abs(sum) > largest) {
+                largest = abs(sum)
+                max = row
+            }
+        }
+
+        // Singularity check
+        check(!checkSingular(abs(lu[max, col]))) { "The matrix is singular" }
+
+        // Pivot if necessary
+        if (max != col) {
+            for (i in 0 until m) {
+                val tmp = lu[max, i]
+                lu[max, i] = lu[col, i]
+                lu[col, i] = tmp
+            }
+
+            val temp = pivot[max]
+            pivot[max] = pivot[col]
+            pivot[col] = temp
+            even = !even
+        }
+
+        // Divide the lower elements by the "winning" diagonal elt.
+        val luDiag = lu[col, col]
+        for (row in col + 1 until m) lu[row, col] /= luDiag
     }
+
+
+    return GenericLupDecomposition(this, lu.as2D(), pivot.asBuffer(), even)
+
 }
 
 
-public fun LinearSpace<Double, Float64Field>.lup(
+public fun Field<Float64>.lup(
     matrix: Matrix<Double>,
     singularityThreshold: Double = 1e-11,
-): LupDecomposition<Double> = lup(matrix) { it < singularityThreshold }
+): GenericLupDecomposition<Double> = lup(matrix) { it < singularityThreshold }
 
-internal fun <T> LinearSpace<T, Field<T>>.solve(
+private fun <T> Field<T>.solve(
     lup: LupDecomposition<T>,
     matrix: Matrix<T>,
 ): Matrix<T> {
     require(matrix.rowNum == lup.l.rowNum) { "Matrix dimension mismatch. Expected ${lup.l.rowNum}, but got ${matrix.colNum}" }
 
-    BufferAccessor2D(matrix.rowNum, matrix.colNum, elementAlgebra.bufferFactory).run {
-        elementAlgebra {
-            // Apply permutations to b
-            val bp = create { _, _ -> zero }
+    val strides = RowStrides(ShapeND(matrix.rowNum, matrix.colNum))
 
-            for (row in 0 until rowNum) {
-                val bpRow = bp.row(row)
-                val pRow = lup.pivot[row]
-                for (col in 0 until matrix.colNum) bpRow[col] = matrix[pRow, col]
-            }
+    // Apply permutations to b
+    val bp = MutableBufferND(
+        strides,
+        bufferAlgebra.buffer(strides.linearSize) { _ -> zero }
+    )
 
-            // Solve LY = b
-            for (col in 0 until colNum) {
-                val bpCol = bp.row(col)
 
-                for (i in col + 1 until colNum) {
-                    val bpI = bp.row(i)
-                    val luICol = lup.l[i, col]
-                    for (j in 0 until matrix.colNum) {
-                        bpI[j] -= bpCol[j] * luICol
-                    }
-                }
-            }
-
-            // Solve UX = Y
-            for (col in colNum - 1 downTo 0) {
-                val bpCol = bp.row(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 = lup.u[i, col]
-                    for (j in 0 until matrix.colNum) bpI[j] -= bpCol[j] * luICol
-                }
-            }
-
-            return buildMatrix(matrix.rowNum, matrix.colNum) { i, j -> bp[i, j] }
+    for (row in 0 until matrix.rowNum) {
+        val pRow = lup.pivot[row]
+        for (col in 0 until matrix.colNum) {
+            bp[row, col] = matrix[pRow, col]
         }
     }
+
+    // Solve LY = b
+    for (col in 0 until matrix.colNum) {
+
+        for (i in col + 1 until matrix.colNum) {
+            val luICol = lup.l[i, col]
+            for (j in 0 until matrix.colNum) {
+                bp[i, j] -= bp[col, j] * luICol
+            }
+        }
+    }
+
+    // Solve UX = Y
+    for (col in matrix.colNum - 1 downTo 0) {
+        val luDiag = lup.u[col, col]
+        for (j in 0 until matrix.colNum) {
+            bp[col, j] /= luDiag
+        }
+
+        for (i in 0 until col) {
+            val luICol = lup.u[i, col]
+            for (j in 0 until matrix.colNum) {
+                bp[i, j] -= bp[col, j] * luICol
+            }
+        }
+    }
+
+    return bp.as2D()
 }
 
+
 /**
  * Produce a generic solver based on LUP decomposition
  */
@@ -213,7 +226,7 @@ internal fun <T> LinearSpace<T, Field<T>>.solve(
 public fun <T : Comparable<T>> LinearSpace<T, Field<T>>.lupSolver(
     singularityCheck: (T) -> Boolean,
 ): LinearSolver<T> = object : LinearSolver<T> {
-    override fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
+    override fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T>  = elementAlgebra{
         // Use existing decomposition if it is provided by matrix or linear space itself
         val decomposition = a.getOrComputeAttribute(LUP) ?: lup(a, singularityCheck)
         return solve(decomposition, b)

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/Transposed.kt

@@ -25,5 +25,4 @@ public class TransposedMatrix<T>(public val origin: Matrix<T>) : Matrix<T> {
 /**
  * Create a virtual transposed matrix without copying anything. `A.transpose().transpose() === A`
  */
-public val <T> Matrix<T>.transposed: Matrix<T>
-    get() = (this as? TransposedMatrix<T>)?.origin ?: TransposedMatrix(this)
+public fun <T> Matrix<T>.transposed(): Matrix<T> = (this as? TransposedMatrix<T>)?.origin ?: TransposedMatrix(this)

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

@@ -5,9 +5,11 @@
 
 package space.kscience.kmath.nd
 
+import space.kscience.attributes.SafeType
 import space.kscience.kmath.PerformancePitfall
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.Float64
 import space.kscience.kmath.structures.Float64Buffer
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract
@@ -20,7 +22,10 @@ import kotlin.math.pow as kpow
 public class Float64BufferND(
     indexes: ShapeIndexer,
     override val buffer: Float64Buffer,
-) : MutableBufferND<Double>(indexes, buffer), MutableStructureNDOfDouble{
+) : MutableBufferND<Double>(indexes, buffer), MutableStructureNDOfDouble {
+
+    override val type: SafeType<Float64> get() = Float64Field.type
+
     override fun getDouble(index: IntArray): Double = buffer[indices.offset(index)]
 
     override fun setDouble(index: IntArray, value: Double) {

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

@@ -69,6 +69,29 @@ public interface Structure2D<out T> : StructureND<T> {
     public companion object
 }
 
+/**
+ * A linear accessor for a [MutableStructureND]
+ */
+@OptIn(PerformancePitfall::class)
+public class MutableStructureNDAccessorBuffer<T>(
+    public val structure: MutableStructureND<T>,
+    override val size: Int,
+    private val indexer: (Int) -> IntArray,
+) : MutableBuffer<T> {
+
+    override val type: SafeType<T> get() = structure.type
+
+    override fun set(index: Int, value: T) {
+        structure[indexer(index)] = value
+    }
+
+    override fun get(index: Int): T = structure[indexer(index)]
+
+    override fun toString(): String = "AccessorBuffer(structure=$structure, size=$size)"
+
+    override fun copy(): MutableBuffer<T> = MutableBuffer(type, size, ::get)
+}
+
 /**
  * Represents mutable [Structure2D].
  */
@@ -87,14 +110,18 @@ public interface MutableStructure2D<T> : Structure2D<T>, MutableStructureND<T> {
      */
     @PerformancePitfall
     override val rows: List<MutableBuffer<T>>
-        get() = List(rowNum) { i -> MutableBuffer(type, colNum) { j -> get(i, j) } }
+        get() = List(rowNum) { i ->
+            MutableStructureNDAccessorBuffer(this, colNum) { j -> intArrayOf(i, j) }
+        }
 
     /**
      * The buffer of columns for this structure. It gets elements from the structure dynamically.
      */
     @PerformancePitfall
     override val columns: List<MutableBuffer<T>>
-        get() = List(colNum) { j -> MutableBuffer(type, rowNum) { i -> get(i, j) } }
+        get() = List(colNum) { j ->
+            MutableStructureNDAccessorBuffer(this, rowNum) { i -> intArrayOf(i, j) }
+        }
 }
 
 /**

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

@@ -5,9 +5,15 @@
 
 package space.kscience.kmath.nd
 
+import space.kscience.attributes.SafeType
 import space.kscience.kmath.PerformancePitfall
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.Int32Field
+import space.kscience.kmath.structures.Float64
 
 public interface StructureNDOfDouble : StructureND<Double> {
+    override val type: SafeType<Float64> get() = Float64Field.type
+
     /**
      * Guaranteed non-blocking access to content
      */
@@ -22,6 +28,7 @@ public fun StructureND<Double>.getDouble(index: IntArray): Double =
     if (this is StructureNDOfDouble) getDouble(index) else get(index)
 
 public interface MutableStructureNDOfDouble : StructureNDOfDouble, MutableStructureND<Double> {
+
     /**
      * Guaranteed non-blocking access to content
      */
@@ -34,6 +41,9 @@ public fun MutableStructureND<Double>.getDouble(index: IntArray): Double =
 
 
 public interface StructureNDOfInt : StructureND<Int> {
+
+    override val type: SafeType<Int> get() = Int32Field.type
+
     /**
      * Guaranteed non-blocking access to content
      */

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt

@@ -77,13 +77,11 @@ private inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.zipInline(
     return elementBufferFactory(l.size) { elementAlgebra.block(l[it], r[it]) }
 }
 
-public fun <T> BufferAlgebra<T, *>.buffer(size: Int, initializer: (Int) -> T): MutableBuffer<T> {
-    return elementBufferFactory(size, initializer)
-}
+public fun <T> BufferAlgebra<T, *>.buffer(size: Int, initializer: (Int) -> T): MutableBuffer<T> =
+    elementBufferFactory(size, initializer)
 
-public fun <T, A> A.buffer(initializer: (Int) -> T): MutableBuffer<T> where A : BufferAlgebra<T, *>, A : WithSize {
-    return elementBufferFactory(size, initializer)
-}
+public fun <T, A> A.buffer(initializer: (Int) -> T): MutableBuffer<T> where A : BufferAlgebra<T, *>, A : WithSize =
+    elementBufferFactory(size, initializer)
 
 public fun <T, A : TrigonometricOperations<T>> BufferAlgebra<T, A>.sin(arg: Buffer<T>): Buffer<T> =
     mapInline(arg) { sin(it) }

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

@@ -42,6 +42,8 @@ public inline fun <reified T> BufferFactory(): BufferFactory<T> = BufferFactory(
  */
 public interface MutableBufferFactory<T> : BufferFactory<T> {
     override fun invoke(size: Int, builder: (Int) -> T): MutableBuffer<T>
+
+    public companion object
 }
 
 /**

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

@@ -6,12 +6,7 @@
 package space.kscience.kmath.structures
 
 import space.kscience.attributes.SafeType
-import space.kscience.kmath.nd.BufferND
-import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.nd.Structure2D
-import space.kscience.kmath.nd.as2D
-import kotlin.collections.component1
-import kotlin.collections.component2
 
 /**
  * A context that allows to operate on a [MutableBuffer] as on 2d array
@@ -32,10 +27,10 @@ internal class BufferAccessor2D<T>(
 
     fun create(mat: Structure2D<T>): MutableBuffer<T> = create { i, j -> mat[i, j] }
 
-    //TODO optimize wrapper
-    fun MutableBuffer<T>.toStructure2D(): Structure2D<T> = BufferND(
-        type, ShapeND(rowNum, colNum)
-    ) { (i, j) -> get(i, j) }.as2D()
+//    //TODO optimize wrapper
+//    fun MutableBuffer<T>.toStructure2D(): Structure2D<T> = BufferND(
+//        type, ShapeND(rowNum, colNum)
+//    ) { (i, j) -> get(i, j) }.as2D()
 
     inner class Row(val buffer: MutableBuffer<T>, val rowIndex: Int) : MutableBuffer<T> {
         override val type: SafeType<T> get() = buffer.type

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

@@ -39,6 +39,7 @@ public value class Float64Buffer(public val array: DoubleArray) : PrimitiveBuffe
     }
 }
 
+@Deprecated("Use Float64Buffer instead", ReplaceWith("Float64Buffer"))
 public typealias DoubleBuffer = Float64Buffer
 
 /**

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

@@ -77,7 +77,6 @@ public inline fun <T> MutableBuffer(
     size: Int,
     initializer: (Int) -> T,
 ): MutableBuffer<T> = when (type.kType) {
-    typeOf<Boolean>() -> TODO()
     typeOf<Int8>() -> Int8Buffer(size) { initializer(it) as Int8 } as MutableBuffer<T>
     typeOf<Int16>() -> MutableBuffer.short(size) { initializer(it) as Int16 } as MutableBuffer<T>
     typeOf<Int32>() -> MutableBuffer.int(size) { initializer(it) as Int32 } as MutableBuffer<T>

kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/DoubleLUSolverTest.kt

@@ -10,7 +10,6 @@ import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.operations.algebra
 import kotlin.test.Test
-import kotlin.test.assertEquals
 import kotlin.test.assertTrue
 
 @OptIn(PerformancePitfall::class)
@@ -29,18 +28,18 @@ class DoubleLUSolverTest {
     }
 
     @Test
-    fun testDecomposition() = Double.algebra.linearSpace.run {
+    fun testDecomposition() = with(Double.algebra.linearSpace){
         val matrix = matrix(2, 2)(
             3.0, 1.0,
             2.0, 3.0
         )
 
-        val lup = lup(matrix)
+        val lup = elementAlgebra.lup(matrix)
 
         //Check determinant
-        assertEquals(7.0, lup.determinant)
+//        assertEquals(7.0, lup.determinant)
 
-        assertMatrixEquals(lup.pivotMatrix dot matrix, lup.l dot lup.u)
+        assertMatrixEquals(lup.pivotMatrix(this) dot matrix, lup.l dot lup.u)
     }
 
     @Test

kmath-core/src/commonTest/kotlin/space/kscience/kmath/linear/MatrixTest.kt

@@ -8,6 +8,7 @@ package space.kscience.kmath.linear
 import space.kscience.kmath.PerformancePitfall
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.algebra
 import kotlin.test.Test
 import kotlin.test.assertEquals
@@ -21,7 +22,7 @@ class MatrixTest {
     @Test
     fun testTranspose() = Double.algebra.linearSpace.run {
         val matrix = one(3, 3)
-        val transposed = matrix.transposed
+        val transposed = matrix.transposed()
         assertTrue { StructureND.contentEquals(matrix, transposed) }
     }
 
@@ -58,7 +59,7 @@ class MatrixTest {
     }
 
     @Test
-    fun test2DDot() = Double.algebra.linearSpace.run {
+    fun test2DDot() = Float64Field.linearSpace {
         val firstMatrix = buildMatrix(2, 3) { i, j -> (i + j).toDouble() }
         val secondMatrix = buildMatrix(3, 2) { i, j -> (i + j).toDouble() }
 
@@ -70,6 +71,5 @@ class MatrixTest {
         assertEquals(8.0, result[0, 1])
         assertEquals(8.0, result[1, 0])
         assertEquals(14.0, result[1, 1])
-
     }
 }

kmath-core/src/jvmMain/kotlin/space/kscience/kmath/linear/Float64ParallelLinearSpace.kt

@@ -0,0 +1,125 @@
+/*
+ * Copyright 2018-2024 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.kmath.linear
+
+import space.kscience.kmath.PerformancePitfall
+import space.kscience.kmath.nd.*
+import space.kscience.kmath.operations.Float64BufferOps
+import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.Float64Buffer
+import space.kscience.kmath.structures.asBuffer
+import java.util.stream.IntStream
+
+
+public object Float64ParallelLinearSpace : LinearSpace<Double, Float64Field> {
+
+    override val elementAlgebra: Float64Field = Float64Field
+
+
+    override fun buildMatrix(
+        rows: Int,
+        columns: Int,
+        initializer: Float64Field.(i: Int, j: Int) -> Double,
+    ): Matrix<Double> {
+        val shape = ShapeND(rows, columns)
+        val indexer = BufferAlgebraND.defaultIndexerBuilder(shape)
+
+        val buffer = IntStream.range(0, indexer.linearSize).parallel().mapToDouble { offset ->
+            val (i, j) = indexer.index(offset)
+            elementAlgebra.initializer(i, j)
+        }.toArray().asBuffer()
+
+        return MutableBufferND(
+            indexer,
+            buffer
+        ).as2D()
+    }
+
+    override fun buildVector(size: Int, initializer: Float64Field.(Int) -> Double): Float64Buffer =
+        IntStream.range(0, size).parallel().mapToDouble{ Float64Field.initializer(it) }.toArray().asBuffer()
+
+    override fun Matrix<Double>.unaryMinus(): Matrix<Double> = Floa64FieldOpsND {
+        asND().map { -it }.as2D()
+    }
+
+    override fun Matrix<Double>.plus(other: Matrix<Double>): Matrix<Double> = Floa64FieldOpsND {
+        require(shape.contentEquals(other.shape)) { "Shape mismatch on Matrix::plus. Expected $shape but found ${other.shape}" }
+        asND().plus(other.asND()).as2D()
+    }
+
+    override fun Matrix<Double>.minus(other: Matrix<Double>): Matrix<Double> = Floa64FieldOpsND {
+        require(shape.contentEquals(other.shape)) { "Shape mismatch on Matrix::minus. Expected $shape but found ${other.shape}" }
+        asND().minus(other.asND()).as2D()
+    }
+
+    // Create a continuous in-memory representation of this vector for better memory layout handling
+    private fun Buffer<Double>.linearize() = if (this is Float64Buffer) {
+        this.array
+    } else {
+        DoubleArray(size) { get(it) }
+    }
+
+    @OptIn(PerformancePitfall::class)
+    override fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> {
+        require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
+        val rows = this@dot.rows.map { it.linearize() }
+        val columns = other.columns.map { it.linearize() }
+        val indices = 0 until this.colNum
+        return buildMatrix(rowNum, other.colNum) { i, j ->
+            val r = rows[i]
+            val c = columns[j]
+            var res = 0.0
+            for (l in indices) {
+                res += r[l] * c[l]
+            }
+            res
+        }
+    }
+
+    @OptIn(PerformancePitfall::class)
+    override fun Matrix<Double>.dot(vector: Point<Double>): Float64Buffer {
+        require(colNum == vector.size) { "Matrix dot vector operation dimension mismatch: ($rowNum, $colNum) x (${vector.size})" }
+        val rows = this@dot.rows.map { it.linearize() }
+        val indices = 0 until this.colNum
+        return Float64Buffer(rowNum) { i ->
+            val r = rows[i]
+            var res = 0.0
+            for (j in indices) {
+                res += r[j] * vector[j]
+            }
+            res
+        }
+
+    }
+
+    override fun Matrix<Double>.times(value: Double): Matrix<Double> = Floa64FieldOpsND {
+        asND().map { it * value }.as2D()
+    }
+
+    public override fun Point<Double>.plus(other: Point<Double>): Float64Buffer = Float64BufferOps.run {
+        this@plus + other
+    }
+
+    public override fun Point<Double>.minus(other: Point<Double>): Float64Buffer = Float64BufferOps.run {
+        this@minus - other
+    }
+
+    public override fun Point<Double>.times(value: Double): Float64Buffer = Float64BufferOps.run {
+        scale(this@times, value)
+    }
+
+    public operator fun Point<Double>.div(value: Double): Float64Buffer = Float64BufferOps.run {
+        scale(this@div, 1.0 / value)
+    }
+
+    public override fun Double.times(v: Point<Double>): Float64Buffer = v * this
+
+
+}
+
+public val Float64LinearSpace.parallel: Float64ParallelLinearSpace get() = Float64ParallelLinearSpace

kmath-core/src/jvmMain/kotlin/space/kscience/kmath/structures/parallelMutableBuffer.kt

@@ -0,0 +1,58 @@
+/*
+ * Copyright 2018-2024 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.kmath.structures
+
+import space.kscience.attributes.SafeType
+import java.util.stream.Collectors
+import java.util.stream.IntStream
+import kotlin.reflect.typeOf
+
+/**
+ * A Java stream-based parallel version of [MutableBuffer].
+ * There is no parallelization for [Int8], [Int16] and [Float32] types.
+ * They are processed sequentially.
+ */
+@Suppress("UNCHECKED_CAST")
+public fun <T> MutableBuffer.Companion.parallel(
+    type: SafeType<T>,
+    size: Int,
+    initializer: (Int) -> T,
+): MutableBuffer<T> = when (type.kType) {
+    typeOf<Int8>() -> Int8Buffer(size) { initializer(it) as Int8 } as MutableBuffer<T>
+    typeOf<Int16>() -> Int16Buffer(size) { initializer(it) as Int16 } as MutableBuffer<T>
+    typeOf<Int32>() -> IntStream.range(0, size).parallel().map { initializer(it) as Int32 }.toArray()
+        .asBuffer() as MutableBuffer<T>
+
+    typeOf<Int64>() -> IntStream.range(0, size).parallel().mapToLong { initializer(it) as Int64 }.toArray()
+        .asBuffer() as MutableBuffer<T>
+
+    typeOf<Float>() -> Float32Buffer(size) { initializer(it) as Float } as MutableBuffer<T>
+    typeOf<Double>() -> IntStream.range(0, size).parallel().mapToDouble { initializer(it) as Float64 }.toArray()
+        .asBuffer() as MutableBuffer<T>
+    //TODO add unsigned types
+    else -> IntStream.range(0, size).parallel().mapToObj { initializer(it) }.collect(Collectors.toList<T>()).asMutableBuffer(type)
+}
+
+public class ParallelBufferFactory<T>(override val type: SafeType<T>) : MutableBufferFactory<T> {
+    override fun invoke(size: Int, builder: (Int) -> T): MutableBuffer<T> {
+        TODO("Not yet implemented")
+    }
+
+}
+
+/**
+ * A Java stream-based parallel alternative to [MutableBufferFactory].
+ * See [MutableBuffer.Companion.parallel] for details.
+ */
+public fun <T> MutableBufferFactory.Companion.parallel(
+    type: SafeType<T>,
+): MutableBufferFactory<T> = object : MutableBufferFactory<T> {
+
+    override val type: SafeType<T> = type
+
+    override fun invoke(size: Int, builder: (Int) -> T): MutableBuffer<T> = MutableBuffer.parallel(type, size, builder)
+
+}

kmath-core/src/jvmTest/kotlin/space/kscience/kmath/linear/ParallelMatrixTest.kt

@@ -0,0 +1,74 @@
+/*
+ * Copyright 2018-2024 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.kmath.linear
+
+import space.kscience.kmath.PerformancePitfall
+import space.kscience.kmath.UnstableKMathAPI
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.operations.Float64Field
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertTrue
+
+@UnstableKMathAPI
+@OptIn(PerformancePitfall::class)
+@Suppress("UNUSED_VARIABLE")
+class ParallelMatrixTest {
+
+    @Test
+    fun testTranspose() =  Float64Field.linearSpace.parallel{
+        val matrix = one(3, 3)
+        val transposed = matrix.transposed()
+        assertTrue { StructureND.contentEquals(matrix, transposed) }
+    }
+
+    @Test
+    fun testBuilder() = Float64Field.linearSpace.parallel{
+        val matrix = matrix(2, 3)(
+            1.0, 0.0, 0.0,
+            0.0, 1.0, 2.0
+        )
+
+        assertEquals(2.0, matrix[1, 2])
+    }
+
+    @Test
+    fun testMatrixExtension() = Float64Field.linearSpace.parallel{
+        val transitionMatrix: Matrix<Double> = VirtualMatrix(type,6, 6) { row, col ->
+            when {
+                col == 0 -> .50
+                row + 1 == col -> .50
+                row == 5 && col == 5 -> 1.0
+                else -> 0.0
+            }
+        }
+
+        infix fun Matrix<Double>.pow(power: Int): Matrix<Double> {
+            var res = this
+            repeat(power - 1) {
+                res =  res dot this@pow
+            }
+            return res
+        }
+
+        val toTenthPower = transitionMatrix pow 10
+    }
+
+    @Test
+    fun test2DDot() = Float64Field.linearSpace.parallel {
+        val firstMatrix = buildMatrix(2, 3) { i, j -> (i + j).toDouble() }
+        val secondMatrix = buildMatrix(3, 2) { i, j -> (i + j).toDouble() }
+
+//            val firstMatrix = produce(2, 3) { i, j -> (i + j).toDouble() }
+//            val secondMatrix = produce(3, 2) { i, j -> (i + j).toDouble() }
+        val result = firstMatrix dot secondMatrix
+        assertEquals(2, result.rowNum)
+        assertEquals(2, result.colNum)
+        assertEquals(8.0, result[0, 1])
+        assertEquals(8.0, result[1, 0])
+        assertEquals(14.0, result[1, 1])
+    }
+}

kmath-coroutines/README.md

@@ -6,19 +6,8 @@
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-coroutines:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-coroutines:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-coroutines:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -27,6 +16,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-coroutines:0.4.0-dev-1")
+    implementation("space.kscience:kmath-coroutines:0.4.0-dev-3")
 }
 ```

kmath-coroutines/build.gradle.kts

@@ -6,6 +6,7 @@ kscience {
     jvm()
     js()
     native()
+    wasm()
 
     dependencies {
         api(project(":kmath-core"))

kmath-coroutines/src/jvmTest/kotlin/space/kscience/kmath/streaming/RingBufferTest.kt

@@ -7,6 +7,7 @@ package space.kscience.kmath.streaming
 
 import kotlinx.coroutines.flow.*
 import kotlinx.coroutines.runBlocking
+import space.kscience.kmath.operations.Int32Ring
 import space.kscience.kmath.operations.asSequence
 import kotlin.test.Test
 import kotlin.test.assertEquals
@@ -14,7 +15,7 @@ import kotlin.test.assertEquals
 internal class RingBufferTest {
     @Test
     fun push() {
-        val buffer = RingBuffer.build(20, Double.NaN)
+        val buffer = RingBuffer(20, Double.NaN)
         runBlocking {
             for (i in 1..30) {
                 buffer.push(i.toDouble())
@@ -30,7 +31,7 @@ internal class RingBufferTest {
             while (true) emit(i++)
         }
 
-        val windowed = flow.windowed(10)
+        val windowed = flow.windowed(10, Int32Ring)
 
         runBlocking {
             @Suppress("UNUSED_VARIABLE") val first = windowed.take(1).single()

kmath-dimensions/README.md

@@ -6,19 +6,8 @@ A proof of concept module for adding type-safe dimensions to structures
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-dimensions:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-dimensions:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-dimensions:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -27,6 +16,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-dimensions:0.4.0-dev-1")
+    implementation("space.kscience:kmath-dimensions:0.4.0-dev-3")
 }
 ```

kmath-dimensions/build.gradle.kts

@@ -6,6 +6,7 @@ kscience{
     jvm()
     js()
     native()
+    wasm()
 
     dependencies{
         api(projects.kmathCore)

kmath-dimensions/src/commonMain/kotlin/space/kscience/kmath/dimensions/Wrappers.kt

@@ -153,7 +153,7 @@ public value class DMatrixContext<T : Any, out A : Ring<T>>(public val context:
         context.run { this@unaryMinus.unaryMinus() }.coerce()
 
     public inline fun <reified R : Dimension, reified C : Dimension> DMatrix<T, C, R>.transposed(): DMatrix<T, R, C> =
-        context.run { (this@transposed as Matrix<T>).transposed }.coerce()
+        context.run { (this@transposed as Matrix<T>).transposed() }.coerce()
 
     public companion object {
         public val real: DMatrixContext<Double, Float64Field> = DMatrixContext(Double.algebra.linearSpace)

kmath-dimensions/src/wasmJsMain/kotlin/space/kscience/kmath/dimensions/Dimension.wasmJs.kt

@@ -0,0 +1,23 @@
+/*
+ * Copyright 2018-2024 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.kmath.dimensions
+
+import kotlin.reflect.KClass
+
+private val dimensionMap: MutableMap<Int, Dimension> = hashMapOf(1 to D1, 2 to D2, 3 to D3)
+
+@Suppress("UNCHECKED_CAST")
+public actual fun <D : Dimension> Dimension.Companion.resolve(type: KClass<D>): D = dimensionMap
+    .entries
+    .map(MutableMap.MutableEntry<Int, Dimension>::value)
+    .find { it::class == type } as? D
+    ?: error("Can't resolve dimension $type")
+
+public actual fun Dimension.Companion.of(dim: Int): Dimension = dimensionMap.getOrPut(dim) {
+    object : Dimension {
+        override val dim: Int get() = dim
+    }
+}

kmath-ejml/README.md

@@ -9,19 +9,8 @@ EJML based linear algebra implementation.
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-ejml:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-ejml:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-ejml:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -30,6 +19,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-ejml:0.4.0-dev-1")
+    implementation("space.kscience:kmath-ejml:0.4.0-dev-3")
 }
 ```

kmath-ejml/build.gradle.kts

@@ -1,15 +1,15 @@
-import space.kscience.kmath.ejml.codegen.ejmlCodegen
-
 plugins {
     id("space.kscience.gradle.jvm")
 }
 
+val ejmlVerision = "0.43.1"
+
 dependencies {
-    api("org.ejml:ejml-ddense:0.41")
-    api("org.ejml:ejml-fdense:0.41")
-    api("org.ejml:ejml-dsparse:0.41")
-    api("org.ejml:ejml-fsparse:0.41")
-    api(project(":kmath-core"))
+    api("org.ejml:ejml-ddense:$ejmlVerision")
+    api("org.ejml:ejml-fdense:$ejmlVerision")
+    api("org.ejml:ejml-dsparse:$ejmlVerision")
+    api("org.ejml:ejml-fsparse:$ejmlVerision")
+    api(projects.kmathCore)
 }
 
 readme {
@@ -32,10 +32,10 @@ readme {
     ) { "LinearSpace implementations." }
 }
 
-kotlin.sourceSets.main {
-    val codegen by tasks.creating {
-        ejmlCodegen(kotlin.srcDirs.first().absolutePath + "/space/kscience/kmath/ejml/_generated.kt")
-    }
-
-    kotlin.srcDirs(files().builtBy(codegen))
-}
+//kotlin.sourceSets.main {
+//    val codegen by tasks.creating {
+//        ejmlCodegen(kotlin.srcDirs.first().absolutePath + "/space/kscience/kmath/ejml/_generated.kt")
+//    }
+//
+//    kotlin.srcDirs(files().builtBy(codegen))
+//}

kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/implementations.kt

@@ -17,21 +17,18 @@ 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.attributes.SafeType
 import space.kscience.attributes.safeTypeOf
-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.nd.Structure2D
+import space.kscience.kmath.nd.StructureAttribute
 import space.kscience.kmath.operations.Float32Field
 import space.kscience.kmath.operations.Float64Field
 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
+import space.kscience.kmath.structures.Float32
+import space.kscience.kmath.structures.IntBuffer
+import space.kscience.kmath.structures.asBuffer
 
 /**
  * [EjmlVector] specialization for [Double].
@@ -78,7 +75,6 @@ public class EjmlFloatMatrix<out M : FMatrix>(override val origin: M) : EjmlMatr
 }
 
 
-
 /**
  * [EjmlLinearSpace] implementation based on [CommonOps_DDRM], [DecompositionFactory_DDRM] operations and
  * [DMatrixRMaj] matrices.
@@ -143,10 +139,10 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
         val out = DMatrixRMaj(1, 1)
 
         CommonOps_DDRM.add(
-            elementAlgebra.one, 
-            toEjml().origin, 
+            elementAlgebra.one,
+            toEjml().origin,
             elementAlgebra { -one },
-            other.toEjml().origin, 
+            other.toEjml().origin,
             out,
         )
 
@@ -167,12 +163,12 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
 
     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, 
+            other.toEjml().origin,
             out,
         )
 
@@ -217,77 +213,65 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
 
     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 }
+    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Double>, attribute: A): V? {
         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()
-                }
+        val raw: Any? = when (attribute) {
+            Inverted -> {
+                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)
+            Determinant -> CommonOps_DDRM.det(origin)
+            SVD -> object : SingularValueDecomposition<Double> {
+                val ejmlSvd by lazy {
+                    DecompositionFactory_DDRM
+                        .svd(origin.numRows, origin.numCols, true, true, false)
                         .apply { decompose(origin.copy()) }
                 }
+                override val u: Matrix<Double> get() = ejmlSvd.getU(null, false).wrapMatrix()
+
+                override val s: Matrix<Double> get() = ejmlSvd.getW(null).wrapMatrix()
+                override val v: Matrix<Double> get() = ejmlSvd.getV(null, false).wrapMatrix()
+                override val singularValues: Point<Double> get() = ejmlSvd.singularValues.asBuffer()
 
-                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) }
+            QR -> object : QRDecomposition<Double> {
+                val ejmlQr by lazy { DecompositionFactory_DDRM.qr().apply { decompose(origin.copy()) } }
+                override val q: Matrix<Double> get() = ejmlQr.getQ(null, false).wrapMatrix()
+                override val r: Matrix<Double> get() = ejmlQr.getR(null, false).wrapMatrix()
             }
 
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
+            Cholesky -> object : CholeskyDecomposition<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)
+                    cholesky.getT(null).wrapMatrix().withAttribute(LowerTriangular)
                 }
             }
 
-            LupDecompositionFeature::class -> object : LupDecompositionFeature<Double> {
+            LUP -> object : LupDecomposition<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 l: Matrix<Double>
+                    get() = lup.getLower(null).wrapMatrix().withAttribute(LowerTriangular)
 
-                override val u: Matrix<Double> by lazy {
-                    lup.getUpper(null).wrapMatrix().withFeature(UFeature)
-                }
 
-                override val p: Matrix<Double> by lazy { lup.getRowPivot(null).wrapMatrix() }
+                override val u: Matrix<Double>
+                    get() = lup.getUpper(null).wrapMatrix().withAttribute(UpperTriangular)
+                override val pivot: IntBuffer get() = lup.getRowPivotV(null).asBuffer()
             }
 
             else -> null
-        }?.let{
-            type.cast(it)
         }
+
+        @Suppress("UNCHECKED_CAST")
+        return raw as V?
     }
 
     /**
@@ -318,7 +302,6 @@ public object EjmlLinearSpaceDDRM : EjmlLinearSpace<Double, Float64Field, DMatri
 }
 
 
-
 /**
  * [EjmlLinearSpace] implementation based on [CommonOps_FDRM], [DecompositionFactory_FDRM] operations and
  * [FMatrixRMaj] matrices.
@@ -383,10 +366,10 @@ public object EjmlLinearSpaceFDRM : EjmlLinearSpace<Float, Float32Field, FMatrix
         val out = FMatrixRMaj(1, 1)
 
         CommonOps_FDRM.add(
-            elementAlgebra.one, 
-            toEjml().origin, 
+            elementAlgebra.one,
+            toEjml().origin,
             elementAlgebra { -one },
-            other.toEjml().origin, 
+            other.toEjml().origin,
             out,
         )
 
@@ -407,12 +390,12 @@ public object EjmlLinearSpaceFDRM : EjmlLinearSpace<Float, Float32Field, FMatrix
 
     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, 
+            other.toEjml().origin,
             out,
         )
 
@@ -457,77 +440,65 @@ public object EjmlLinearSpaceFDRM : EjmlLinearSpace<Float, Float32Field, FMatrix
 
     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 }
+    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Float32>, attribute: A): V? {
         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()
-                }
+        val raw: Any? = when (attribute) {
+            Inverted -> {
+                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)
+            Determinant -> CommonOps_FDRM.det(origin)
+            SVD -> object : SingularValueDecomposition<Float32> {
+                val ejmlSvd by lazy {
+                    DecompositionFactory_FDRM
+                        .svd(origin.numRows, origin.numCols, true, true, false)
                         .apply { decompose(origin.copy()) }
                 }
+                override val u: Matrix<Float32> get() = ejmlSvd.getU(null, false).wrapMatrix()
+
+                override val s: Matrix<Float32> get() = ejmlSvd.getW(null).wrapMatrix()
+                override val v: Matrix<Float32> get() = ejmlSvd.getV(null, false).wrapMatrix()
+                override val singularValues: Point<Float32> get() = ejmlSvd.singularValues.asBuffer()
 
-                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) }
+            QR -> object : QRDecomposition<Float32> {
+                val ejmlQr by lazy { DecompositionFactory_FDRM.qr().apply { decompose(origin.copy()) } }
+                override val q: Matrix<Float32> get() = ejmlQr.getQ(null, false).wrapMatrix()
+                override val r: Matrix<Float32> get() = ejmlQr.getR(null, false).wrapMatrix()
             }
 
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Float> {
-                override val l: Matrix<Float> by lazy {
+            Cholesky -> object : CholeskyDecomposition<Float32> {
+                override val l: Matrix<Float32> by lazy {
                     val cholesky =
                         DecompositionFactory_FDRM.chol(structure.rowNum, true).apply { decompose(origin.copy()) }
 
-                    cholesky.getT(null).wrapMatrix().withFeature(LFeature)
+                    cholesky.getT(null).wrapMatrix().withAttribute(LowerTriangular)
                 }
             }
 
-            LupDecompositionFeature::class -> object : LupDecompositionFeature<Float> {
+            LUP -> object : LupDecomposition<Float32> {
                 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 l: Matrix<Float32>
+                    get() = lup.getLower(null).wrapMatrix().withAttribute(LowerTriangular)
 
-                override val u: Matrix<Float> by lazy {
-                    lup.getUpper(null).wrapMatrix().withFeature(UFeature)
-                }
 
-                override val p: Matrix<Float> by lazy { lup.getRowPivot(null).wrapMatrix() }
+                override val u: Matrix<Float32>
+                    get() = lup.getUpper(null).wrapMatrix().withAttribute(UpperTriangular)
+                override val pivot: IntBuffer get() = lup.getRowPivotV(null).asBuffer()
             }
 
             else -> null
-        }?.let{
-            type.cast(it)
         }
+
+        @Suppress("UNCHECKED_CAST")
+        return raw as V?
     }
 
     /**
@@ -558,7 +529,6 @@ public object EjmlLinearSpaceFDRM : EjmlLinearSpace<Float, Float32Field, FMatrix
 }
 
 
-
 /**
  * [EjmlLinearSpace] implementation based on [CommonOps_DSCC], [DecompositionFactory_DSCC] operations and
  * [DMatrixSparseCSC] matrices.
@@ -623,12 +593,12 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
         val out = DMatrixSparseCSC(1, 1)
 
         CommonOps_DSCC.add(
-            elementAlgebra.one, 
-            toEjml().origin, 
+            elementAlgebra.one,
+            toEjml().origin,
             elementAlgebra { -one },
-            other.toEjml().origin, 
+            other.toEjml().origin,
             out,
-            null, 
+            null,
             null,
         )
 
@@ -649,14 +619,14 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
 
     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, 
+            other.toEjml().origin,
             out,
-            null, 
+            null,
             null,
         )
 
@@ -672,7 +642,7 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
             elementAlgebra.one,
             other.toEjml().origin,
             out,
-            null, 
+            null,
             null,
         )
 
@@ -688,7 +658,7 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
             elementAlgebra { -one },
             other.toEjml().origin,
             out,
-            null, 
+            null,
             null,
         )
 
@@ -705,64 +675,52 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
 
     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 }
+    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Double>, attribute: A): V? {
         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) }
+        val raw: Any? = when (attribute) {
+            Inverted -> {
+                val res = DMatrixRMaj(origin.numRows,origin.numCols)
+                CommonOps_DSCC.invert(origin,res)
+                res.wrapMatrix()
             }
 
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Double> {
+            Determinant -> CommonOps_DSCC.det(origin)
+
+            QR -> object : QRDecomposition<Double> {
+                val ejmlQr by lazy { DecompositionFactory_DSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) } }
+                override val q: Matrix<Double> get() = ejmlQr.getQ(null, false).wrapMatrix()
+                override val r: Matrix<Double> get() = ejmlQr.getR(null, false).wrapMatrix()
+            }
+
+            Cholesky -> object : CholeskyDecomposition<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)
+                    (cholesky.getT(null) as DMatrix).wrapMatrix().withAttribute(LowerTriangular)
                 }
             }
 
-            LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
-                LUDecompositionFeature<Double>, DeterminantFeature<Double>, InverseMatrixFeature<Double> {
-                private val lu by lazy {
+            LUP -> object : LupDecomposition<Double> {
+                private val lup 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 l: Matrix<Double>
+                    get() = lup.getLower(null).wrapMatrix().withAttribute(LowerTriangular)
 
-                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) }
+                override val u: Matrix<Double>
+                    get() = lup.getUpper(null).wrapMatrix().withAttribute(UpperTriangular)
+                override val pivot: IntBuffer get() = lup.getRowPivotV(null).asBuffer()
             }
 
             else -> null
-        }?.let{
-            type.cast(it)
         }
+
+        @Suppress("UNCHECKED_CAST")
+        return raw as V?
     }
 
     /**
@@ -793,7 +751,6 @@ public object EjmlLinearSpaceDSCC : EjmlLinearSpace<Double, Float64Field, DMatri
 }
 
 
-
 /**
  * [EjmlLinearSpace] implementation based on [CommonOps_FSCC], [DecompositionFactory_FSCC] operations and
  * [FMatrixSparseCSC] matrices.
@@ -858,12 +815,12 @@ public object EjmlLinearSpaceFSCC : EjmlLinearSpace<Float, Float32Field, FMatrix
         val out = FMatrixSparseCSC(1, 1)
 
         CommonOps_FSCC.add(
-            elementAlgebra.one, 
-            toEjml().origin, 
+            elementAlgebra.one,
+            toEjml().origin,
             elementAlgebra { -one },
-            other.toEjml().origin, 
+            other.toEjml().origin,
             out,
-            null, 
+            null,
             null,
         )
 
@@ -884,14 +841,14 @@ public object EjmlLinearSpaceFSCC : EjmlLinearSpace<Float, Float32Field, FMatrix
 
     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, 
+            other.toEjml().origin,
             out,
-            null, 
+            null,
             null,
         )
 
@@ -907,7 +864,7 @@ public object EjmlLinearSpaceFSCC : EjmlLinearSpace<Float, Float32Field, FMatrix
             elementAlgebra.one,
             other.toEjml().origin,
             out,
-            null, 
+            null,
             null,
         )
 
@@ -923,7 +880,7 @@ public object EjmlLinearSpaceFSCC : EjmlLinearSpace<Float, Float32Field, FMatrix
             elementAlgebra { -one },
             other.toEjml().origin,
             out,
-            null, 
+            null,
             null,
         )
 
@@ -939,65 +896,52 @@ public object EjmlLinearSpaceFSCC : EjmlLinearSpace<Float, Float32Field, FMatrix
     }
 
     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 }
+    override fun <V, A : StructureAttribute<V>> computeAttribute(structure: Structure2D<Float32>, attribute: A): V? {
         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) }
+        val raw: Any? = when (attribute) {
+            Inverted -> {
+                val res = FMatrixRMaj(origin.numRows,origin.numCols)
+                CommonOps_FSCC.invert(origin,res)
+                res.wrapMatrix()
             }
 
-            CholeskyDecompositionFeature::class -> object : CholeskyDecompositionFeature<Float> {
-                override val l: Matrix<Float> by lazy {
+            Determinant -> CommonOps_FSCC.det(origin)
+
+            QR -> object : QRDecomposition<Float32> {
+                val ejmlQr by lazy { DecompositionFactory_FSCC.qr(FillReducing.NONE).apply { decompose(origin.copy()) } }
+                override val q: Matrix<Float32> get() = ejmlQr.getQ(null, false).wrapMatrix()
+                override val r: Matrix<Float32> get() = ejmlQr.getR(null, false).wrapMatrix()
+            }
+
+            Cholesky -> object : CholeskyDecomposition<Float32> {
+                override val l: Matrix<Float32> by lazy {
                     val cholesky =
                         DecompositionFactory_FSCC.cholesky().apply { decompose(origin.copy()) }
 
-                    (cholesky.getT(null) as FMatrix).wrapMatrix().withFeature(LFeature)
+                    (cholesky.getT(null) as FMatrix).wrapMatrix().withAttribute(LowerTriangular)
                 }
             }
 
-            LUDecompositionFeature::class, DeterminantFeature::class, InverseMatrixFeature::class -> object :
-                LUDecompositionFeature<Float>, DeterminantFeature<Float>, InverseMatrixFeature<Float> {
-                private val lu by lazy {
+            LUP -> object : LupDecomposition<Float32> {
+                private val lup 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 l: Matrix<Float32>
+                    get() = lup.getLower(null).wrapMatrix().withAttribute(LowerTriangular)
 
-                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) }
+                override val u: Matrix<Float32>
+                    get() = lup.getUpper(null).wrapMatrix().withAttribute(UpperTriangular)
+                override val pivot: IntBuffer get() = lup.getRowPivotV(null).asBuffer()
             }
 
             else -> null
-        }?.let{
-            type.cast(it)
         }
+
+        @Suppress("UNCHECKED_CAST")
+        return raw as V?
     }
 
     /**

kmath-ejml/src/test/kotlin/space/kscience/kmath/ejml/EjmlMatrixTest.kt

@@ -58,19 +58,19 @@ internal class EjmlMatrixTest {
 
     @OptIn(UnstableKMathAPI::class)
     @Test
-    fun features() {
+    fun features() = EjmlLinearSpaceDDRM {
         val m = randomMatrix
         val w = EjmlDoubleMatrix(m)
-        val det: Determinant<Double> = EjmlLinearSpaceDDRM.attributeForOrNull(w) ?: fail()
-        assertEquals(CommonOps_DDRM.det(m), det.determinant)
-        val lup: LupDecompositionAttribute<Double> = EjmlLinearSpaceDDRM.attributeForOrNull(w) ?: fail()
+        val det: Double = w.getOrComputeAttribute(Determinant) ?: fail()
+        assertEquals(CommonOps_DDRM.det(m), det)
+        val lup: LupDecomposition<Double> = w.getOrComputeAttribute(LUP) ?: fail()
 
         val ludecompositionF64 = DecompositionFactory_DDRM.lu(m.numRows, m.numCols)
             .also { it.decompose(m.copy()) }
 
         assertMatrixEquals(EjmlDoubleMatrix(ludecompositionF64.getLower(null)), lup.l)
         assertMatrixEquals(EjmlDoubleMatrix(ludecompositionF64.getUpper(null)), lup.u)
-        assertMatrixEquals(EjmlDoubleMatrix(ludecompositionF64.getRowPivot(null)), lup.p)
+        assertMatrixEquals(EjmlDoubleMatrix(ludecompositionF64.getRowPivot(null)), lup.pivotMatrix(this))
     }
 
     @Test

kmath-for-real/README.md

@@ -9,19 +9,8 @@ Specialization of KMath APIs for Double numbers.
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-for-real:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-for-real:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-for-real:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -30,6 +19,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-for-real:0.4.0-dev-1")
+    implementation("space.kscience:kmath-for-real:0.4.0-dev-3")
 }
 ```

kmath-for-real/build.gradle.kts

@@ -6,6 +6,7 @@ kscience {
     jvm()
     js()
     native()
+    wasm()
 
     dependencies {
         api(projects.kmathCore)

kmath-for-real/src/commonTest/kotlin/space/kscience/kmath/real/DoubleVectorTest.kt

@@ -7,7 +7,7 @@ package space.kscience.kmath.real
 
 import space.kscience.kmath.linear.asMatrix
 import space.kscience.kmath.linear.linearSpace
-import space.kscience.kmath.linear.transpose
+import space.kscience.kmath.linear.transposed
 import space.kscience.kmath.operations.algebra
 import space.kscience.kmath.structures.Float64Buffer
 import kotlin.test.Test
@@ -34,7 +34,7 @@ internal class DoubleVectorTest {
         val vector1 = Float64Buffer(5) { it.toDouble() }
         val vector2 = Float64Buffer(5) { 5 - it.toDouble() }
         val matrix1 = vector1.asMatrix()
-        val matrix2 = vector2.asMatrix().transpose()
+        val matrix2 = vector2.asMatrix().transposed()
         val product = matrix1 dot matrix2
         assertEquals(5.0, product[1, 0])
         assertEquals(6.0, product[2, 2])

kmath-functions/README.md

@@ -11,19 +11,8 @@ Functions and interpolations.
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-functions:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-functions:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-functions:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -32,6 +21,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-functions:0.4.0-dev-1")
+    implementation("space.kscience:kmath-functions:0.4.0-dev-3")
 }
 ```

kmath-functions/build.gradle.kts

@@ -6,7 +6,6 @@ kscience{
     jvm()
     js()
     native()
-
     wasm()
 
     dependencies {

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/GaussIntegrator.kt

@@ -92,7 +92,7 @@ public inline fun <reified T : Any> GaussIntegrator<T>.integrate(
     range: ClosedRange<Double>,
     order: Int = 10,
     intervals: Int = 10,
-    attributesBuilder: AttributesBuilder<UnivariateIntegrand<T>>.() -> Unit,
+    attributesBuilder: AttributesBuilder<UnivariateIntegrand<T>>.() -> Unit = {},
     noinline function: (Double) -> T,
 ): UnivariateIntegrand<T> {
     require(range.endInclusive > range.start) { "The range upper bound should be higher than lower bound" }

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt

@@ -23,7 +23,7 @@ import space.kscience.kmath.structures.MutableBufferFactory
  */
 public class SplineInterpolator<T : Comparable<T>>(
     override val algebra: Field<T>,
-    public val bufferFactory: MutableBufferFactory<T>,
+    public val bufferFactory: MutableBufferFactory<T> = algebra.bufferFactory,
 ) : PolynomialInterpolator<T> {
     //TODO possibly optimize zeroed buffers
 

kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/testUtils/IntModulo.kt

@@ -9,6 +9,7 @@ package space.kscience.kmath.functions.testUtils
 
 import space.kscience.kmath.operations.Ring
 import space.kscience.kmath.operations.ScaleOperations
+import space.kscience.kmath.structures.MutableBufferFactory
 
 
 class IntModulo {
@@ -109,15 +110,17 @@ class IntModulo {
 }
 
 @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE")
-class IntModuloRing : Ring<IntModulo>, ScaleOperations<IntModulo> {
+class IntModuloRing(modulus: Int) : Ring<IntModulo>, ScaleOperations<IntModulo> {
 
     val modulus: Int
 
-    constructor(modulus: Int) {
+    init {
         require(modulus != 0) { "modulus can not be zero" }
         this.modulus = if (modulus < 0) -modulus else modulus
     }
 
+    override val bufferFactory: MutableBufferFactory<IntModulo> = MutableBufferFactory()
+
     override inline val zero: IntModulo get() = IntModulo(0, modulus, toCheckInput = false)
     override inline val one: IntModulo get() = IntModulo(1, modulus, toCheckInput = false)
 

kmath-functions/src/commonTest/kotlin/space/kscience/kmath/functions/testUtils/Rational.kt

@@ -10,6 +10,7 @@ package space.kscience.kmath.functions.testUtils
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.operations.Field
 import space.kscience.kmath.operations.NumbersAddOps
+import space.kscience.kmath.structures.MutableBufferFactory
 
 @Suppress("NAME_SHADOWING")
 class Rational {
@@ -159,6 +160,7 @@ class Rational {
 @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE")
 @OptIn(UnstableKMathAPI::class)
 object RationalField : Field<Rational>, NumbersAddOps<Rational> {
+    override val bufferFactory: MutableBufferFactory<Rational> = MutableBufferFactory()
     override inline val zero: Rational get() = Rational.ZERO
     override inline val one: Rational get() = Rational.ONE
 

kmath-functions/src/commonTest/kotlin/space/kscience/kmath/integration/SplineIntegralTest.kt

@@ -18,15 +18,15 @@ import kotlin.test.assertEquals
 class SplineIntegralTest {
 
     @Test
-    fun integratePolynomial(){
+    fun integratePolynomial() {
         val polynomial = Polynomial(1.0, 2.0, 3.0)
-        val integral = polynomial.integrate(Float64Field,1.0..2.0)
+        val integral = polynomial.integrate(Float64Field, 1.0..2.0)
         assertEquals(11.0, integral, 0.001)
     }
 
     @Test
     fun gaussSin() {
-        val res = Float64Field.splineIntegrator.integrate(0.0..2 * PI, IntegrandMaxCalls(5)) { x ->
+        val res = Float64Field.splineIntegrator.integrate(0.0..2 * PI, { IntegrandMaxCalls(5) }) { x ->
             sin(x)
         }
         assertEquals(0.0, res.value, 1e-2)
@@ -34,8 +34,8 @@ class SplineIntegralTest {
 
     @Test
     fun gaussUniform() {
-        val res = Float64Field.splineIntegrator.integrate(35.0..100.0, IntegrandMaxCalls(20)) { x ->
-            if(x in 30.0..50.0){
+        val res = Float64Field.splineIntegrator.integrate(35.0..100.0, { IntegrandMaxCalls(20) }) { x ->
+            if (x in 30.0..50.0) {
                 1.0
             } else {
                 0.0

kmath-geometry/README.md

@@ -6,19 +6,8 @@
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-geometry:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-geometry:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-geometry:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -27,6 +16,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-geometry:0.4.0-dev-1")
+    implementation("space.kscience:kmath-geometry:0.4.0-dev-3")
 }
 ```

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/Line.kt

@@ -7,8 +7,8 @@ package space.kscience.kmath.geometry
 
 import kotlinx.serialization.SerialName
 import kotlinx.serialization.Serializable
-import space.kscience.kmath.geometry.euclidean2d.DoubleVector2D
-import space.kscience.kmath.geometry.euclidean3d.DoubleVector3D
+import space.kscience.kmath.geometry.euclidean2d.Float64Vector2D
+import space.kscience.kmath.geometry.euclidean3d.Float64Vector3D
 
 /**
  * A line formed by [start] vector of start and a [direction] vector. Direction vector is not necessarily normalized,
@@ -25,8 +25,8 @@ private data class LineImpl<out V : Any>(override val start: V, override val dir
 
 public fun <V : Any> Line(base: V, direction: V): Line<V> = LineImpl(base, direction)
 
-public typealias Line2D = Line<DoubleVector2D>
-public typealias Line3D = Line<DoubleVector3D>
+public typealias Line2D = Line<Float64Vector2D>
+public typealias Line3D = Line<Float64Vector3D>
 
 /**
  * A directed line segment between [begin] and [end]
@@ -49,5 +49,5 @@ public fun <V : Any> LineSegment<V>.line(algebra: GeometrySpace<V, *>): Line<V>
     Line(begin, end - begin)
 }
 
-public typealias LineSegment2D = LineSegment<DoubleVector2D>
-public typealias LineSegment3D = LineSegment<DoubleVector3D>
+public typealias LineSegment2D = LineSegment<Float64Vector2D>
+public typealias LineSegment3D = LineSegment<Float64Vector3D>

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/Vector3D.kt

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.geometry
 
+import space.kscience.attributes.SafeType
 import space.kscience.kmath.linear.Point
 import space.kscience.kmath.structures.Buffer
 
@@ -33,6 +34,9 @@ public fun <T> Buffer<T>.asVector3D(): Vector3D<T> = object : Vector3D<T> {
         require(this@asVector3D.size == 3) { "Buffer of size 3 is required for Vector3D" }
     }
 
+    override val type: SafeType<T> = this@asVector3D.type
+
+
     override val x: T get() = this@asVector3D[0]
     override val y: T get() = this@asVector3D[1]
     override val z: T get() = this@asVector3D[2]

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/angles.kt

@@ -11,6 +11,8 @@ import kotlinx.serialization.builtins.serializer
 import kotlinx.serialization.descriptors.SerialDescriptor
 import kotlinx.serialization.encoding.Decoder
 import kotlinx.serialization.encoding.Encoder
+import space.kscience.kmath.operations.Group
+import space.kscience.kmath.structures.MutableBufferFactory
 import kotlin.jvm.JvmInline
 import kotlin.math.PI
 import kotlin.math.floor
@@ -28,11 +30,19 @@ public sealed interface Angle : Comparable<Angle> {
     public operator fun div(other: Angle): Double
     public operator fun unaryMinus(): Angle
 
-    public companion object {
-        public val zero: Radians = Radians(0.0)
+    public companion object: Group<Angle> {
+        override val zero: Radians = Radians(0.0)
         public val pi: Radians = Radians(PI)
         public val piTimes2: Radians = Radians(PI * 2)
         public val piDiv2: Radians = Radians(PI / 2)
+
+
+        override fun add(left: Angle, right: Angle): Angle = left + right
+
+        @Suppress("EXTENSION_SHADOWED_BY_MEMBER")
+        override fun Angle.unaryMinus(): Angle = -this
+
+        override val bufferFactory: MutableBufferFactory<Angle> = MutableBufferFactory()
     }
 }
 
@@ -43,7 +53,7 @@ public object AngleSerializer : KSerializer<Angle> {
     override fun deserialize(decoder: Decoder): Angle = decoder.decodeDouble().degrees
 
     override fun serialize(encoder: Encoder, value: Angle) {
-        encoder.encodeDouble(value.degrees)
+        encoder.encodeDouble(value.toDegrees().value)
     }
 }
 
@@ -56,16 +66,16 @@ public value class Radians(public val value: Double) : Angle {
     override fun toRadians(): Radians = this
     override fun toDegrees(): Degrees = Degrees(value * 180 / PI)
 
-    public override fun plus(other: Angle): Radians = Radians(value + other.radians)
-    public override fun minus(other: Angle): Radians = Radians(value - other.radians)
+    public override fun plus(other: Angle): Radians = Radians(value + other.toRadians().value)
+    public override fun minus(other: Angle): Radians = Radians(value - other.toRadians().value)
 
     public override fun times(other: Number): Radians = Radians(value * other.toDouble())
     public override fun div(other: Number): Radians = Radians(value / other.toDouble())
-    override fun div(other: Angle): Double = value / other.radians
+    override fun div(other: Angle): Double = value / other.toRadians().value
 
     public override fun unaryMinus(): Radians = Radians(-value)
 
-    override fun compareTo(other: Angle): Int = value.compareTo(other.radians)
+    override fun compareTo(other: Angle): Int = value.compareTo(other.toRadians().value)
 }
 
 public fun sin(angle: Angle): Double = kotlin.math.sin(angle.toRadians().value)
@@ -85,16 +95,16 @@ public value class Degrees(public val value: Double) : Angle {
     override fun toRadians(): Radians = Radians(value * PI / 180)
     override fun toDegrees(): Degrees = this
 
-    public override fun plus(other: Angle): Degrees = Degrees(value + other.degrees)
-    public override fun minus(other: Angle): Degrees = Degrees(value - other.degrees)
+    public override fun plus(other: Angle): Degrees = Degrees(value + other.toDegrees().value)
+    public override fun minus(other: Angle): Degrees = Degrees(value - other.toDegrees().value)
 
     public override fun times(other: Number): Degrees = Degrees(value * other.toDouble())
     public override fun div(other: Number): Degrees = Degrees(value / other.toDouble())
-    override fun div(other: Angle): Double = value / other.degrees
+    override fun div(other: Angle): Double = value / other.toDegrees().value
 
     public override fun unaryMinus(): Degrees = Degrees(-value)
 
-    override fun compareTo(other: Angle): Int = value.compareTo(other.degrees)
+    override fun compareTo(other: Angle): Int = value.compareTo(other.toDegrees().value)
 }
 
 public val Number.degrees: Degrees get() = Degrees(toDouble())
@@ -106,6 +116,6 @@ public val Angle.degrees: Double get() = toDegrees().value
  * Normalized angle 2 PI range symmetric around [center]. By default, uses (0, 2PI) range.
  */
 public fun Angle.normalized(center: Angle = Angle.pi): Angle =
-    this - Angle.piTimes2 * floor((radians + PI - center.radians) / PI / 2)
+    this - Angle.piTimes2 * floor((toRadians().value + PI - center.toRadians().value) / PI / 2)
 
 public fun abs(angle: Angle): Angle = if (angle < Angle.zero) -angle else angle

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/euclidean2d/Float32Space2D.kt

@@ -11,19 +11,22 @@ import kotlinx.serialization.Serializable
 import kotlinx.serialization.descriptors.SerialDescriptor
 import kotlinx.serialization.encoding.Decoder
 import kotlinx.serialization.encoding.Encoder
+import space.kscience.attributes.SafeType
 import space.kscience.kmath.geometry.GeometrySpace
 import space.kscience.kmath.geometry.Vector2D
 import space.kscience.kmath.operations.Float32Field
 import space.kscience.kmath.structures.Float32
+import space.kscience.kmath.structures.MutableBufferFactory
 import kotlin.math.pow
 import kotlin.math.sqrt
 
 @Serializable(Float32Space2D.VectorSerializer::class)
-public interface Float32Vector2D : Vector2D<Float>
+public interface Float32Vector2D : Vector2D<Float32>{
+    override val type: SafeType<Float32> get() = Float32Field.type
+}
 
 
 public object Float32Space2D : GeometrySpace<Float32Vector2D, Float32> {
-
     @Serializable
     @SerialName("Float32Vector2D")
     private data class Vector2DImpl(
@@ -72,6 +75,8 @@ public object Float32Space2D : GeometrySpace<Float32Vector2D, Float32> {
     public val yAxis: Float32Vector2D = vector(0.0, 1.0)
 
     override val defaultPrecision: Float32 = 1e-3f
+
+    override val bufferFactory: MutableBufferFactory<Float32Vector2D> = MutableBufferFactory()
 }
 
 public fun Float32Vector2D(x: Number, y: Number): Float32Vector2D = Float32Space2D.vector(x, y)

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/euclidean2d/Float64Space2D.kt

@@ -11,66 +11,77 @@ import kotlinx.serialization.Serializable
 import kotlinx.serialization.descriptors.SerialDescriptor
 import kotlinx.serialization.encoding.Decoder
 import kotlinx.serialization.encoding.Encoder
+import space.kscience.attributes.SafeType
 import space.kscience.kmath.geometry.GeometrySpace
 import space.kscience.kmath.geometry.Vector2D
-import space.kscience.kmath.linear.Float64LinearSpace
 import space.kscience.kmath.operations.Float64Field
 import space.kscience.kmath.operations.ScaleOperations
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.structures.MutableBufferFactory
 import kotlin.math.pow
 import kotlin.math.sqrt
 
 
-public typealias DoubleVector2D = Vector2D<Double>
-public typealias Float64Vector2D = Vector2D<Double>
+@Serializable(Float64Space2D.VectorSerializer::class)
+public interface Float64Vector2D : Vector2D<Float64> {
+    override val type: SafeType<Float64> get() = Float64Field.type
+}
+
+@Deprecated("Use Float64Vector2D", ReplaceWith("Float64Vector2D"))
+public typealias DoubleVector2D = Float64Vector2D
 
-public val Vector2D<Double>.r: Double get() = Float64Space2D.norm(this)
 
 
 /**
  * 2D Euclidean space
  */
-public object Float64Space2D : GeometrySpace<DoubleVector2D>, ScaleOperations<DoubleVector2D> {
+public object Float64Space2D : GeometrySpace<Float64Vector2D, Float64>, ScaleOperations<Float64Vector2D> {
+
 
     @Serializable
     @SerialName("Float64Vector2D")
     private data class Vector2DImpl(
         override val x: Double,
         override val y: Double,
-    ) : DoubleVector2D
+    ) : Float64Vector2D
 
-    public object VectorSerializer : KSerializer<DoubleVector2D> {
+    public object VectorSerializer : KSerializer<Float64Vector2D> {
         private val proxySerializer = Vector2DImpl.serializer()
         override val descriptor: SerialDescriptor get() = proxySerializer.descriptor
 
-        override fun deserialize(decoder: Decoder): DoubleVector2D = decoder.decodeSerializableValue(proxySerializer)
+        override fun deserialize(decoder: Decoder): Float64Vector2D = decoder.decodeSerializableValue(proxySerializer)
 
-        override fun serialize(encoder: Encoder, value: DoubleVector2D) {
+        override fun serialize(encoder: Encoder, value: Float64Vector2D) {
             val vector = value as? Vector2DImpl ?: Vector2DImpl(value.x, value.y)
             encoder.encodeSerializableValue(proxySerializer, vector)
         }
     }
 
-    public fun vector(x: Number, y: Number): DoubleVector2D = Vector2DImpl(x.toDouble(), y.toDouble())
+    public fun vector(x: Number, y: Number): Float64Vector2D = Vector2DImpl(x.toDouble(), y.toDouble())
 
-    override val zero: DoubleVector2D by lazy { vector(0.0, 0.0) }
+    override val zero: Float64Vector2D by lazy { vector(0.0, 0.0) }
 
-    override fun norm(arg: DoubleVector2D): Double = sqrt(arg.x.pow(2) + arg.y.pow(2))
+    override fun norm(arg: Float64Vector2D): Double = sqrt(arg.x.pow(2) + arg.y.pow(2))
 
-    override fun DoubleVector2D.unaryMinus(): DoubleVector2D = vector(-x, -y)
+    override fun Float64Vector2D.unaryMinus(): Float64Vector2D = vector(-x, -y)
 
-    override fun DoubleVector2D.distanceTo(other: DoubleVector2D): Double = norm(this - other)
-    override fun add(left: DoubleVector2D, right: DoubleVector2D): DoubleVector2D =
+    override fun Float64Vector2D.distanceTo(other: Float64Vector2D): Double = norm(this - other)
+    override fun add(left: Float64Vector2D, right: Float64Vector2D): Float64Vector2D =
         vector(left.x + right.x, left.y + right.y)
 
-    override fun scale(a: DoubleVector2D, value: Double): DoubleVector2D = vector(a.x * value, a.y * value)
-    override fun DoubleVector2D.dot(other: DoubleVector2D): Double = x * other.x + y * other.y
+    override fun scale(a: Float64Vector2D, value: Double): Float64Vector2D = vector(a.x * value, a.y * value)
+    override fun Float64Vector2D.dot(other: Float64Vector2D): Double = x * other.x + y * other.y
 
-    public val xAxis: DoubleVector2D = vector(1.0, 0.0)
-    public val yAxis: DoubleVector2D = vector(0.0, 1.0)
+    public val xAxis: Float64Vector2D = vector(1.0, 0.0)
+    public val yAxis: Float64Vector2D = vector(0.0, 1.0)
 
     override val defaultPrecision: Double = 1e-6
+
+    override val bufferFactory: MutableBufferFactory<Float64Vector2D> = MutableBufferFactory()
 }
 
 public fun Float64Vector2D(x: Number, y: Number): Float64Vector2D = Float64Space2D.vector(x, y)
 
+public val Float64Vector2D.r: Float64 get() = Float64Space2D.norm(this)
+
 public val Float64Field.euclidean2D: Float64Space2D get() = Float64Space2D

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/euclidean3d/Float32Space3D.kt

@@ -11,15 +11,19 @@ import kotlinx.serialization.Serializable
 import kotlinx.serialization.descriptors.SerialDescriptor
 import kotlinx.serialization.encoding.Decoder
 import kotlinx.serialization.encoding.Encoder
+import space.kscience.attributes.SafeType
 import space.kscience.kmath.geometry.GeometrySpace
 import space.kscience.kmath.geometry.Vector3D
 import space.kscience.kmath.operations.Float32Field
 import space.kscience.kmath.structures.Float32
+import space.kscience.kmath.structures.MutableBufferFactory
 import kotlin.math.pow
 import kotlin.math.sqrt
 
 @Serializable(Float32Space3D.VectorSerializer::class)
-public interface Float32Vector3D : Vector3D<Float>
+public interface Float32Vector3D : Vector3D<Float>{
+    override val type: SafeType<Float32> get() = Float32Field.type
+}
 
 
 public object Float32Space3D : GeometrySpace<Float32Vector3D, Float32> {
@@ -101,6 +105,8 @@ public object Float32Space3D : GeometrySpace<Float32Vector3D, Float32> {
     public val zAxis: Float32Vector3D = vector(0.0, 0.0, 1.0)
 
     override val defaultPrecision: Float32 = 1e-3f
+
+    override val bufferFactory: MutableBufferFactory<Float32Vector3D> = MutableBufferFactory()
 }
 
 public fun Float32Vector3D(x: Number, y: Number, z: Number): Float32Vector3D = Float32Space3D.vector(x, y, z)

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/euclidean3d/Float64Space3D.kt

@@ -11,10 +11,13 @@ import kotlinx.serialization.Serializable
 import kotlinx.serialization.descriptors.SerialDescriptor
 import kotlinx.serialization.encoding.Decoder
 import kotlinx.serialization.encoding.Encoder
+import space.kscience.attributes.SafeType
 import space.kscience.kmath.geometry.GeometrySpace
 import space.kscience.kmath.geometry.Vector3D
 import space.kscience.kmath.linear.Float64LinearSpace
 import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
+import space.kscience.kmath.structures.MutableBufferFactory
 import kotlin.math.pow
 import kotlin.math.sqrt
 
@@ -31,12 +34,16 @@ internal fun leviCivita(i: Int, j: Int, k: Int): Int = when {
     else -> 0
 }
 
-public typealias DoubleVector3D = Vector3D<Double>
-public typealias Float64Vector3D = Vector3D<Double>
+@Serializable(Float64Space3D.VectorSerializer::class)
+public interface Float64Vector3D : Vector3D<Float64> {
+    override val type: SafeType<Float64> get() = Float64Field.type
+}
 
-public val DoubleVector3D.r: Double get() = Float64Space3D.norm(this)
+@Deprecated("Use Float64Vector3D", ReplaceWith("Float64Vector3D"))
+public typealias DoubleVector3D = Float64Vector3D
 
-public object Float64Space3D : GeometrySpace<DoubleVector3D, Double>{
+
+public object Float64Space3D : GeometrySpace<Vector3D<Float64>, Double> {
 
     public val linearSpace: Float64LinearSpace = Float64LinearSpace
 
@@ -46,52 +53,52 @@ public object Float64Space3D : GeometrySpace<DoubleVector3D, Double>{
         override val x: Double,
         override val y: Double,
         override val z: Double,
-    ) : DoubleVector3D
+    ) : Float64Vector3D
 
-    public object VectorSerializer : KSerializer<DoubleVector3D> {
+    public object VectorSerializer : KSerializer<Float64Vector3D> {
         private val proxySerializer = Vector3DImpl.serializer()
         override val descriptor: SerialDescriptor get() = proxySerializer.descriptor
 
-        override fun deserialize(decoder: Decoder): DoubleVector3D = decoder.decodeSerializableValue(proxySerializer)
+        override fun deserialize(decoder: Decoder): Float64Vector3D = decoder.decodeSerializableValue(proxySerializer)
 
-        override fun serialize(encoder: Encoder, value: DoubleVector3D) {
+        override fun serialize(encoder: Encoder, value: Float64Vector3D) {
             val vector = value as? Vector3DImpl ?: Vector3DImpl(value.x, value.y, value.z)
             encoder.encodeSerializableValue(proxySerializer, vector)
         }
     }
 
-    public fun vector(x: Double, y: Double, z: Double): DoubleVector3D =
+    public fun vector(x: Double, y: Double, z: Double): Float64Vector3D =
         Vector3DImpl(x, y, z)
 
-    public fun vector(x: Number, y: Number, z: Number): DoubleVector3D =
+    public fun vector(x: Number, y: Number, z: Number): Float64Vector3D =
         vector(x.toDouble(), y.toDouble(), z.toDouble())
 
-    override val zero: DoubleVector3D by lazy { vector(0.0, 0.0, 0.0) }
+    override val zero: Float64Vector3D by lazy { vector(0.0, 0.0, 0.0) }
 
-    override fun norm(arg: DoubleVector3D): Double = sqrt(arg.x.pow(2) + arg.y.pow(2) + arg.z.pow(2))
+    override fun norm(arg: Vector3D<Float64>): Double = sqrt(arg.x.pow(2) + arg.y.pow(2) + arg.z.pow(2))
 
-    public fun DoubleVector3D.norm(): Double = norm(this)
+    public fun Vector3D<Float64>.norm(): Double = norm(this)
 
-    override fun DoubleVector3D.unaryMinus(): DoubleVector3D = vector(-x, -y, -z)
+    override fun Vector3D<Float64>.unaryMinus(): Float64Vector3D = vector(-x, -y, -z)
 
-    override fun DoubleVector3D.distanceTo(other: DoubleVector3D): Double = (this - other).norm()
+    override fun Vector3D<Float64>.distanceTo(other: Vector3D<Float64>): Double = (this - other).norm()
 
-    override fun add(left: DoubleVector3D, right: DoubleVector3D): DoubleVector3D =
+    override fun add(left: Vector3D<Float64>, right: Vector3D<Float64>): Float64Vector3D =
         vector(left.x + right.x, left.y + right.y, left.z + right.z)
 
-    override fun scale(a: DoubleVector3D, value: Double): DoubleVector3D =
+    override fun scale(a: Vector3D<Float64>, value: Double): Float64Vector3D =
         vector(a.x * value, a.y * value, a.z * value)
 
-    override fun DoubleVector3D.dot(other: DoubleVector3D): Double =
+    override fun Vector3D<Float64>.dot(other: Vector3D<Float64>): Double =
         x * other.x + y * other.y + z * other.z
 
     /**
      * Compute vector product of [first] and [second]. The basis is assumed to be right-handed.
      */
     public fun vectorProduct(
-        first: DoubleVector3D,
-        second: DoubleVector3D,
-    ): DoubleVector3D {
+        first: Vector3D<Float64>,
+        second: Vector3D<Float64>,
+    ): Float64Vector3D {
         var x = 0.0
         var y = 0.0
         var z = 0.0
@@ -110,13 +117,18 @@ public object Float64Space3D : GeometrySpace<DoubleVector3D, Double>{
     /**
      * Vector product with the right basis
      */
-    public infix fun DoubleVector3D.cross(other: DoubleVector3D): Vector3D<Double> = vectorProduct(this, other)
+    public infix fun Vector3D<Float64>.cross(other: Vector3D<Float64>): Vector3D<Double> = vectorProduct(this, other)
 
-    public val xAxis: DoubleVector3D = vector(1.0, 0.0, 0.0)
-    public val yAxis: DoubleVector3D = vector(0.0, 1.0, 0.0)
-    public val zAxis: DoubleVector3D = vector(0.0, 0.0, 1.0)
+    public val xAxis: Float64Vector3D = vector(1.0, 0.0, 0.0)
+    public val yAxis: Float64Vector3D = vector(0.0, 1.0, 0.0)
+    public val zAxis: Float64Vector3D = vector(0.0, 0.0, 1.0)
 
     override val defaultPrecision: Double = 1e-6
+
+    override val bufferFactory: MutableBufferFactory<Vector3D<Float64>> = MutableBufferFactory()
 }
 
 public val Float64Field.euclidean3D: Float64Space3D get() = Float64Space3D
+
+
+public val Float64Vector3D.r: Double get() = Float64Space3D.norm(this)

kmath-geometry/src/commonMain/kotlin/space/kscience/kmath/geometry/euclidean3d/rotations3D.kt

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.geometry.euclidean3d
 
+import space.kscience.attributes.SafeType
 import space.kscience.kmath.UnstableKMathAPI
 import space.kscience.kmath.complex.*
 import space.kscience.kmath.geometry.*
@@ -13,6 +14,7 @@ import space.kscience.kmath.linear.Matrix
 import space.kscience.kmath.linear.linearSpace
 import space.kscience.kmath.linear.matrix
 import space.kscience.kmath.operations.Float64Field
+import space.kscience.kmath.structures.Float64
 import kotlin.math.*
 
 public operator fun Quaternion.times(other: Quaternion): Quaternion = QuaternionAlgebra.multiply(this, other)
@@ -35,7 +37,7 @@ public infix fun Quaternion.dot(other: Quaternion): Double = w * other.w + x * o
 /**
  * Represent a vector as quaternion with zero a rotation angle.
  */
-internal fun DoubleVector3D.asQuaternion(): Quaternion = Quaternion(0.0, x, y, z)
+internal fun Float64Vector3D.asQuaternion(): Quaternion = Quaternion(0.0, x, y, z)
 
 /**
  * Angle in radians denoted by this quaternion rotation
@@ -45,7 +47,7 @@ public val Quaternion.theta: Radians get() = (kotlin.math.acos(normalized().w) *
 /**
  * Create a normalized Quaternion from rotation angle and rotation vector
  */
-public fun Quaternion.Companion.fromRotation(theta: Angle, vector: DoubleVector3D): Quaternion {
+public fun Quaternion.Companion.fromRotation(theta: Angle, vector: Float64Vector3D): Quaternion {
     val s = sin(theta / 2)
     val c = cos(theta / 2)
     val norm = with(Float64Space3D) { vector.norm() }
@@ -55,9 +57,9 @@ public fun Quaternion.Companion.fromRotation(theta: Angle, vector: DoubleVector3
 /**
  * An axis of quaternion rotation
  */
-public val Quaternion.vector: DoubleVector3D
+public val Quaternion.vector: Float64Vector3D
     get() {
-        return object : DoubleVector3D {
+        return object : Float64Vector3D {
             private val sint2 = sqrt(1 - w * w)
             override val x: Double get() = this@vector.x / sint2
             override val y: Double get() = this@vector.y / sint2
@@ -69,7 +71,7 @@ public val Quaternion.vector: DoubleVector3D
 /**
  * Rotate a vector in a [Float64Space3D] with [quaternion]
  */
-public fun Float64Space3D.rotate(vector: DoubleVector3D, quaternion: Quaternion): DoubleVector3D =
+public fun Float64Space3D.rotate(vector: Float64Vector3D, quaternion: Quaternion): Float64Vector3D =
     with(QuaternionAlgebra) {
         val p = vector.asQuaternion()
         (quaternion * p * quaternion.reciprocal).vector
@@ -80,15 +82,15 @@ public fun Float64Space3D.rotate(vector: DoubleVector3D, quaternion: Quaternion)
  */
 @UnstableKMathAPI
 public fun Float64Space3D.rotate(
-    vector: DoubleVector3D,
+    vector: Float64Vector3D,
     composition: QuaternionAlgebra.() -> Quaternion,
-): DoubleVector3D =
+): Float64Vector3D =
     rotate(vector, QuaternionAlgebra.composition())
 
 /**
  * Rotate a [Float64] vector in 3D space with a rotation matrix
  */
-public fun Float64Space3D.rotate(vector: DoubleVector3D, matrix: Matrix<Double>): DoubleVector3D {
+public fun Float64Space3D.rotate(vector: Float64Vector3D, matrix: Matrix<Double>):  Vector3D<Float64> {
     require(matrix.colNum == 3 && matrix.rowNum == 3) { "Square 3x3 rotation matrix is required" }
     return with(linearSpace) { (matrix dot vector).asVector3D() }
 }
@@ -242,6 +244,8 @@ public fun Quaternion.Companion.fromEuler(
  * A vector consisting of angles
  */
 public data class AngleVector(override val x: Angle, override val y: Angle, override val z: Angle) : Vector3D<Angle> {
+    override val type: SafeType<Angle> get() = Angle.type
+
     public companion object
 }
 

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/testUtils.kt

@@ -5,8 +5,7 @@
 
 package space.kscience.kmath.geometry
 
-import space.kscience.kmath.geometry.euclidean2d.DoubleVector2D
-import space.kscience.kmath.geometry.euclidean3d.DoubleVector3D
+import space.kscience.kmath.structures.Float64
 import kotlin.math.abs
 import kotlin.test.assertEquals
 
@@ -27,12 +26,12 @@ fun grid(
     return xs.flatMap { x -> ys.map { y -> x to y } }
 }
 
-fun assertVectorEquals(expected: DoubleVector2D, actual: DoubleVector2D, absoluteTolerance: Double = 1e-3) {
+fun assertVectorEquals(expected:  Vector2D<Float64>, actual:  Vector2D<Float64>, absoluteTolerance: Double = 1e-3) {
     assertEquals(expected.x, actual.x, absoluteTolerance)
     assertEquals(expected.y, actual.y, absoluteTolerance)
 }
 
-fun assertVectorEquals(expected: DoubleVector3D, actual: DoubleVector3D, absoluteTolerance: Double = 1e-6) {
+fun assertVectorEquals(expected: Vector3D<Float64>, actual:  Vector3D<Float64>, absoluteTolerance: Double = 1e-6) {
     assertEquals(expected.x, actual.x, absoluteTolerance)
     assertEquals(expected.y, actual.y, absoluteTolerance)
     assertEquals(expected.z, actual.z, absoluteTolerance)

kmath-histograms/README.md

@@ -6,19 +6,8 @@
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-histograms:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-histograms:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-histograms:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -27,6 +16,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-histograms:0.4.0-dev-1")
+    implementation("space.kscience:kmath-histograms:0.4.0-dev-3")
 }
 ```

kmath-histograms/build.gradle.kts

@@ -6,6 +6,8 @@ kscience{
     jvm()
     js()
     native()
+    wasm()
+    useCoroutines()
 }
 
 //apply(plugin = "kotlinx-atomicfu")
@@ -21,7 +23,6 @@ kotlin.sourceSets {
         dependencies {
             implementation(project(":kmath-for-real"))
             implementation(projects.kmath.kmathStat)
-            implementation("org.jetbrains.kotlinx:kotlinx-coroutines-test:1.6.0")
         }
     }
 }

kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/Histogram.kt

@@ -47,7 +47,7 @@ public interface Histogram<in T : Any, out V, out B : Bin<T, V>> {
     }
 }
 
-public interface HistogramBuilder<in T : Any, V : Any> {
+public interface HistogramBuilder<in T, V> {
 
     /**
      * The default value increment for a bin
@@ -61,9 +61,9 @@ public interface HistogramBuilder<in T : Any, V : Any> {
 
 }
 
-public fun <T : Any> HistogramBuilder<T, *>.put(point: Point<out T>): Unit = putValue(point)
+public fun <T> HistogramBuilder<T, *>.put(point: Point<T>): Unit = putValue(point)
 
-public fun <T : Any> HistogramBuilder<T, *>.put(vararg point: T): Unit = put(point.asBuffer())
+public inline fun <reified T> HistogramBuilder<T, *>.put(vararg point: T): Unit = put(point.asList().asBuffer())
 
 public fun HistogramBuilder<Double, *>.put(vararg point: Number): Unit =
     put(Float64Buffer(point.map { it.toDouble() }.toDoubleArray()))

kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/UniformHistogram1D.kt

@@ -12,6 +12,7 @@ import space.kscience.kmath.operations.Ring
 import space.kscience.kmath.operations.ScaleOperations
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.MutableBufferFactory
 import kotlin.math.floor
 
 @OptIn(UnstableKMathAPI::class)
@@ -46,6 +47,8 @@ public class UniformHistogram1DGroup<V : Any, A>(
     public val startPoint: Double = 0.0,
 ) : Group<Histogram1D<Double, V>>, ScaleOperations<Histogram1D<Double, V>> where A : Ring<V>, A : ScaleOperations<V> {
 
+    override val bufferFactory: MutableBufferFactory<Histogram1D<Double, V>> = MutableBufferFactory()
+
     override val zero: UniformHistogram1D<V> = UniformHistogram1D(this, emptyMap())
 
     /**

kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/UniformHistogramGroupND.kt

@@ -37,6 +37,8 @@ public class UniformHistogramGroupND<V : Any, A : Field<V>>(
         require(!lower.indices.any { upper[it] - lower[it] < 0 }) { "Range for one of axis is not strictly positive" }
     }
 
+    override val bufferFactory: MutableBufferFactory<HistogramND<Double, HyperSquareDomain, V>> = MutableBufferFactory()
+
     public val dimension: Int get() = lower.size
 
     override val shape: ShapeND = ShapeND(IntArray(binNums.size) { binNums[it] + 2 })
@@ -87,7 +89,7 @@ public class UniformHistogramGroupND<V : Any, A : Field<V>>(
         builder: HistogramBuilder<Double, V>.() -> Unit,
     ): HistogramND<Double, HyperSquareDomain, V> {
         val ndCounter: BufferND<ObjectCounter<V>> =
-            StructureND.buffered(shape) { Counter.of(valueAlgebraND.elementAlgebra) }
+            BufferND(shape) { Counter.of(valueAlgebraND.elementAlgebra) }
         val hBuilder = object : HistogramBuilder<Double, V> {
             override val defaultValue: V get() = valueAlgebraND.elementAlgebra.one
 
@@ -97,7 +99,8 @@ public class UniformHistogramGroupND<V : Any, A : Field<V>>(
             }
         }
         hBuilder.apply(builder)
-        val values: BufferND<V> = BufferND(ndCounter.indices, ndCounter.buffer.mapToBuffer(valueBufferFactory) { it.value })
+        val values: BufferND<V> =
+            BufferND(ndCounter.indices, ndCounter.buffer.mapToBuffer(valueBufferFactory) { it.value })
 
         return HistogramND(this, values)
     }
@@ -128,8 +131,7 @@ public fun <V : Any, A : Field<V>> Histogram.Companion.uniformNDFromRanges(
 
 public fun Histogram.Companion.uniformDoubleNDFromRanges(
     vararg ranges: ClosedFloatingPointRange<Double>,
-): UniformHistogramGroupND<Double, Float64Field> =
-    uniformNDFromRanges(Floa64FieldOpsND, *ranges, bufferFactory = ::Float64Buffer)
+): UniformHistogramGroupND<Double, Float64Field> = uniformNDFromRanges(Floa64FieldOpsND, *ranges)
 
 
 /**
@@ -147,21 +149,18 @@ public fun <V : Any, A : Field<V>> Histogram.Companion.uniformNDFromRanges(
     bufferFactory: BufferFactory<V> = valueAlgebraND.elementAlgebra.bufferFactory,
 ): UniformHistogramGroupND<V, A> = UniformHistogramGroupND(
     valueAlgebraND,
-    DoubleBuffer(
-        ranges
-            .map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
-            .map(ClosedFloatingPointRange<Double>::start)
-    ),
-    DoubleBuffer(
-        ranges
-            .map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
-            .map(ClosedFloatingPointRange<Double>::endInclusive)
-    ),
+    ranges
+        .map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
+        .map(ClosedFloatingPointRange<Double>::start)
+        .asBuffer(),
+    ranges
+        .map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
+        .map(ClosedFloatingPointRange<Double>::endInclusive)
+        .asBuffer(),
     ranges.map(Pair<ClosedFloatingPointRange<Double>, Int>::second).toIntArray(),
     valueBufferFactory = bufferFactory
 )
 
 public fun Histogram.Companion.uniformDoubleNDFromRanges(
     vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>,
-): UniformHistogramGroupND<Double, Float64Field> =
-    uniformNDFromRanges(Floa64FieldOpsND, *ranges, bufferFactory = ::Float64Buffer)
+): UniformHistogramGroupND<Double, Float64Field> = uniformNDFromRanges(Floa64FieldOpsND, *ranges)

kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/TreeHistogramGroup.kt

@@ -14,10 +14,7 @@ import space.kscience.kmath.misc.sorted
 import space.kscience.kmath.operations.Group
 import space.kscience.kmath.operations.Ring
 import space.kscience.kmath.operations.ScaleOperations
-import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.first
-import space.kscience.kmath.structures.indices
-import space.kscience.kmath.structures.last
+import space.kscience.kmath.structures.*
 import java.util.*
 
 private fun <B : ClosedRange<Double>> TreeMap<Double, B>.getBin(value: Double): B? {
@@ -53,6 +50,8 @@ public class TreeHistogramGroup<V : Any, A>(
     @PublishedApi internal val binFactory: (Double) -> DoubleDomain1D,
 ) : Group<TreeHistogram<V>>, ScaleOperations<TreeHistogram<V>> where A : Ring<V>, A : ScaleOperations<V> {
 
+    override val bufferFactory: MutableBufferFactory<TreeHistogram<V>> = MutableBufferFactory()
+
     internal inner class DomainCounter(val domain: DoubleDomain1D, val counter: Counter<V> = Counter.of(valueAlgebra)) :
         ClosedRange<Double> by domain.range
 

kmath-jafama/README.md

@@ -7,19 +7,8 @@ Integration with [Jafama](https://github.com/jeffhain/jafama).
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-jafama:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-jafama:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-jafama:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -28,7 +17,7 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-jafama:0.4.0-dev-1")
+    implementation("space.kscience:kmath-jafama:0.4.0-dev-3")
 }
 ```
 

kmath-jafama/api/kmath-jafama.api

@@ -0,0 +1,128 @@
+public final class space/kscience/kmath/jafama/JafamaDoubleField : space/kscience/kmath/operations/ExtendedField, space/kscience/kmath/operations/Norm, space/kscience/kmath/operations/ScaleOperations {
+	public static final field INSTANCE Lspace/kscience/kmath/jafama/JafamaDoubleField;
+	public fun acos (D)Ljava/lang/Double;
+	public synthetic fun acos (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun acosh (D)Ljava/lang/Double;
+	public synthetic fun acosh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun add (DD)Ljava/lang/Double;
+	public synthetic fun add (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun asin (D)Ljava/lang/Double;
+	public synthetic fun asin (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun asinh (D)Ljava/lang/Double;
+	public synthetic fun asinh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun atan (D)Ljava/lang/Double;
+	public synthetic fun atan (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun atanh (D)Ljava/lang/Double;
+	public synthetic fun atanh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun binaryOperationFunction (Ljava/lang/String;)Lkotlin/jvm/functions/Function2;
+	public fun cos (D)Ljava/lang/Double;
+	public synthetic fun cos (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun cosh (D)Ljava/lang/Double;
+	public synthetic fun cosh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun div (DD)Ljava/lang/Double;
+	public synthetic fun div (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun divide (DD)Ljava/lang/Double;
+	public synthetic fun divide (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun exp (D)Ljava/lang/Double;
+	public synthetic fun exp (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun getBufferFactory ()Lspace/kscience/kmath/structures/MutableBufferFactory;
+	public fun getOne ()Ljava/lang/Double;
+	public synthetic fun getOne ()Ljava/lang/Object;
+	public fun getZero ()Ljava/lang/Double;
+	public synthetic fun getZero ()Ljava/lang/Object;
+	public fun ln (D)Ljava/lang/Double;
+	public synthetic fun ln (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun minus (DD)Ljava/lang/Double;
+	public synthetic fun minus (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun multiply (DD)Ljava/lang/Double;
+	public synthetic fun multiply (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun norm (D)Ljava/lang/Double;
+	public synthetic fun norm (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun number (Ljava/lang/Number;)Ljava/lang/Double;
+	public synthetic fun number (Ljava/lang/Number;)Ljava/lang/Object;
+	public fun plus (DD)Ljava/lang/Double;
+	public synthetic fun plus (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun power (DLjava/lang/Number;)Ljava/lang/Double;
+	public synthetic fun power (Ljava/lang/Object;Ljava/lang/Number;)Ljava/lang/Object;
+	public fun scale (DD)Ljava/lang/Double;
+	public synthetic fun scale (Ljava/lang/Object;D)Ljava/lang/Object;
+	public fun sin (D)Ljava/lang/Double;
+	public synthetic fun sin (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun sinh (D)Ljava/lang/Double;
+	public synthetic fun sinh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun sqrt (D)Ljava/lang/Double;
+	public synthetic fun sqrt (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun tan (D)Ljava/lang/Double;
+	public synthetic fun tan (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun tanh (D)Ljava/lang/Double;
+	public synthetic fun tanh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun times (DD)Ljava/lang/Double;
+	public synthetic fun times (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun unaryMinus (D)Ljava/lang/Double;
+	public synthetic fun unaryMinus (Ljava/lang/Object;)Ljava/lang/Object;
+}
+
+public final class space/kscience/kmath/jafama/StrictJafamaDoubleField : space/kscience/kmath/operations/ExtendedField, space/kscience/kmath/operations/Norm, space/kscience/kmath/operations/ScaleOperations {
+	public static final field INSTANCE Lspace/kscience/kmath/jafama/StrictJafamaDoubleField;
+	public fun acos (D)Ljava/lang/Double;
+	public synthetic fun acos (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun acosh (D)Ljava/lang/Double;
+	public synthetic fun acosh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun add (DD)Ljava/lang/Double;
+	public synthetic fun add (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun asin (D)Ljava/lang/Double;
+	public synthetic fun asin (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun asinh (D)Ljava/lang/Double;
+	public synthetic fun asinh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun atan (D)Ljava/lang/Double;
+	public synthetic fun atan (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun atanh (D)Ljava/lang/Double;
+	public synthetic fun atanh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun binaryOperationFunction (Ljava/lang/String;)Lkotlin/jvm/functions/Function2;
+	public fun cos (D)Ljava/lang/Double;
+	public synthetic fun cos (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun cosh (D)Ljava/lang/Double;
+	public synthetic fun cosh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun div (DD)Ljava/lang/Double;
+	public synthetic fun div (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun divide (DD)Ljava/lang/Double;
+	public synthetic fun divide (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun exp (D)Ljava/lang/Double;
+	public synthetic fun exp (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun getBufferFactory ()Lspace/kscience/kmath/structures/MutableBufferFactory;
+	public fun getOne ()Ljava/lang/Double;
+	public synthetic fun getOne ()Ljava/lang/Object;
+	public fun getZero ()Ljava/lang/Double;
+	public synthetic fun getZero ()Ljava/lang/Object;
+	public fun ln (D)Ljava/lang/Double;
+	public synthetic fun ln (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun minus (DD)Ljava/lang/Double;
+	public synthetic fun minus (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun multiply (DD)Ljava/lang/Double;
+	public synthetic fun multiply (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun norm (D)Ljava/lang/Double;
+	public synthetic fun norm (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun number (Ljava/lang/Number;)Ljava/lang/Double;
+	public synthetic fun number (Ljava/lang/Number;)Ljava/lang/Object;
+	public fun plus (DD)Ljava/lang/Double;
+	public synthetic fun plus (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun power (DLjava/lang/Number;)Ljava/lang/Double;
+	public synthetic fun power (Ljava/lang/Object;Ljava/lang/Number;)Ljava/lang/Object;
+	public fun scale (DD)Ljava/lang/Double;
+	public synthetic fun scale (Ljava/lang/Object;D)Ljava/lang/Object;
+	public fun sin (D)Ljava/lang/Double;
+	public synthetic fun sin (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun sinh (D)Ljava/lang/Double;
+	public synthetic fun sinh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun sqrt (D)Ljava/lang/Double;
+	public synthetic fun sqrt (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun tan (D)Ljava/lang/Double;
+	public synthetic fun tan (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun tanh (D)Ljava/lang/Double;
+	public synthetic fun tanh (Ljava/lang/Object;)Ljava/lang/Object;
+	public fun times (DD)Ljava/lang/Double;
+	public synthetic fun times (Ljava/lang/Object;Ljava/lang/Object;)Ljava/lang/Object;
+	public fun unaryMinus (D)Ljava/lang/Double;
+	public synthetic fun unaryMinus (Ljava/lang/Object;)Ljava/lang/Object;
+}
+

kmath-jupyter/README.md

@@ -6,19 +6,8 @@
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-jupyter:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-jupyter:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-jupyter:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -27,6 +16,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-jupyter:0.4.0-dev-1")
+    implementation("space.kscience:kmath-jupyter:0.4.0-dev-3")
 }
 ```

kmath-kotlingrad/README.md

@@ -8,19 +8,8 @@
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-kotlingrad:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-kotlingrad:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-kotlingrad:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -29,6 +18,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-kotlingrad:0.4.0-dev-1")
+    implementation("space.kscience:kmath-kotlingrad:0.4.0-dev-3")
 }
 ```

kmath-kotlingrad/src/main/kotlin/space/kscience/kmath/kotlingrad/KotlingradExpression.kt

@@ -7,6 +7,7 @@ package space.kscience.kmath.kotlingrad
 
 import ai.hypergraph.kotlingrad.api.SFun
 import ai.hypergraph.kotlingrad.api.SVar
+import space.kscience.attributes.SafeType
 import space.kscience.kmath.expressions.*
 import space.kscience.kmath.operations.NumericAlgebra
 
@@ -25,6 +26,8 @@ public class KotlingradExpression<T : Number, A : NumericAlgebra<T>>(
     public val algebra: A,
     public val mst: MST,
 ) : SpecialDifferentiableExpression<T, KotlingradExpression<T, A>> {
+    override val type: SafeType<T> = algebra.type
+
     override fun invoke(arguments: Map<Symbol, T>): T = mst.interpret(algebra, arguments)
 
     override fun derivativeOrNull(

kmath-memory/README.md

@@ -6,19 +6,8 @@
 
 ## Artifact:
 
-The Maven coordinates of this project are `space.kscience:kmath-memory:0.4.0-dev-1`.
+The Maven coordinates of this project are `space.kscience:kmath-memory:0.4.0-dev-3`.
 
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-memory:0.4.0-dev-1'
-}
-```
 **Gradle Kotlin DSL:**
 ```kotlin
 repositories {
@@ -27,6 +16,6 @@ repositories {
 }
 
 dependencies {
-    implementation("space.kscience:kmath-memory:0.4.0-dev-1")
+    implementation("space.kscience:kmath-memory:0.4.0-dev-3")
 }
 ```