diff --git a/CHANGELOG.md b/CHANGELOG.md
index 259c56c34..e542d210c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -16,17 +16,22 @@
 - ND4J support module submitting `NDStructure` and `NDAlgebra` over `INDArray`.
 - Coroutine-deterministic Monte-Carlo scope with a random number generator.
 - Some minor utilities to `kmath-for-real`.
+- Generic operation result parameter to `MatrixContext`
 
 ### Changed
 - Package changed from `scientifik` to `kscience.kmath`.
-- Gradle version: 6.6 -> 6.7
+- Gradle version: 6.6 -> 6.7.1
 - Minor exceptions refactor (throwing `IllegalArgumentException` by argument checks instead of `IllegalStateException`)
 - `Polynomial` secondary constructor made function.
-- Kotlin version: 1.3.72 -> 1.4.20-M1
+- Kotlin version: 1.3.72 -> 1.4.20
 - `kmath-ast` doesn't depend on heavy `kotlin-reflect` library.
 - Full autodiff refactoring based on `Symbol`
 - `kmath-prob` renamed to `kmath-stat`
 - Grid generators moved to `kmath-for-real`
+- Use `Point<Double>` instead of specialized type in `kmath-for-real`
+- Optimized dot product for buffer matrices moved to `kmath-for-real`
+- EjmlMatrix context is an object
+- Matrix LUP `inverse` renamed to `inverseWithLUP`
 
 ### Deprecated
 
@@ -34,6 +39,7 @@
 - `kmath-koma` module because it doesn't support Kotlin 1.4.
 - Support of `legacy` JS backend (we will support only IR)
 - `toGrid` method.
+- Public visibility of `BufferAccessor2D`
 
 ### Fixed
 - `symbol` method in `MstExtendedField` (https://github.com/mipt-npm/kmath/pull/140)
diff --git a/README.md b/README.md
index c4e3e5374..50a916d2c 100644
--- a/README.md
+++ b/README.md
@@ -132,9 +132,17 @@ submit a feature request if you want something to be implemented first.
 <hr/>
 
 * ### [kmath-for-real](kmath-for-real)
-> 
+> Extension module that should be used to achieve numpy-like behavior. 
+All operations are specialized to work with `Double` numbers without declaring algebraic contexts.
+One can still use generic algebras though.
 >
 > **Maturity**: EXPERIMENTAL
+>
+> **Features:**
+> - [RealVector](kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt) : Numpy-like operations for Buffers/Points
+> - [RealMatrix](kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt) : Numpy-like operations for 2d real structures
+> - [grids](kmath-for-real/src/commonMain/kotlin/kscience/kmath/structures/grids.kt) : Uniform grid generators
+
 <hr/>
 
 * ### [kmath-functions](kmath-functions)
@@ -155,6 +163,12 @@ submit a feature request if you want something to be implemented first.
 > **Maturity**: EXPERIMENTAL
 <hr/>
 
+* ### [kmath-kotlingrad](kmath-kotlingrad)
+> 
+>
+> **Maturity**: EXPERIMENTAL
+<hr/>
+
 * ### [kmath-memory](kmath-memory)
 > 
 >
@@ -167,7 +181,7 @@ submit a feature request if you want something to be implemented first.
 > **Maturity**: EXPERIMENTAL
 >
 > **Features:**
-> - [nd4jarraystrucure](kmath-nd4j/src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt) : NDStructure wrapper for INDArray
+> - [nd4jarraystructure](kmath-nd4j/src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt) : NDStructure wrapper for INDArray
 > - [nd4jarrayrings](kmath-nd4j/src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt) : Rings over Nd4jArrayStructure of Int and Long
 > - [nd4jarrayfields](kmath-nd4j/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : Fields over Nd4jArrayStructure of Float and Double
 
@@ -211,20 +225,12 @@ Release artifacts are accessible from bintray with following configuration (see
 
 ```kotlin
 repositories {
-    jcenter()
-    maven("https://clojars.org/repo")
-    maven("https://dl.bintray.com/egor-bogomolov/astminer/")
-    maven("https://dl.bintray.com/hotkeytlt/maven")
-    maven("https://dl.bintray.com/kotlin/kotlin-eap")
-    maven("https://dl.bintray.com/kotlin/kotlinx")
     maven("https://dl.bintray.com/mipt-npm/kscience")
-    maven("https://jitpack.io")
-    mavenCentral()
 }
 
 dependencies {
-    api("kscience.kmath:kmath-core:0.2.0-dev-3")
-    // api("kscience.kmath:kmath-core-jvm:0.2.0-dev-3") for jvm-specific version
+    api("kscience.kmath:kmath-core:0.2.0-dev-4")
+    // api("kscience.kmath:kmath-core-jvm:0.2.0-dev-4") for jvm-specific version
 }
 ```
 
@@ -236,15 +242,7 @@ Development builds are uploaded to the separate repository:
 
 ```kotlin
 repositories {
-    jcenter()
-    maven("https://clojars.org/repo")
-    maven("https://dl.bintray.com/egor-bogomolov/astminer/")
-    maven("https://dl.bintray.com/hotkeytlt/maven")
-    maven("https://dl.bintray.com/kotlin/kotlin-eap")
-    maven("https://dl.bintray.com/kotlin/kotlinx")
     maven("https://dl.bintray.com/mipt-npm/dev")
-    maven("https://jitpack.io")
-    mavenCentral()
 }
 ```
 
diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts
index d42627ff0..c079eaa84 100644
--- a/examples/build.gradle.kts
+++ b/examples/build.gradle.kts
@@ -35,6 +35,9 @@ dependencies {
     implementation(project(":kmath-dimensions"))
     implementation(project(":kmath-ejml"))
     implementation(project(":kmath-nd4j"))
+
+    implementation(project(":kmath-for-real"))
+
     implementation("org.deeplearning4j:deeplearning4j-core:1.0.0-beta7")
     implementation("org.nd4j:nd4j-native:1.0.0-beta7")
 
@@ -51,7 +54,11 @@ dependencies {
     implementation("org.jetbrains.kotlinx:kotlinx-io:0.2.0-npm-dev-11")
     implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-20")
     implementation("org.slf4j:slf4j-simple:1.7.30")
-    "benchmarksImplementation"("org.jetbrains.kotlinx:kotlinx.benchmark.runtime-jvm:0.2.0-dev-8")
+
+    // plotting
+    implementation("kscience.plotlykt:plotlykt-server:0.3.1-dev")
+
+    "benchmarksImplementation"("org.jetbrains.kotlinx:kotlinx.benchmark.runtime-jvm:0.2.0-dev-20")
     "benchmarksImplementation"(sourceSets.main.get().output + sourceSets.main.get().runtimeClasspath)
 }
 
@@ -62,7 +69,7 @@ benchmark {
     // This one matches sourceSet name above
 
     configurations.register("fast") {
-        warmups = 5 // number of warmup iterations
+        warmups = 1 // number of warmup iterations
         iterations = 3 // number of iterations
         iterationTime = 500 // time in seconds per iteration
         iterationTimeUnit = "ms" // time unity for iterationTime, default is seconds
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/structures/ArrayBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ArrayBenchmark.kt
similarity index 96%
rename from examples/src/benchmarks/kotlin/kscience/kmath/structures/ArrayBenchmark.kt
rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ArrayBenchmark.kt
index f7e0d05c6..ebf31a590 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/structures/ArrayBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ArrayBenchmark.kt
@@ -1,4 +1,4 @@
-package kscience.kmath.structures
+package kscience.kmath.benchmarks
 
 import org.openjdk.jmh.annotations.Benchmark
 import org.openjdk.jmh.annotations.Scope
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/structures/BufferBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/BufferBenchmark.kt
similarity index 85%
rename from examples/src/benchmarks/kotlin/kscience/kmath/structures/BufferBenchmark.kt
rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/BufferBenchmark.kt
index 8b6fd4a51..4c64517f1 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/structures/BufferBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/BufferBenchmark.kt
@@ -1,7 +1,9 @@
-package kscience.kmath.structures
+package kscience.kmath.benchmarks
 
 import kscience.kmath.operations.Complex
 import kscience.kmath.operations.complex
+import kscience.kmath.structures.MutableBuffer
+import kscience.kmath.structures.RealBuffer
 import org.openjdk.jmh.annotations.Benchmark
 import org.openjdk.jmh.annotations.Scope
 import org.openjdk.jmh.annotations.State
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
new file mode 100644
index 000000000..ec8714617
--- /dev/null
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
@@ -0,0 +1,50 @@
+package kscience.kmath.linear
+
+
+import kotlinx.benchmark.Benchmark
+import kscience.kmath.commons.linear.CMMatrixContext
+import kscience.kmath.commons.linear.CMMatrixContext.dot
+import kscience.kmath.commons.linear.inverse
+import kscience.kmath.commons.linear.toCM
+import kscience.kmath.ejml.EjmlMatrixContext
+import kscience.kmath.ejml.inverse
+import kscience.kmath.ejml.toEjml
+import kscience.kmath.operations.invoke
+import kscience.kmath.structures.Matrix
+import org.openjdk.jmh.annotations.Scope
+import org.openjdk.jmh.annotations.State
+import kotlin.random.Random
+
+@State(Scope.Benchmark)
+class LinearAlgebraBenchmark {
+    companion object {
+        val random = Random(1224)
+        val dim = 100
+
+        //creating invertible matrix
+        val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+        val l = Matrix.real(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 }
+        val matrix = l dot u
+    }
+
+    @Benchmark
+    fun kmathLUPInversion() {
+        MatrixContext.real.inverseWithLUP(matrix)
+    }
+
+    @Benchmark
+    fun cmLUPInversion() {
+        CMMatrixContext {
+            val cm = matrix.toCM()             //avoid overhead on conversion
+            inverse(cm)
+        }
+    }
+
+    @Benchmark
+    fun ejmlInverse() {
+        EjmlMatrixContext {
+            val km = matrix.toEjml()      //avoid overhead on conversion
+            inverse(km)
+        }
+    }
+}
\ No newline at end of file
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt
new file mode 100644
index 000000000..9d2b02245
--- /dev/null
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/MultiplicationBenchmark.kt
@@ -0,0 +1,60 @@
+package kscience.kmath.benchmarks
+
+import kotlinx.benchmark.Benchmark
+import kscience.kmath.commons.linear.CMMatrixContext
+import kscience.kmath.commons.linear.CMMatrixContext.dot
+import kscience.kmath.commons.linear.toCM
+import kscience.kmath.ejml.EjmlMatrixContext
+import kscience.kmath.ejml.toEjml
+import kscience.kmath.linear.real
+import kscience.kmath.operations.invoke
+import kscience.kmath.structures.Matrix
+import org.openjdk.jmh.annotations.Scope
+import org.openjdk.jmh.annotations.State
+import kotlin.random.Random
+
+@State(Scope.Benchmark)
+class MultiplicationBenchmark {
+    companion object {
+        val random = Random(12224)
+        val dim = 1000
+
+        //creating invertible matrix
+        val matrix1 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+        val matrix2 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
+
+        val cmMatrix1 = matrix1.toCM()
+        val cmMatrix2 = matrix2.toCM()
+
+        val ejmlMatrix1 = matrix1.toEjml()
+        val ejmlMatrix2 = matrix2.toEjml()
+    }
+
+    @Benchmark
+    fun commonsMathMultiplication() {
+        CMMatrixContext.invoke {
+            cmMatrix1 dot cmMatrix2
+        }
+    }
+
+    @Benchmark
+    fun ejmlMultiplication() {
+        EjmlMatrixContext.invoke {
+            ejmlMatrix1 dot ejmlMatrix2
+        }
+    }
+
+    @Benchmark
+    fun ejmlMultiplicationwithConversion() {
+        val ejmlMatrix1 = matrix1.toEjml()
+        val ejmlMatrix2 = matrix2.toEjml()
+        EjmlMatrixContext.invoke {
+            ejmlMatrix1 dot ejmlMatrix2
+        }
+    }
+
+    @Benchmark
+    fun bufferedMultiplication() {
+        matrix1 dot matrix2
+    }
+}
\ No newline at end of file
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/structures/NDFieldBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt
similarity index 94%
rename from examples/src/benchmarks/kotlin/kscience/kmath/structures/NDFieldBenchmark.kt
rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt
index 8ec47ae81..1be8e7236 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/structures/NDFieldBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt
@@ -1,7 +1,8 @@
-package kscience.kmath.structures
+package kscience.kmath.benchmarks
 
 import kscience.kmath.operations.RealField
 import kscience.kmath.operations.invoke
+import kscience.kmath.structures.*
 import org.openjdk.jmh.annotations.Benchmark
 import org.openjdk.jmh.annotations.Scope
 import org.openjdk.jmh.annotations.State
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/structures/ViktorBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt
similarity index 90%
rename from examples/src/benchmarks/kotlin/kscience/kmath/structures/ViktorBenchmark.kt
rename to examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt
index df35b0b78..8663e353c 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/structures/ViktorBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt
@@ -1,7 +1,10 @@
-package kscience.kmath.structures
+package kscience.kmath.benchmarks
 
 import kscience.kmath.operations.RealField
 import kscience.kmath.operations.invoke
+import kscience.kmath.structures.BufferedNDField
+import kscience.kmath.structures.NDField
+import kscience.kmath.structures.RealNDField
 import kscience.kmath.viktor.ViktorNDField
 import org.jetbrains.bio.viktor.F64Array
 import org.openjdk.jmh.annotations.Benchmark
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/utils/utils.kt b/examples/src/benchmarks/kotlin/kscience/kmath/utils/utils.kt
deleted file mode 100644
index 329b7b17b..000000000
--- a/examples/src/benchmarks/kotlin/kscience/kmath/utils/utils.kt
+++ /dev/null
@@ -1,11 +0,0 @@
-package kscience.kmath.utils
-
-import kotlin.contracts.InvocationKind
-import kotlin.contracts.contract
-import kotlin.system.measureTimeMillis
-
-internal inline fun measureAndPrint(title: String, block: () -> Unit) {
-    contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
-    val time = measureTimeMillis(block)
-    println("$title completed in $time millis")
-}
diff --git a/examples/src/main/kotlin/kscience/kmath/commons/fit/fitWithAutoDiff.kt b/examples/src/main/kotlin/kscience/kmath/commons/fit/fitWithAutoDiff.kt
new file mode 100644
index 000000000..c0cd9dc5c
--- /dev/null
+++ b/examples/src/main/kotlin/kscience/kmath/commons/fit/fitWithAutoDiff.kt
@@ -0,0 +1,102 @@
+package kscience.kmath.commons.fit
+
+import kotlinx.html.br
+import kotlinx.html.h3
+import kscience.kmath.commons.optimization.chiSquared
+import kscience.kmath.commons.optimization.minimize
+import kscience.kmath.expressions.symbol
+import kscience.kmath.real.RealVector
+import kscience.kmath.real.map
+import kscience.kmath.real.step
+import kscience.kmath.stat.*
+import kscience.kmath.structures.asIterable
+import kscience.kmath.structures.toList
+import kscience.plotly.*
+import kscience.plotly.models.ScatterMode
+import kscience.plotly.models.TraceValues
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+//Forward declaration of symbols that will be used in expressions.
+// This declaration is required for
+private val a by symbol
+private val b by symbol
+private val c by symbol
+
+/**
+ * Shortcut to use buffers in plotly
+ */
+operator fun TraceValues.invoke(vector: RealVector) {
+    numbers = vector.asIterable()
+}
+
+/**
+ * Least squares fie with auto-differentiation. Uses `kmath-commons` and `kmath-for-real` modules.
+ */
+fun main() {
+
+    //A generator for a normally distributed values
+    val generator = Distribution.normal()
+
+    //A chain/flow of random values with the given seed
+    val chain = generator.sample(RandomGenerator.default(112667))
+
+
+    //Create a uniformly distributed x values like numpy.arrange
+    val x = 1.0..100.0 step 1.0
+
+
+    //Perform an operation on each x value (much more effective, than numpy)
+    val y = x.map {
+        val value = it.pow(2) + it + 1
+        value + chain.nextDouble() * sqrt(value)
+    }
+    // this will also work, but less effective:
+    // val y = x.pow(2)+ x + 1 + chain.nextDouble()
+
+    // create same errors for all xs
+    val yErr = y.map { sqrt(it) }//RealVector.same(x.size, sigma)
+
+    // compute differentiable chi^2 sum for given model ax^2 + bx + c
+    val chi2 = Fitting.chiSquared(x, y, yErr) { x1 ->
+        //bind variables to autodiff context
+        val a = bind(a)
+        val b = bind(b)
+        //Include default value for c if it is not provided as a parameter
+        val c = bindOrNull(c) ?: one
+        a * x1.pow(2) + b * x1 + c
+    }
+
+    //minimize the chi^2 in given starting point. Derivatives are not required, they are already included.
+    val result: OptimizationResult<Double> = chi2.minimize(a to 1.5, b to 0.9, c to 1.0)
+
+    //display a page with plot and numerical results
+    val page = Plotly.page {
+        plot {
+            scatter {
+                mode = ScatterMode.markers
+                x(x)
+                y(y)
+                error_y {
+                    array = yErr.toList()
+                }
+                name = "data"
+            }
+            scatter {
+                mode = ScatterMode.lines
+                x(x)
+                y(x.map { result.point[a]!! * it.pow(2) + result.point[b]!! * it + 1 })
+                name = "fit"
+            }
+        }
+        br()
+        h3{
+            +"Fit result: $result"
+        }
+        h3{
+            +"Chi2/dof = ${result.value / (x.size - 3)}"
+        }
+    }
+
+    page.makeFile()
+}
\ No newline at end of file
diff --git a/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt b/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt
deleted file mode 100644
index 3316f3236..000000000
--- a/examples/src/main/kotlin/kscience/kmath/linear/LinearAlgebraBenchmark.kt
+++ /dev/null
@@ -1,50 +0,0 @@
-package kscience.kmath.linear
-
-import kscience.kmath.commons.linear.CMMatrixContext
-import kscience.kmath.commons.linear.inverse
-import kscience.kmath.commons.linear.toCM
-import kscience.kmath.ejml.EjmlMatrixContext
-import kscience.kmath.ejml.inverse
-import kscience.kmath.operations.RealField
-import kscience.kmath.operations.invoke
-import kscience.kmath.structures.Matrix
-import kotlin.random.Random
-import kotlin.system.measureTimeMillis
-
-fun main() {
-    val random = Random(1224)
-    val dim = 100
-    //creating invertible matrix
-    val u = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
-    val l = Matrix.real(dim, dim) { i, j -> if (i >= j) random.nextDouble() else 0.0 }
-    val matrix = l dot u
-
-    val n = 5000 // iterations
-
-    MatrixContext.real {
-        repeat(50) { inverse(matrix) }
-        val inverseTime = measureTimeMillis { repeat(n) { inverse(matrix) } }
-        println("[kmath] Inversion of $n matrices $dim x $dim finished in $inverseTime millis")
-    }
-
-    //commons-math
-
-    val commonsTime = measureTimeMillis {
-        CMMatrixContext {
-            val cm = matrix.toCM()             //avoid overhead on conversion
-            repeat(n) { inverse(cm) }
-        }
-    }
-
-
-    println("[commons-math] Inversion of $n matrices $dim x $dim finished in $commonsTime millis")
-
-    val ejmlTime = measureTimeMillis {
-        (EjmlMatrixContext(RealField)) {
-            val km = matrix.toEjml()      //avoid overhead on conversion
-            repeat(n) { inverse(km) }
-        }
-    }
-
-    println("[ejml] Inversion of $n matrices $dim x $dim finished in $ejmlTime millis")
-}
diff --git a/examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt b/examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt
deleted file mode 100644
index d1011e8f5..000000000
--- a/examples/src/main/kotlin/kscience/kmath/linear/MultiplicationBenchmark.kt
+++ /dev/null
@@ -1,38 +0,0 @@
-package kscience.kmath.linear
-
-import kscience.kmath.commons.linear.CMMatrixContext
-import kscience.kmath.commons.linear.toCM
-import kscience.kmath.ejml.EjmlMatrixContext
-import kscience.kmath.operations.RealField
-import kscience.kmath.operations.invoke
-import kscience.kmath.structures.Matrix
-import kotlin.random.Random
-import kotlin.system.measureTimeMillis
-
-fun main() {
-    val random = Random(12224)
-    val dim = 1000
-    //creating invertible matrix
-    val matrix1 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
-    val matrix2 = Matrix.real(dim, dim) { i, j -> if (i <= j) random.nextDouble() else 0.0 }
-
-//    //warmup
-//    matrix1 dot matrix2
-
-    CMMatrixContext {
-        val cmMatrix1 = matrix1.toCM()
-        val cmMatrix2 = matrix2.toCM()
-        val cmTime = measureTimeMillis { cmMatrix1 dot cmMatrix2 }
-        println("CM implementation time: $cmTime")
-    }
-
-    (EjmlMatrixContext(RealField)) {
-        val ejmlMatrix1 = matrix1.toEjml()
-        val ejmlMatrix2 = matrix2.toEjml()
-        val ejmlTime = measureTimeMillis { ejmlMatrix1 dot ejmlMatrix2 }
-        println("EJML implementation time: $ejmlTime")
-    }
-
-    val genericTime = measureTimeMillis { val res = matrix1 dot matrix2 }
-    println("Generic implementation time: $genericTime")
-}
diff --git a/examples/src/main/kotlin/kscience/kmath/commons/prob/DistributionBenchmark.kt b/examples/src/main/kotlin/kscience/kmath/stat/DistributionBenchmark.kt
similarity index 100%
rename from examples/src/main/kotlin/kscience/kmath/commons/prob/DistributionBenchmark.kt
rename to examples/src/main/kotlin/kscience/kmath/stat/DistributionBenchmark.kt
diff --git a/examples/src/main/kotlin/kscience/kmath/commons/prob/DistributionDemo.kt b/examples/src/main/kotlin/kscience/kmath/stat/DistributionDemo.kt
similarity index 83%
rename from examples/src/main/kotlin/kscience/kmath/commons/prob/DistributionDemo.kt
rename to examples/src/main/kotlin/kscience/kmath/stat/DistributionDemo.kt
index 6146e17af..24a4cb1a7 100644
--- a/examples/src/main/kotlin/kscience/kmath/commons/prob/DistributionDemo.kt
+++ b/examples/src/main/kotlin/kscience/kmath/stat/DistributionDemo.kt
@@ -1,14 +1,17 @@
-package kscience.kmath.commons.prob
+package kscience.kmath.stat
 
 import kotlinx.coroutines.runBlocking
 import kscience.kmath.chains.Chain
 import kscience.kmath.chains.collectWithState
-import kscience.kmath.stat.Distribution
-import kscience.kmath.stat.RandomGenerator
-import kscience.kmath.stat.normal
 
+/**
+ * The state of distribution averager
+ */
 private data class AveragingChainState(var num: Int = 0, var value: Double = 0.0)
 
+/**
+ * Averaging
+ */
 private fun Chain<Double>.mean(): Chain<Double> = collectWithState(AveragingChainState(), { it.copy() }) { chain ->
     val next = chain.next()
     num++
diff --git a/kmath-ast/build.gradle.kts b/kmath-ast/build.gradle.kts
index 2f95d427d..3e3c0475f 100644
--- a/kmath-ast/build.gradle.kts
+++ b/kmath-ast/build.gradle.kts
@@ -17,3 +17,7 @@ kotlin.sourceSets {
         }
     }
 }
+
+readme{
+    maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE
+}
\ No newline at end of file
diff --git a/kmath-ast/src/jvmMain/kotlin/kscience/kmath/ast/parser.kt b/kmath-ast/src/jvmMain/kotlin/kscience/kmath/ast/parser.kt
index 94cd3b321..79cf77a6a 100644
--- a/kmath-ast/src/jvmMain/kotlin/kscience/kmath/ast/parser.kt
+++ b/kmath-ast/src/jvmMain/kotlin/kscience/kmath/ast/parser.kt
@@ -83,11 +83,10 @@ public object ArithmeticsEvaluator : Grammar<MST>() {
 }
 
 /**
- * Tries to parse the string into [MST].
+ * Tries to parse the string into [MST]. Returns [ParseResult] representing expression or error.
  *
  * @receiver the string to parse.
  * @return the [MST] node.
- * @author Alexander Nozik
  */
 public fun String.tryParseMath(): ParseResult<MST> = ArithmeticsEvaluator.tryParseToEnd(this)
 
@@ -96,6 +95,5 @@ public fun String.tryParseMath(): ParseResult<MST> = ArithmeticsEvaluator.tryPar
  *
  * @receiver the string to parse.
  * @return the [MST] node.
- * @author Alexander Nozik
  */
 public fun String.parseMath(): MST = ArithmeticsEvaluator.parseToEnd(this)
diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
index 377d167bc..712927400 100644
--- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
+++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
@@ -5,8 +5,7 @@ import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.NDStructure
 import org.apache.commons.math3.linear.*
 
-public class CMMatrix(public val origin: RealMatrix, features: Set<MatrixFeature>? = null) :
-    FeaturedMatrix<Double> {
+public class CMMatrix(public val origin: RealMatrix, features: Set<MatrixFeature>? = null) : FeaturedMatrix<Double> {
     public override val rowNum: Int get() = origin.rowDimension
     public override val colNum: Int get() = origin.columnDimension
 
@@ -55,7 +54,7 @@ public fun Point<Double>.toCM(): CMVector = if (this is CMVector) this else {
 
 public fun RealVector.toPoint(): CMVector = CMVector(this)
 
-public object CMMatrixContext : MatrixContext<Double> {
+public object CMMatrixContext : MatrixContext<Double, CMMatrix> {
     public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): CMMatrix {
         val array = Array(rows) { i -> DoubleArray(columns) { j -> initializer(i, j) } }
         return CMMatrix(Array2DRowRealMatrix(array))
@@ -79,7 +78,7 @@ public object CMMatrixContext : MatrixContext<Double> {
     public override fun multiply(a: Matrix<Double>, k: Number): CMMatrix =
         CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble()))
 
-    public override operator fun Matrix<Double>.times(value: Double): Matrix<Double> =
+    public override operator fun Matrix<Double>.times(value: Double): CMMatrix =
         produce(rowNum, colNum) { i, j -> get(i, j) * value }
 }
 
diff --git a/kmath-core/README.md b/kmath-core/README.md
index 42a513a10..3a919e85a 100644
--- a/kmath-core/README.md
+++ b/kmath-core/README.md
@@ -12,7 +12,7 @@ The core features of KMath:
 
 > #### Artifact:
 >
-> This module artifact: `kscience.kmath:kmath-core:0.2.0-dev-3`.
+> This module artifact: `kscience.kmath:kmath-core:0.2.0-dev-4`.
 >
 > Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-core/_latestVersion)
 >
@@ -30,7 +30,7 @@ The core features of KMath:
 > }
 > 
 > dependencies {
->     implementation 'kscience.kmath:kmath-core:0.2.0-dev-3'
+>     implementation 'kscience.kmath:kmath-core:0.2.0-dev-4'
 > }
 > ```
 > **Gradle Kotlin DSL:**
@@ -44,6 +44,6 @@ The core features of KMath:
 > }
 > 
 > dependencies {
->     implementation("kscience.kmath:kmath-core:0.2.0-dev-3")
+>     implementation("kscience.kmath:kmath-core:0.2.0-dev-4")
 > }
 > ```
diff --git a/kmath-core/build.gradle.kts b/kmath-core/build.gradle.kts
index 7f889d9b4..9ed7e690b 100644
--- a/kmath-core/build.gradle.kts
+++ b/kmath-core/build.gradle.kts
@@ -1,5 +1,3 @@
-import ru.mipt.npm.gradle.Maturity
-
 plugins {
     id("ru.mipt.npm.mpp")
     id("ru.mipt.npm.native")
@@ -13,7 +11,7 @@ kotlin.sourceSets.commonMain {
 
 readme {
     description = "Core classes, algebra definitions, basic linear algebra"
-    maturity = Maturity.DEVELOPMENT
+    maturity = ru.mipt.npm.gradle.Maturity.DEVELOPMENT
     propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md"))
 
     feature(
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
index d51f40890..8b50bbe33 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
@@ -9,8 +9,8 @@ import kscience.kmath.structures.*
  */
 public class BufferMatrixContext<T : Any, R : Ring<T>>(
     public override val elementContext: R,
-    private val bufferFactory: BufferFactory<T>
-) : GenericMatrixContext<T, R> {
+    private val bufferFactory: BufferFactory<T>,
+) : GenericMatrixContext<T, R, BufferMatrix<T>> {
     public override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): BufferMatrix<T> {
         val buffer = bufferFactory(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
         return BufferMatrix(rows, columns, buffer)
@@ -22,15 +22,15 @@ public class BufferMatrixContext<T : Any, R : Ring<T>>(
 }
 
 @Suppress("OVERRIDE_BY_INLINE")
-public object RealMatrixContext : GenericMatrixContext<Double, RealField> {
+public object RealMatrixContext : GenericMatrixContext<Double, RealField, BufferMatrix<Double>> {
     public override val elementContext: RealField
         get() = RealField
 
     public override inline fun produce(
         rows: Int,
         columns: Int,
-        initializer: (i: Int, j: Int) -> Double
-    ): Matrix<Double> {
+        initializer: (i: Int, j: Int) -> Double,
+    ): BufferMatrix<Double> {
         val buffer = RealBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
         return BufferMatrix(rows, columns, buffer)
     }
@@ -43,15 +43,15 @@ public class BufferMatrix<T : Any>(
     public override val rowNum: Int,
     public override val colNum: Int,
     public val buffer: Buffer<out T>,
-    public override val features: Set<MatrixFeature> = emptySet()
+    public override val features: Set<MatrixFeature> = emptySet(),
 ) : FeaturedMatrix<T> {
-    override val shape: IntArray
-        get() = intArrayOf(rowNum, colNum)
 
     init {
         require(buffer.size == rowNum * colNum) { "Dimension mismatch for matrix structure" }
     }
 
+    override val shape: IntArray get() = intArrayOf(rowNum, colNum)
+
     public override fun suggestFeature(vararg features: MatrixFeature): BufferMatrix<T> =
         BufferMatrix(rowNum, colNum, buffer, this.features + features)
 
@@ -86,28 +86,3 @@ public class BufferMatrix<T : Any>(
         else "Matrix(rowsNum = $rowNum, colNum = $colNum, features=$features)"
     }
 }
-
-/**
- * Optimized dot product for real matrices
- */
-public infix fun BufferMatrix<Double>.dot(other: BufferMatrix<Double>): BufferMatrix<Double> {
-    require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
-    val array = DoubleArray(this.rowNum * other.colNum)
-
-    //convert to array to insure there is not memory indirection
-    fun Buffer<out Double>.unsafeArray() = if (this is RealBuffer)
-        array
-    else
-        DoubleArray(size) { get(it) }
-
-    val a = this.buffer.unsafeArray()
-    val b = other.buffer.unsafeArray()
-
-    for (i in (0 until rowNum))
-        for (j in (0 until other.colNum))
-            for (k in (0 until colNum))
-                array[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j]
-
-    val buffer = RealBuffer(array)
-    return BufferMatrix(rowNum, other.colNum, buffer)
-}
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt
index 5d9af8608..68272203c 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/FeaturedMatrix.kt
@@ -27,9 +27,8 @@ public interface FeaturedMatrix<T : Any> : Matrix<T> {
 public inline fun Structure2D.Companion.real(
     rows: Int,
     columns: Int,
-    initializer: (Int, Int) -> Double
-): Matrix<Double> =
-    MatrixContext.real.produce(rows, columns, initializer)
+    initializer: (Int, Int) -> Double,
+): BufferMatrix<Double> = MatrixContext.real.produce(rows, columns, initializer)
 
 /**
  * Build a square matrix from given elements.
@@ -58,7 +57,7 @@ public inline fun <reified T : Any> Matrix<*>.getFeature(): T? =
 /**
  * Diagonal matrix of ones. The matrix is virtual no actual matrix is created
  */
-public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.one(rows: Int, columns: Int): FeaturedMatrix<T> =
+public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R, *>.one(rows: Int, columns: Int): FeaturedMatrix<T> =
     VirtualMatrix(rows, columns, DiagonalFeature) { i, j ->
         if (i == j) elementContext.one else elementContext.zero
     }
@@ -67,7 +66,7 @@ public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.one(rows: Int, colu
 /**
  * A virtual matrix of zeroes
  */
-public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R>.zero(rows: Int, columns: Int): FeaturedMatrix<T> =
+public fun <T : Any, R : Ring<T>> GenericMatrixContext<T, R, *>.zero(rows: Int, columns: Int): FeaturedMatrix<T> =
     VirtualMatrix(rows, columns) { _, _ -> elementContext.zero }
 
 public class TransposedFeature<T : Any>(public val original: Matrix<T>) : MatrixFeature
@@ -81,6 +80,4 @@ public fun <T : Any> Matrix<T>.transpose(): Matrix<T> {
         rowNum,
         setOf(TransposedFeature(this))
     ) { i, j -> get(j, i) }
-}
-
-public infix fun Matrix<Double>.dot(other: Matrix<Double>): Matrix<Double> = with(MatrixContext.real) { dot(other) }
+}
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt
index bb80bcafd..099fa1909 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LUPDecomposition.kt
@@ -1,25 +1,18 @@
 package kscience.kmath.linear
 
-import kscience.kmath.operations.Field
-import kscience.kmath.operations.RealField
-import kscience.kmath.operations.Ring
-import kscience.kmath.operations.invoke
-import kscience.kmath.structures.BufferAccessor2D
-import kscience.kmath.structures.Matrix
-import kscience.kmath.structures.Structure2D
-import kotlin.reflect.KClass
+import kscience.kmath.operations.*
+import kscience.kmath.structures.*
 
 /**
  * Common implementation of [LUPDecompositionFeature]
  */
 public class LUPDecomposition<T : Any>(
-    public val context: GenericMatrixContext<T, out Field<T>>,
+    public val context: MatrixContext<T, FeaturedMatrix<T>>,
+    public val elementContext: Field<T>,
     public val lu: Structure2D<T>,
     public val pivot: IntArray,
-    private val even: Boolean
+    private val even: Boolean,
 ) : LUPDecompositionFeature<T>, DeterminantFeature<T> {
-    public val elementContext: Field<T>
-        get() = context.elementContext
 
     /**
      * Returns the matrix L of the decomposition.
@@ -64,23 +57,25 @@ public class LUPDecomposition<T : Any>(
 
 }
 
-public fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.abs(value: T): T =
+@PublishedApi
+internal fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, *>.abs(value: T): T =
     if (value > elementContext.zero) value else elementContext { -value }
 
 /**
- * Create a lup decomposition of generic matrix
+ * Create a lup decomposition of generic matrix.
  */
-public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
-    type: KClass<T>,
+public fun <T : Comparable<T>> MatrixContext<T, FeaturedMatrix<T>>.lup(
+    factory: MutableBufferFactory<T>,
+    elementContext: Field<T>,
     matrix: Matrix<T>,
-    checkSingular: (T) -> Boolean
+    checkSingular: (T) -> Boolean,
 ): LUPDecomposition<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 KEEP-176
-    BufferAccessor2D(type, matrix.rowNum, matrix.colNum).run {
+    BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {
         elementContext {
             val lu = create(matrix)
 
@@ -112,14 +107,14 @@ public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.l
                     luRow[col] = sum
 
                     // maintain best permutation choice
-                    if (this@lup.abs(sum) > largest) {
-                        largest = this@lup.abs(sum)
+                    if (abs(sum) > largest) {
+                        largest = abs(sum)
                         max = row
                     }
                 }
 
                 // Singularity check
-                check(!checkSingular(this@lup.abs(lu[max, col]))) { "The matrix is singular" }
+                check(!checkSingular(abs(lu[max, col]))) { "The matrix is singular" }
 
                 // Pivot if necessary
                 if (max != col) {
@@ -143,23 +138,23 @@ public inline fun <T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.l
                 for (row in col + 1 until m) lu[row, col] /= luDiag
             }
 
-            return LUPDecomposition(this@lup, lu.collect(), pivot, even)
+            return LUPDecomposition(this@lup, elementContext, lu.collect(), pivot, even)
         }
     }
 }
 
-public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.lup(
+public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.lup(
     matrix: Matrix<T>,
-    checkSingular: (T) -> Boolean
-): LUPDecomposition<T> = lup(T::class, matrix, checkSingular)
+    noinline checkSingular: (T) -> Boolean,
+): LUPDecomposition<T> = lup(MutableBuffer.Companion::auto, elementContext, matrix, checkSingular)
 
-public fun GenericMatrixContext<Double, RealField>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> =
-    lup(Double::class, matrix) { it < 1e-11 }
+public fun MatrixContext<Double, FeaturedMatrix<Double>>.lup(matrix: Matrix<Double>): LUPDecomposition<Double> =
+    lup(Buffer.Companion::real, RealField, matrix) { it < 1e-11 }
 
-public fun <T : Any> LUPDecomposition<T>.solve(type: KClass<T>, matrix: Matrix<T>): Matrix<T> {
+public fun <T : Any> LUPDecomposition<T>.solveWithLUP(factory: MutableBufferFactory<T>, matrix: Matrix<T>): FeaturedMatrix<T> {
     require(matrix.rowNum == pivot.size) { "Matrix dimension mismatch. Expected ${pivot.size}, but got ${matrix.colNum}" }
 
-    BufferAccessor2D(type, matrix.rowNum, matrix.colNum).run {
+    BufferAccessor2D(matrix.rowNum, matrix.colNum, factory).run {
         elementContext {
             // Apply permutations to b
             val bp = create { _, _ -> zero }
@@ -201,27 +196,34 @@ public fun <T : Any> LUPDecomposition<T>.solve(type: KClass<T>, matrix: Matrix<T
     }
 }
 
-public inline fun <reified T : Any> LUPDecomposition<T>.solve(matrix: Matrix<T>): Matrix<T> = solve(T::class, matrix)
+public inline fun <reified T : Any> LUPDecomposition<T>.solveWithLUP(matrix: Matrix<T>): Matrix<T> =
+    solveWithLUP(MutableBuffer.Companion::auto, matrix)
 
 /**
- * Solve a linear equation **a*x = b**
+ * Solve a linear equation **a*x = b** using LUP decomposition
  */
-public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.solve(
+public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.solveWithLUP(
     a: Matrix<T>,
     b: Matrix<T>,
-    checkSingular: (T) -> Boolean
-): Matrix<T> {
+    noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
+    noinline checkSingular: (T) -> Boolean,
+): FeaturedMatrix<T> {
     // Use existing decomposition if it is provided by matrix
-    val decomposition = a.getFeature() ?: lup(T::class, a, checkSingular)
-    return decomposition.solve(T::class, b)
+    val decomposition = a.getFeature() ?: lup(bufferFactory, elementContext, a, checkSingular)
+    return decomposition.solveWithLUP(bufferFactory, b)
 }
 
-public fun RealMatrixContext.solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> = solve(a, b) { it < 1e-11 }
+public fun RealMatrixContext.solveWithLUP(a: Matrix<Double>, b: Matrix<Double>): FeaturedMatrix<Double> =
+    solveWithLUP(a, b) { it < 1e-11 }
 
-public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F>.inverse(
+public inline fun <reified T : Comparable<T>, F : Field<T>> GenericMatrixContext<T, F, FeaturedMatrix<T>>.inverseWithLUP(
     matrix: Matrix<T>,
-    checkSingular: (T) -> Boolean
-): Matrix<T> = solve(matrix, one(matrix.rowNum, matrix.colNum), checkSingular)
+    noinline bufferFactory: MutableBufferFactory<T> = MutableBuffer.Companion::auto,
+    noinline checkSingular: (T) -> Boolean,
+): FeaturedMatrix<T> = solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), bufferFactory, checkSingular)
 
-public fun RealMatrixContext.inverse(matrix: Matrix<Double>): Matrix<Double> =
-    solve(matrix, one(matrix.rowNum, matrix.colNum)) { it < 1e-11 }
+/**
+ * Inverses a square matrix using LUP decomposition. Non square matrix will throw a error.
+ */
+public fun RealMatrixContext.inverseWithLUP(matrix: Matrix<Double>): FeaturedMatrix<Double> =
+    solveWithLUP(matrix, one(matrix.rowNum, matrix.colNum), Buffer.Companion::real) { it < 1e-11 }
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
index 174f62c12..912c3f950 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
@@ -12,16 +12,17 @@ import kscience.kmath.structures.asSequence
 /**
  * Basic operations on matrices. Operates on [Matrix]
  */
-public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
+public interface MatrixContext<T : Any, out M : Matrix<T>> : SpaceOperations<Matrix<T>> {
     /**
      * Produce a matrix with this context and given dimensions
      */
-    public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T>
+    public fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): M
 
-    public override fun binaryOperation(operation: String): (left: Matrix<T>, right: Matrix<T>) -> Matrix<T> =
+    @Suppress("UNCHECKED_CAST")
+    public override fun binaryOperation(operation: String): (left: Matrix<T>, right: Matrix<T>) -> M =
         when (operation) {
             "dot" -> { left, right -> left dot right }
-            else -> super.binaryOperation(operation)
+            else -> super.binaryOperation(operation) as (Matrix<T>, Matrix<T>) -> M
         }
 
     /**
@@ -31,7 +32,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
      * @param other the multiplier.
      * @return the dot product.
      */
-    public infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T>
+    public infix fun Matrix<T>.dot(other: Matrix<T>): M
 
     /**
      * Computes the dot product of this matrix and a vector.
@@ -49,7 +50,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
      * @param value the multiplier.
      * @receiver the product.
      */
-    public operator fun Matrix<T>.times(value: T): Matrix<T>
+    public operator fun Matrix<T>.times(value: T): M
 
     /**
      * Multiplies an element by a matrix of it.
@@ -58,7 +59,7 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
      * @param value the multiplier.
      * @receiver the product.
      */
-    public operator fun T.times(m: Matrix<T>): Matrix<T> = m * this
+    public operator fun T.times(m: Matrix<T>): M = m * this
 
     public companion object {
         /**
@@ -71,18 +72,18 @@ public interface MatrixContext<T : Any> : SpaceOperations<Matrix<T>> {
          */
         public fun <T : Any, R : Ring<T>> buffered(
             ring: R,
-            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
-        ): GenericMatrixContext<T, R> = BufferMatrixContext(ring, bufferFactory)
+            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
+        ): GenericMatrixContext<T, R, BufferMatrix<T>> = BufferMatrixContext(ring, bufferFactory)
 
         /**
          * Automatic buffered matrix, unboxed if it is possible
          */
-        public inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R> =
+        public inline fun <reified T : Any, R : Ring<T>> auto(ring: R): GenericMatrixContext<T, R, BufferMatrix<T>> =
             buffered(ring, Buffer.Companion::auto)
     }
 }
 
-public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
+public interface GenericMatrixContext<T : Any, R : Ring<T>, out M : Matrix<T>> : MatrixContext<T, M> {
     /**
      * The ring context for matrix elements
      */
@@ -93,7 +94,7 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
      */
     public fun point(size: Int, initializer: (Int) -> T): Point<T>
 
-    public override infix fun Matrix<T>.dot(other: Matrix<T>): Matrix<T> {
+    public override infix fun Matrix<T>.dot(other: Matrix<T>): M {
         //TODO add typed error
         require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
 
@@ -114,10 +115,10 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
         }
     }
 
-    public override operator fun Matrix<T>.unaryMinus(): Matrix<T> =
+    public override operator fun Matrix<T>.unaryMinus(): M =
         produce(rowNum, colNum) { i, j -> elementContext { -get(i, j) } }
 
-    public override fun add(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
+    public override fun add(a: Matrix<T>, b: Matrix<T>): M {
         require(a.rowNum == b.rowNum && a.colNum == b.colNum) {
             "Matrix operation dimension mismatch. [${a.rowNum},${a.colNum}] + [${b.rowNum},${b.colNum}]"
         }
@@ -125,7 +126,7 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
         return produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] + b[i, j] } }
     }
 
-    public override operator fun Matrix<T>.minus(b: Matrix<T>): Matrix<T> {
+    public override operator fun Matrix<T>.minus(b: Matrix<T>): M {
         require(rowNum == b.rowNum && colNum == b.colNum) {
             "Matrix operation dimension mismatch. [$rowNum,$colNum] - [${b.rowNum},${b.colNum}]"
         }
@@ -133,11 +134,11 @@ public interface GenericMatrixContext<T : Any, R : Ring<T>> : MatrixContext<T> {
         return produce(rowNum, colNum) { i, j -> elementContext { get(i, j) + b[i, j] } }
     }
 
-    public override fun multiply(a: Matrix<T>, k: Number): Matrix<T> =
+    public override fun multiply(a: Matrix<T>, k: Number): M =
         produce(a.rowNum, a.colNum) { i, j -> elementContext { a[i, j] * k } }
 
-    public operator fun Number.times(matrix: FeaturedMatrix<T>): Matrix<T> = matrix * this
+    public operator fun Number.times(matrix: FeaturedMatrix<T>): M = multiply(matrix, this)
 
-    public override operator fun Matrix<T>.times(value: T): Matrix<T> =
+    public override operator fun Matrix<T>.times(value: T): M =
         produce(rowNum, colNum) { i, j -> elementContext { get(i, j) * value } }
 }
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/VectorSpace.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/VectorSpace.kt
index 67056c6b2..2a3b8f5d1 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/VectorSpace.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/VectorSpace.kt
@@ -15,7 +15,7 @@ public interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
     public val space: S
     override val zero: Point<T> get() = produce { space.zero }
 
-    public fun produce(initializer: (Int) -> T): Point<T>
+    public fun produce(initializer: S.(Int) -> T): Point<T>
 
     /**
      * Produce a space-element of this vector space for expressions
@@ -48,7 +48,7 @@ public interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
         public fun <T : Any, S : Space<T>> buffered(
             size: Int,
             space: S,
-            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
+            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
         ): BufferVectorSpace<T, S> = BufferVectorSpace(size, space, bufferFactory)
 
         /**
@@ -63,8 +63,8 @@ public interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
 public class BufferVectorSpace<T : Any, S : Space<T>>(
     override val size: Int,
     override val space: S,
-    public val bufferFactory: BufferFactory<T>
+    public val bufferFactory: BufferFactory<T>,
 ) : VectorSpace<T, S> {
-    override fun produce(initializer: (Int) -> T): Buffer<T> = bufferFactory(size, initializer)
+    override fun produce(initializer: S.(Int) -> T): Buffer<T> = bufferFactory(size) { space.initializer(it) }
     //override fun produceElement(initializer: (Int) -> T): Vector<T, S> = BufferVector(this, produce(initializer))
 }
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/misc/annotations.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/misc/annotations.kt
new file mode 100644
index 000000000..d70ac7b39
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/misc/annotations.kt
@@ -0,0 +1,4 @@
+package kscience.kmath.misc
+
+@RequiresOptIn("This API is unstable and could change in future", RequiresOptIn.Level.WARNING)
+public annotation class UnstableKMathAPI
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt
index 657da6b4e..4527a2a42 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraExtensions.kt
@@ -38,6 +38,11 @@ public fun <T> Space<T>.average(data: Iterable<T>): T = sum(data) / data.count()
  */
 public fun <T> Space<T>.average(data: Sequence<T>): T = sum(data) / data.count()
 
+/**
+ * Absolute of the comparable [value]
+ */
+public fun <T : Comparable<T>> Space<T>.abs(value: T): T = if (value > zero) value else -value
+
 /**
  * Returns the sum of all elements in the iterable in provided space.
  *
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt
index 00fc4aef0..5d7ba611f 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt
@@ -1,25 +1,31 @@
 package kscience.kmath.structures
 
-import kotlin.reflect.KClass
-
 /**
  * A context that allows to operate on a [MutableBuffer] as on 2d array
  */
-public class BufferAccessor2D<T : Any>(public val type: KClass<T>, public val rowNum: Int, public val colNum: Int) {
+internal class BufferAccessor2D<T : Any>(
+    public val rowNum: Int,
+    public val colNum: Int,
+    val factory: MutableBufferFactory<T>,
+) {
     public operator fun Buffer<T>.get(i: Int, j: Int): T = get(i + colNum * j)
 
     public operator fun MutableBuffer<T>.set(i: Int, j: Int, value: T) {
         set(i + colNum * j, value)
     }
 
-    public inline fun create(init: (i: Int, j: Int) -> T): MutableBuffer<T> =
-        MutableBuffer.auto(type, rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
+    public inline fun create(crossinline init: (i: Int, j: Int) -> T): MutableBuffer<T> =
+        factory(rowNum * colNum) { offset -> init(offset / colNum, offset % colNum) }
 
     public fun create(mat: Structure2D<T>): MutableBuffer<T> = create { i, j -> mat[i, j] }
 
     //TODO optimize wrapper
-    public fun MutableBuffer<T>.collect(): Structure2D<T> =
-        NDStructure.auto(type, rowNum, colNum) { (i, j) -> get(i, j) }.as2D()
+    public fun MutableBuffer<T>.collect(): Structure2D<T> = NDStructure.build(
+        DefaultStrides(intArrayOf(rowNum, colNum)),
+        factory
+    ) { (i, j) ->
+        get(i, j)
+    }.as2D()
 
     public inner class Row(public val buffer: MutableBuffer<T>, public val rowIndex: Int) : MutableBuffer<T> {
         override val size: Int get() = colNum
@@ -30,7 +36,7 @@ public class BufferAccessor2D<T : Any>(public val type: KClass<T>, public val ro
             buffer[rowIndex, index] = value
         }
 
-        override fun copy(): MutableBuffer<T> = MutableBuffer.auto(type, colNum) { get(it) }
+        override fun copy(): MutableBuffer<T> = factory(colNum) { get(it) }
         override operator fun iterator(): Iterator<T> = (0 until colNum).map(::get).iterator()
 
     }
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt
index 5174eb314..bfec6f871 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt
@@ -102,6 +102,11 @@ public fun <T> Buffer<T>.asSequence(): Sequence<T> = Sequence(::iterator)
  */
 public fun <T> Buffer<T>.asIterable(): Iterable<T> = Iterable(::iterator)
 
+/**
+ * Converts this [Buffer] to a new [List]
+ */
+public fun <T> Buffer<T>.toList(): List<T> = asSequence().toList()
+
 /**
  * Returns an [IntRange] of the valid indices for this [Buffer].
  */
diff --git a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt
index 7cfa25a66..0a582e339 100644
--- a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt
+++ b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt
@@ -1,5 +1,6 @@
 package kscience.kmath.linear
 
+import kscience.kmath.operations.invoke
 import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.NDStructure
 import kscience.kmath.structures.as2D
@@ -38,7 +39,7 @@ class MatrixTest {
         infix fun Matrix<Double>.pow(power: Int): Matrix<Double> {
             var res = this
             repeat(power - 1) {
-                res = res dot this
+                res = RealMatrixContext.invoke { res dot this@pow }
             }
             return res
         }
diff --git a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt
index de07d3639..28dfe46ec 100644
--- a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt
+++ b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt
@@ -9,7 +9,7 @@ class RealLUSolverTest {
     @Test
     fun testInvertOne() {
         val matrix = MatrixContext.real.one(2, 2)
-        val inverted = MatrixContext.real.inverse(matrix)
+        val inverted = MatrixContext.real.inverseWithLUP(matrix)
         assertEquals(matrix, inverted)
     }
 
@@ -37,7 +37,7 @@ class RealLUSolverTest {
             1.0, 3.0
         )
 
-        val inverted = MatrixContext.real.inverse(matrix)
+        val inverted = MatrixContext.real.inverseWithLUP(matrix)
 
         val expected = Matrix.square(
             0.375, -0.125,
diff --git a/kmath-coroutines/src/commonMain/kotlin/kscience/kmath/chains/Chain.kt b/kmath-coroutines/src/commonMain/kotlin/kscience/kmath/chains/Chain.kt
index 8c15e52c7..7ff7b7aae 100644
--- a/kmath-coroutines/src/commonMain/kotlin/kscience/kmath/chains/Chain.kt
+++ b/kmath-coroutines/src/commonMain/kotlin/kscience/kmath/chains/Chain.kt
@@ -47,7 +47,7 @@ public fun <T> Iterator<T>.asChain(): Chain<T> = SimpleChain { next() }
 public fun <T> Sequence<T>.asChain(): Chain<T> = iterator().asChain()
 
 /**
- * A simple chain of independent tokens
+ * A simple chain of independent tokens. [fork] returns the same chain.
  */
 public class SimpleChain<out R>(private val gen: suspend () -> R) : Chain<R> {
     public override suspend fun next(): R = gen()
diff --git a/kmath-dimensions/build.gradle.kts b/kmath-dimensions/build.gradle.kts
index 412d7162f..9bf89fc43 100644
--- a/kmath-dimensions/build.gradle.kts
+++ b/kmath-dimensions/build.gradle.kts
@@ -18,3 +18,7 @@ kotlin.sourceSets {
         }
     }
 }
+
+readme{
+    maturity = ru.mipt.npm.gradle.Maturity.PROTOTYPE
+}
diff --git a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
index 4111cb78d..68a5dc262 100644
--- a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
+++ b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
@@ -42,7 +42,7 @@ public interface DMatrix<T, R : Dimension, C : Dimension> : Structure2D<T> {
  * An inline wrapper for a Matrix
  */
 public inline class DMatrixWrapper<T, R : Dimension, C : Dimension>(
-    public val structure: Structure2D<T>
+    private val structure: Structure2D<T>
 ) : DMatrix<T, R, C> {
     override val shape: IntArray get() = structure.shape
     override operator fun get(i: Int, j: Int): T = structure[i, j]
@@ -81,7 +81,7 @@ public inline class DPointWrapper<T, D : Dimension>(public val point: Point<T>)
 /**
  * Basic operations on dimension-safe matrices. Operates on [Matrix]
  */
-public inline class DMatrixContext<T : Any, Ri : Ring<T>>(public val context: GenericMatrixContext<T, Ri>) {
+public inline class DMatrixContext<T : Any, Ri : Ring<T>>(public val context: GenericMatrixContext<T, Ri, Matrix<T>>) {
     public inline fun <reified R : Dimension, reified C : Dimension> Matrix<T>.coerce(): DMatrix<T, R, C> {
         require(rowNum == Dimension.dim<R>().toInt()) {
             "Row number mismatch: expected ${Dimension.dim<R>()} but found $rowNum"
diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
index 52826a7b1..31792e39c 100644
--- a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
+++ b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
@@ -1,23 +1,22 @@
 package kscience.kmath.ejml
 
-import org.ejml.simple.SimpleMatrix
 import kscience.kmath.linear.MatrixContext
 import kscience.kmath.linear.Point
-import kscience.kmath.operations.Space
-import kscience.kmath.operations.invoke
 import kscience.kmath.structures.Matrix
+import org.ejml.simple.SimpleMatrix
+
+/**
+ * Converts this matrix to EJML one.
+ */
+public fun Matrix<Double>.toEjml(): EjmlMatrix =
+    if (this is EjmlMatrix) this else EjmlMatrixContext.produce(rowNum, colNum) { i, j -> get(i, j) }
 
 /**
  * Represents context of basic operations operating with [EjmlMatrix].
  *
  * @author Iaroslav Postovalov
  */
-public class EjmlMatrixContext(private val space: Space<Double>) : MatrixContext<Double> {
-    /**
-     * Converts this matrix to EJML one.
-     */
-    public fun Matrix<Double>.toEjml(): EjmlMatrix =
-        if (this is EjmlMatrix) this else produce(rowNum, colNum) { i, j -> get(i, j) }
+public object EjmlMatrixContext : MatrixContext<Double, EjmlMatrix> {
 
     /**
      * Converts this vector to EJML one.
@@ -47,11 +46,10 @@ public class EjmlMatrixContext(private val space: Space<Double>) : MatrixContext
         EjmlMatrix(toEjml().origin - b.toEjml().origin)
 
     public override fun multiply(a: Matrix<Double>, k: Number): EjmlMatrix =
-        produce(a.rowNum, a.colNum) { i, j -> space { a[i, j] * k } }
+        produce(a.rowNum, a.colNum) { i, j -> a[i, j] * k.toDouble() }
 
-    public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix = EjmlMatrix(toEjml().origin.scale(value))
-
-    public companion object
+    public override operator fun Matrix<Double>.times(value: Double): EjmlMatrix =
+        EjmlMatrix(toEjml().origin.scale(value))
 }
 
 /**
diff --git a/kmath-for-real/README.md b/kmath-for-real/README.md
new file mode 100644
index 000000000..2ddf78e57
--- /dev/null
+++ b/kmath-for-real/README.md
@@ -0,0 +1,44 @@
+# Real number specialization module (`kmath-for-real`)
+
+ - [RealVector](src/commonMain/kotlin/kscience/kmath/real/RealVector.kt) : Numpy-like operations for Buffers/Points
+ - [RealMatrix](src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt) : Numpy-like operations for 2d real structures
+ - [grids](src/commonMain/kotlin/kscience/kmath/structures/grids.kt) : Uniform grid generators
+
+
+> #### Artifact:
+>
+> This module artifact: `kscience.kmath:kmath-for-real:0.2.0-dev-4`.
+>
+> Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-for-real/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-for-real/_latestVersion)
+>
+> Bintray development version:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-for-real/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-for-real/_latestVersion)
+>
+> **Gradle:**
+>
+> ```gradle
+> repositories {
+>     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" }
+>     maven { url 'https://dl.bintray.com/mipt-npm/kscience' }
+>     maven { url 'https://dl.bintray.com/mipt-npm/dev' }
+>     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+
+> }
+> 
+> dependencies {
+>     implementation 'kscience.kmath:kmath-for-real:0.2.0-dev-4'
+> }
+> ```
+> **Gradle Kotlin DSL:**
+>
+> ```kotlin
+> repositories {
+>     maven("https://dl.bintray.com/kotlin/kotlin-eap")
+>     maven("https://dl.bintray.com/mipt-npm/kscience")
+>     maven("https://dl.bintray.com/mipt-npm/dev")
+>     maven("https://dl.bintray.com/hotkeytlt/maven")
+> }
+> 
+> dependencies {
+>     implementation("kscience.kmath:kmath-for-real:0.2.0-dev-4")
+> }
+> ```
diff --git a/kmath-for-real/build.gradle.kts b/kmath-for-real/build.gradle.kts
index 2a4539c10..f26f98c2c 100644
--- a/kmath-for-real/build.gradle.kts
+++ b/kmath-for-real/build.gradle.kts
@@ -7,3 +7,31 @@ kotlin.sourceSets.commonMain {
         api(project(":kmath-core"))
     }
 }
+
+readme {
+    description = """
+        Extension module that should be used to achieve numpy-like behavior. 
+        All operations are specialized to work with `Double` numbers without declaring algebraic contexts.
+        One can still use generic algebras though.
+        """.trimIndent()
+    maturity = ru.mipt.npm.gradle.Maturity.EXPERIMENTAL
+    propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md"))
+
+    feature(
+        id = "RealVector",
+        description = "Numpy-like operations for Buffers/Points",
+        ref = "src/commonMain/kotlin/kscience/kmath/real/RealVector.kt"
+    )
+
+    feature(
+        id = "RealMatrix",
+        description = "Numpy-like operations for 2d real structures",
+        ref = "src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt"
+    )
+
+    feature(
+        id = "grids",
+        description = "Uniform grid generators",
+        ref = "src/commonMain/kotlin/kscience/kmath/structures/grids.kt"
+    )
+}
diff --git a/kmath-for-real/docs/README-TEMPLATE.md b/kmath-for-real/docs/README-TEMPLATE.md
new file mode 100644
index 000000000..670844bd0
--- /dev/null
+++ b/kmath-for-real/docs/README-TEMPLATE.md
@@ -0,0 +1,5 @@
+# Real number specialization module (`kmath-for-real`)
+
+${features}
+
+${artifact}
diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
index 1860b5870..e8ad835e5 100644
--- a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
+++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
@@ -1,12 +1,14 @@
 package kscience.kmath.real
 
+import kscience.kmath.linear.FeaturedMatrix
 import kscience.kmath.linear.MatrixContext
 import kscience.kmath.linear.RealMatrixContext.elementContext
 import kscience.kmath.linear.VirtualMatrix
+import kscience.kmath.linear.inverseWithLUP
+import kscience.kmath.misc.UnstableKMathAPI
 import kscience.kmath.operations.invoke
 import kscience.kmath.operations.sum
 import kscience.kmath.structures.Buffer
-import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.RealBuffer
 import kscience.kmath.structures.asIterable
 import kotlin.math.pow
@@ -23,7 +25,7 @@ import kotlin.math.pow
  *  Functions that help create a real (Double) matrix
  */
 
-public typealias RealMatrix = Matrix<Double>
+public typealias RealMatrix = FeaturedMatrix<Double>
 
 public fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double): RealMatrix =
     MatrixContext.real.produce(rowNum, colNum, initializer)
@@ -36,7 +38,7 @@ public fun Sequence<DoubleArray>.toMatrix(): RealMatrix = toList().let {
     MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] }
 }
 
-public fun Matrix<Double>.repeatStackVertical(n: Int): RealMatrix =
+public fun RealMatrix.repeatStackVertical(n: Int): RealMatrix =
     VirtualMatrix(rowNum * n, colNum) { row, col ->
         get(if (row == 0) 0 else row % rowNum, col)
     }
@@ -45,76 +47,65 @@ public fun Matrix<Double>.repeatStackVertical(n: Int): RealMatrix =
  *  Operations for matrix and real number
  */
 
-public operator fun Matrix<Double>.times(double: Double): RealMatrix =
+public operator fun RealMatrix.times(double: Double): RealMatrix =
     MatrixContext.real.produce(rowNum, colNum) { row, col ->
         this[row, col] * double
     }
 
-public operator fun Matrix<Double>.plus(double: Double): RealMatrix =
+public operator fun RealMatrix.plus(double: Double): RealMatrix =
     MatrixContext.real.produce(rowNum, colNum) { row, col ->
         this[row, col] + double
     }
 
-public operator fun Matrix<Double>.minus(double: Double): RealMatrix =
+public operator fun RealMatrix.minus(double: Double): RealMatrix =
     MatrixContext.real.produce(rowNum, colNum) { row, col ->
         this[row, col] - double
     }
 
-public operator fun Matrix<Double>.div(double: Double): RealMatrix =
+public operator fun RealMatrix.div(double: Double): RealMatrix =
     MatrixContext.real.produce(rowNum, colNum) { row, col ->
         this[row, col] / double
     }
 
-public operator fun Double.times(matrix: Matrix<Double>): RealMatrix =
+public operator fun Double.times(matrix: RealMatrix): RealMatrix =
     MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
         this * matrix[row, col]
     }
 
-public operator fun Double.plus(matrix: Matrix<Double>): RealMatrix =
+public operator fun Double.plus(matrix: RealMatrix): RealMatrix =
     MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
         this + matrix[row, col]
     }
 
-public operator fun Double.minus(matrix: Matrix<Double>): RealMatrix =
+public operator fun Double.minus(matrix: RealMatrix): RealMatrix =
     MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
         this - matrix[row, col]
     }
 
 // TODO: does this operation make sense? Should it be 'this/matrix[row, col]'?
-//operator fun Double.div(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
+//operator fun Double.div(matrix: RealMatrix) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
 //    row, col -> matrix[row, col] / this
 //}
 
-/*
- *  Per-element (!) square and power operations
- */
-
-public fun Matrix<Double>.square(): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { row, col ->
-    this[row, col].pow(2)
-}
-
-public fun Matrix<Double>.pow(n: Int): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { i, j ->
-    this[i, j].pow(n)
-}
-
 /*
  * Operations on two matrices (per-element!)
  */
 
-public operator fun Matrix<Double>.times(other: Matrix<Double>): RealMatrix =
+@UnstableKMathAPI
+public operator fun RealMatrix.times(other: RealMatrix): RealMatrix =
     MatrixContext.real.produce(rowNum, colNum) { row, col -> this[row, col] * other[row, col] }
 
-public operator fun Matrix<Double>.plus(other: Matrix<Double>): RealMatrix =
+public operator fun RealMatrix.plus(other: RealMatrix): RealMatrix =
     MatrixContext.real.add(this, other)
 
-public operator fun Matrix<Double>.minus(other: Matrix<Double>): RealMatrix =
+public operator fun RealMatrix.minus(other: RealMatrix): RealMatrix =
     MatrixContext.real.produce(rowNum, colNum) { row, col -> this[row, col] - other[row, col] }
 
 /*
  *  Operations on columns
  */
 
-public inline fun Matrix<Double>.appendColumn(crossinline mapper: (Buffer<Double>) -> Double): Matrix<Double> =
+public inline fun RealMatrix.appendColumn(crossinline mapper: (Buffer<Double>) -> Double): RealMatrix =
     MatrixContext.real.produce(rowNum, colNum + 1) { row, col ->
         if (col < colNum)
             this[row, col]
@@ -122,28 +113,28 @@ public inline fun Matrix<Double>.appendColumn(crossinline mapper: (Buffer<Double
             mapper(rows[row])
     }
 
-public fun Matrix<Double>.extractColumns(columnRange: IntRange): RealMatrix =
+public fun RealMatrix.extractColumns(columnRange: IntRange): RealMatrix =
     MatrixContext.real.produce(rowNum, columnRange.count()) { row, col ->
         this[row, columnRange.first + col]
     }
 
-public fun Matrix<Double>.extractColumn(columnIndex: Int): RealMatrix =
+public fun RealMatrix.extractColumn(columnIndex: Int): RealMatrix =
     extractColumns(columnIndex..columnIndex)
 
-public fun Matrix<Double>.sumByColumn(): RealBuffer = RealBuffer(colNum) { j ->
+public fun RealMatrix.sumByColumn(): RealBuffer = RealBuffer(colNum) { j ->
     val column = columns[j]
     elementContext { sum(column.asIterable()) }
 }
 
-public fun Matrix<Double>.minByColumn(): RealBuffer = RealBuffer(colNum) { j ->
+public fun RealMatrix.minByColumn(): RealBuffer = RealBuffer(colNum) { j ->
     columns[j].asIterable().minOrNull() ?: error("Cannot produce min on empty column")
 }
 
-public fun Matrix<Double>.maxByColumn(): RealBuffer = RealBuffer(colNum) { j ->
+public fun RealMatrix.maxByColumn(): RealBuffer = RealBuffer(colNum) { j ->
     columns[j].asIterable().maxOrNull() ?: error("Cannot produce min on empty column")
 }
 
-public fun Matrix<Double>.averageByColumn(): RealBuffer = RealBuffer(colNum) { j ->
+public fun RealMatrix.averageByColumn(): RealBuffer = RealBuffer(colNum) { j ->
     columns[j].asIterable().average()
 }
 
@@ -151,7 +142,39 @@ public fun Matrix<Double>.averageByColumn(): RealBuffer = RealBuffer(colNum) { j
  * Operations processing all elements
  */
 
-public fun Matrix<Double>.sum(): Double = elements().map { (_, value) -> value }.sum()
-public fun Matrix<Double>.min(): Double? = elements().map { (_, value) -> value }.minOrNull()
-public fun Matrix<Double>.max(): Double? = elements().map { (_, value) -> value }.maxOrNull()
-public fun Matrix<Double>.average(): Double = elements().map { (_, value) -> value }.average()
+public fun RealMatrix.sum(): Double = elements().map { (_, value) -> value }.sum()
+public fun RealMatrix.min(): Double? = elements().map { (_, value) -> value }.minOrNull()
+public fun RealMatrix.max(): Double? = elements().map { (_, value) -> value }.maxOrNull()
+public fun RealMatrix.average(): Double = elements().map { (_, value) -> value }.average()
+
+public inline fun RealMatrix.map(transform: (Double) -> Double): RealMatrix =
+    MatrixContext.real.produce(rowNum, colNum) { i, j ->
+        transform(get(i, j))
+    }
+
+/**
+ * Inverse a square real matrix using LUP decomposition
+ */
+public fun RealMatrix.inverseWithLUP(): RealMatrix = MatrixContext.real.inverseWithLUP(this)
+
+//extended operations
+
+public fun RealMatrix.pow(p: Double): RealMatrix = map { it.pow(p) }
+
+public fun RealMatrix.pow(p: Int): RealMatrix = map { it.pow(p) }
+
+public fun exp(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.exp(it) }
+
+public fun sqrt(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.sqrt(it) }
+
+public fun RealMatrix.square(): RealMatrix = map { it.pow(2) }
+
+public fun sin(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.sin(it) }
+
+public fun cos(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.cos(it) }
+
+public fun tan(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.tan(it) }
+
+public fun ln(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.ln(it) }
+
+public fun log10(arg: RealMatrix): RealMatrix = arg.map { kotlin.math.log10(it) }
\ No newline at end of file
diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt
index 1ee33ee32..596692782 100644
--- a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt
+++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealVector.kt
@@ -1,47 +1,82 @@
 package kscience.kmath.real
 
-import kscience.kmath.linear.BufferVectorSpace
 import kscience.kmath.linear.Point
-import kscience.kmath.linear.VectorSpace
 import kscience.kmath.operations.Norm
-import kscience.kmath.operations.RealField
-import kscience.kmath.operations.SpaceElement
 import kscience.kmath.structures.Buffer
-import kscience.kmath.structures.RealBuffer
 import kscience.kmath.structures.asBuffer
 import kscience.kmath.structures.asIterable
+import kotlin.math.pow
 import kotlin.math.sqrt
 
-public typealias RealPoint = Point<Double>
-
-public fun RealPoint.asVector(): RealVector = RealVector(this)
-public fun DoubleArray.asVector(): RealVector = asBuffer().asVector()
-public fun List<Double>.asVector(): RealVector = asBuffer().asVector()
+public typealias RealVector = Point<Double>
 
 public object VectorL2Norm : Norm<Point<out Number>, Double> {
     override fun norm(arg: Point<out Number>): Double = sqrt(arg.asIterable().sumByDouble(Number::toDouble))
 }
 
-public inline class RealVector(private val point: Point<Double>) :
-    SpaceElement<RealPoint, RealVector, VectorSpace<Double, RealField>>, RealPoint {
-    public override val size: Int get() = point.size
-    public override val context: VectorSpace<Double, RealField> get() = space(point.size)
+public operator fun Buffer.Companion.invoke(vararg doubles: Double): RealVector = doubles.asBuffer()
 
-    public override fun unwrap(): RealPoint = point
-    public override fun RealPoint.wrap(): RealVector = RealVector(this)
-    public override operator fun get(index: Int): Double = point[index]
-    public override operator fun iterator(): Iterator<Double> = point.iterator()
+/**
+ * Fill the vector of given [size] with given [value]
+ */
+public fun Buffer.Companion.same(size: Int, value: Number): RealVector = real(size) { value.toDouble() }
 
-    public companion object {
-        private val spaceCache: MutableMap<Int, BufferVectorSpace<Double, RealField>> = hashMapOf()
+// Transformation methods
 
-        public inline operator fun invoke(dim: Int, initializer: (Int) -> Double): RealVector =
-            RealVector(RealBuffer(dim, initializer))
+public inline fun RealVector.map(transform: (Double) -> Double): RealVector =
+    Buffer.real(size) { transform(get(it)) }
 
-        public operator fun invoke(vararg values: Double): RealVector = values.asVector()
+public inline fun RealVector.mapIndexed(transform: (index: Int, value: Double) -> Double): RealVector =
+    Buffer.real(size) { transform(it, get(it)) }
 
-        public fun space(dim: Int): BufferVectorSpace<Double, RealField> = spaceCache.getOrPut(dim) {
-            BufferVectorSpace(dim, RealField) { size, init -> Buffer.real(size, init) }
-        }
-    }
-}
+public operator fun RealVector.plus(other: RealVector): RealVector =
+    mapIndexed { index, value -> value + other[index] }
+
+public operator fun RealVector.plus(number: Number): RealVector = map { it + number.toDouble() }
+
+public operator fun Number.plus(vector: RealVector): RealVector = vector + this
+
+public operator fun RealVector.unaryMinus(): Buffer<Double> = map { -it }
+
+public operator fun RealVector.minus(other: RealVector): RealVector =
+    mapIndexed { index, value -> value - other[index] }
+
+public operator fun RealVector.minus(number: Number): RealVector = map { it - number.toDouble() }
+
+public operator fun Number.minus(vector: RealVector): RealVector = vector.map { toDouble() - it }
+
+public operator fun RealVector.times(other: RealVector): RealVector =
+    mapIndexed { index, value -> value * other[index] }
+
+public operator fun RealVector.times(number: Number): RealVector = map { it * number.toDouble() }
+
+public operator fun Number.times(vector: RealVector): RealVector = vector * this
+
+public operator fun RealVector.div(other: RealVector): RealVector =
+    mapIndexed { index, value -> value / other[index] }
+
+public operator fun RealVector.div(number: Number): RealVector = map { it / number.toDouble() }
+
+public operator fun Number.div(vector: RealVector): RealVector = vector.map { toDouble() / it }
+
+//extended operations
+
+public fun RealVector.pow(p: Double): RealVector = map { it.pow(p) }
+
+public fun RealVector.pow(p: Int): RealVector = map { it.pow(p) }
+
+public fun exp(vector: RealVector): RealVector = vector.map { kotlin.math.exp(it) }
+
+public fun sqrt(vector: RealVector): RealVector = vector.map { kotlin.math.sqrt(it) }
+
+public fun RealVector.square(): RealVector = map { it.pow(2) }
+
+public fun sin(vector: RealVector): RealVector = vector.map { kotlin.math.sin(it) }
+
+public fun cos(vector: RealVector): RealVector = vector.map { kotlin.math.cos(it) }
+
+public fun tan(vector: RealVector): RealVector = vector.map { kotlin.math.tan(it) }
+
+public fun ln(vector: RealVector): RealVector = vector.map { kotlin.math.ln(it) }
+
+public fun log10(vector: RealVector): RealVector = vector.map { kotlin.math.log10(it) }
diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/dot.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/dot.kt
new file mode 100644
index 000000000..9beffe6bb
--- /dev/null
+++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/dot.kt
@@ -0,0 +1,31 @@
+package kscience.kmath.real
+
+import kscience.kmath.linear.BufferMatrix
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.RealBuffer
+
+
+/**
+ * Optimized dot product for real matrices
+ */
+public infix fun BufferMatrix<Double>.dot(other: BufferMatrix<Double>): BufferMatrix<Double> {
+    require(colNum == other.rowNum) { "Matrix dot operation dimension mismatch: ($rowNum, $colNum) x (${other.rowNum}, ${other.colNum})" }
+    val resultArray = DoubleArray(this.rowNum * other.colNum)
+
+    //convert to array to insure there is no memory indirection
+    fun Buffer<out Double>.unsafeArray() = if (this is RealBuffer)
+        this.array
+    else
+        DoubleArray(size) { get(it) }
+
+    val a = this.buffer.unsafeArray()
+    val b = other.buffer.unsafeArray()
+
+    for (i in (0 until rowNum))
+        for (j in (0 until other.colNum))
+            for (k in (0 until colNum))
+                resultArray[i * other.colNum + j] += a[i * colNum + k] * b[k * other.colNum + j]
+
+    val buffer = RealBuffer(resultArray)
+    return BufferMatrix(rowNum, other.colNum, buffer)
+}
\ No newline at end of file
diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/grids.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/grids.kt
index bd0e092e0..69a149fb8 100644
--- a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/grids.kt
+++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/grids.kt
@@ -1,6 +1,5 @@
 package kscience.kmath.real
 
-import kscience.kmath.linear.Point
 import kscience.kmath.structures.asBuffer
 import kotlin.math.abs
 
@@ -34,7 +33,7 @@ public fun ClosedFloatingPointRange<Double>.toSequenceWithStep(step: Double): Se
     }
 }
 
-public infix fun ClosedFloatingPointRange<Double>.step(step: Double): Point<Double> =
+public infix fun ClosedFloatingPointRange<Double>.step(step: Double): RealVector =
     toSequenceWithStep(step).toList().asBuffer()
 
 /**
diff --git a/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt b/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt
index 8a9f7a443..6215ba5e8 100644
--- a/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt
+++ b/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealVectorTest.kt
@@ -1,33 +1,33 @@
 package kaceince.kmath.real
 
-import kscience.kmath.linear.MatrixContext
-import kscience.kmath.linear.asMatrix
-import kscience.kmath.linear.transpose
+import kscience.kmath.linear.*
 import kscience.kmath.operations.invoke
 import kscience.kmath.real.RealVector
+import kscience.kmath.real.plus
+import kscience.kmath.structures.Buffer
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
 internal class RealVectorTest {
     @Test
     fun testSum() {
-        val vector1 = RealVector(5) { it.toDouble() }
-        val vector2 = RealVector(5) { 5 - it.toDouble() }
+        val vector1 = Buffer.real(5) { it.toDouble() }
+        val vector2 =  Buffer.real(5) { 5 - it.toDouble() }
         val sum = vector1 + vector2
         assertEquals(5.0, sum[2])
     }
 
     @Test
     fun testVectorToMatrix() {
-        val vector = RealVector(5) { it.toDouble() }
+        val vector =  Buffer.real(5) { it.toDouble() }
         val matrix = vector.asMatrix()
         assertEquals(4.0, matrix[4, 0])
     }
 
     @Test
     fun testDot() {
-        val vector1 = RealVector(5) { it.toDouble() }
-        val vector2 = RealVector(5) { 5 - it.toDouble() }
+        val vector1 =  Buffer.real(5) { it.toDouble() }
+        val vector2 =  Buffer.real(5) { 5 - it.toDouble() }
         val matrix1 = vector1.asMatrix()
         val matrix2 = vector2.asMatrix().transpose()
         val product = MatrixContext.real { matrix1 dot matrix2 }
diff --git a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt
index c58952ab4..f95264ee1 100644
--- a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt
+++ b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt
@@ -3,7 +3,6 @@ package kscience.kmath.histogram
 import kscience.kmath.linear.Point
 import kscience.kmath.operations.SpaceOperations
 import kscience.kmath.operations.invoke
-import kscience.kmath.real.asVector
 import kscience.kmath.structures.*
 import kotlin.math.floor
 
@@ -123,8 +122,8 @@ public class RealHistogram(
          *```
          */
         public fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): RealHistogram = RealHistogram(
-            ranges.map(ClosedFloatingPointRange<Double>::start).asVector(),
-            ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asVector()
+            ranges.map(ClosedFloatingPointRange<Double>::start).asBuffer(),
+            ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asBuffer()
         )
 
         /**
diff --git a/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt b/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt
index eebb41019..af22afc6b 100644
--- a/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt
+++ b/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt
@@ -4,6 +4,8 @@ import kscience.kmath.histogram.RealHistogram
 import kscience.kmath.histogram.fill
 import kscience.kmath.histogram.put
 import kscience.kmath.real.RealVector
+import kscience.kmath.real.invoke
+import kscience.kmath.structures.Buffer
 import kotlin.random.Random
 import kotlin.test.*
 
diff --git a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt b/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt
index 5fada1302..2f3855892 100644
--- a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt
+++ b/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt
@@ -1,8 +1,8 @@
 package kscience.kmath.histogram
 
 import kscience.kmath.real.RealVector
-import kscience.kmath.real.asVector
 import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.asBuffer
 import java.util.*
 import kotlin.math.floor
 
@@ -16,7 +16,7 @@ public class UnivariateBin(
     //TODO add weighting
     public override val value: Number get() = counter.sum()
 
-    public override val center: RealVector get() = doubleArrayOf(position).asVector()
+    public override val center: RealVector get() = doubleArrayOf(position).asBuffer()
     public override val dimension: Int get() = 1
 
     public operator fun contains(value: Double): Boolean = value in (position - size / 2)..(position + size / 2)
diff --git a/kmath-kotlingrad/build.gradle.kts b/kmath-kotlingrad/build.gradle.kts
index 027a03bc9..3925a744c 100644
--- a/kmath-kotlingrad/build.gradle.kts
+++ b/kmath-kotlingrad/build.gradle.kts
@@ -3,7 +3,7 @@ plugins {
 }
 
 dependencies {
-    implementation("com.github.breandan:kaliningraph:0.1.2")
-    implementation("com.github.breandan:kotlingrad:0.3.7")
+    implementation("com.github.breandan:kaliningraph:0.1.4")
+    implementation("com.github.breandan:kotlingrad:0.4.0")
     api(project(":kmath-ast"))
 }
diff --git a/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/DifferentiableMstExpression.kt b/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/DifferentiableMstExpression.kt
index dd5e46f90..abde9e54d 100644
--- a/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/DifferentiableMstExpression.kt
+++ b/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/DifferentiableMstExpression.kt
@@ -1,6 +1,6 @@
 package kscience.kmath.kotlingrad
 
-import edu.umontreal.kotlingrad.experimental.SFun
+import edu.umontreal.kotlingrad.api.SFun
 import kscience.kmath.ast.MST
 import kscience.kmath.ast.MstAlgebra
 import kscience.kmath.ast.MstExpression
diff --git a/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/KMathNumber.kt b/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/KMathNumber.kt
index ce5658137..2a4db4258 100644
--- a/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/KMathNumber.kt
+++ b/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/KMathNumber.kt
@@ -1,7 +1,7 @@
 package kscience.kmath.kotlingrad
 
-import edu.umontreal.kotlingrad.experimental.RealNumber
-import edu.umontreal.kotlingrad.experimental.SConst
+import edu.umontreal.kotlingrad.api.RealNumber
+import edu.umontreal.kotlingrad.api.SConst
 import kscience.kmath.operations.NumericAlgebra
 
 /**
diff --git a/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/ScalarsAdapters.kt b/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/ScalarsAdapters.kt
index b6effab4b..8dc1d3958 100644
--- a/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/ScalarsAdapters.kt
+++ b/kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/ScalarsAdapters.kt
@@ -1,6 +1,6 @@
 package kscience.kmath.kotlingrad
 
-import edu.umontreal.kotlingrad.experimental.*
+import edu.umontreal.kotlingrad.api.*
 import kscience.kmath.ast.MST
 import kscience.kmath.ast.MstAlgebra
 import kscience.kmath.ast.MstExtendedField
diff --git a/kmath-kotlingrad/src/test/kotlin/kscience/kmath/kotlingrad/AdaptingTests.kt b/kmath-kotlingrad/src/test/kotlin/kscience/kmath/kotlingrad/AdaptingTests.kt
index 77902211b..aa4ddd703 100644
--- a/kmath-kotlingrad/src/test/kotlin/kscience/kmath/kotlingrad/AdaptingTests.kt
+++ b/kmath-kotlingrad/src/test/kotlin/kscience/kmath/kotlingrad/AdaptingTests.kt
@@ -1,6 +1,6 @@
 package kscience.kmath.kotlingrad
 
-import edu.umontreal.kotlingrad.experimental.*
+import edu.umontreal.kotlingrad.api.*
 import kscience.kmath.asm.compile
 import kscience.kmath.ast.MstAlgebra
 import kscience.kmath.ast.MstExpression
diff --git a/kmath-nd4j/README.md b/kmath-nd4j/README.md
index b071df499..176dfc09d 100644
--- a/kmath-nd4j/README.md
+++ b/kmath-nd4j/README.md
@@ -2,14 +2,14 @@
 
 This subproject implements the following features:
 
- - [nd4jarraystrucure](src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt) : NDStructure wrapper for INDArray
+ - [nd4jarraystructure](src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt) : NDStructure wrapper for INDArray
  - [nd4jarrayrings](src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt) : Rings over Nd4jArrayStructure of Int and Long
  - [nd4jarrayfields](src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : Fields over Nd4jArrayStructure of Float and Double
 
 
 > #### Artifact:
 >
-> This module artifact: `kscience.kmath:kmath-nd4j:0.2.0-dev-3`.
+> This module artifact: `kscience.kmath:kmath-nd4j:0.2.0-dev-4`.
 >
 > Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-nd4j/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-nd4j/_latestVersion)
 >
@@ -27,7 +27,7 @@ This subproject implements the following features:
 > }
 > 
 > dependencies {
->     implementation 'kscience.kmath:kmath-nd4j:0.2.0-dev-3'
+>     implementation 'kscience.kmath:kmath-nd4j:0.2.0-dev-4'
 > }
 > ```
 > **Gradle Kotlin DSL:**
@@ -41,7 +41,7 @@ This subproject implements the following features:
 > }
 > 
 > dependencies {
->     implementation("kscience.kmath:kmath-nd4j:0.2.0-dev-3")
+>     implementation("kscience.kmath:kmath-nd4j:0.2.0-dev-4")
 > }
 > ```
 
diff --git a/kmath-stat/src/jvmMain/kotlin/kscience/kmath/stat/distributions.kt b/kmath-stat/src/jvmMain/kotlin/kscience/kmath/stat/distributions.kt
index 9a77b0bd2..6cc18a37c 100644
--- a/kmath-stat/src/jvmMain/kotlin/kscience/kmath/stat/distributions.kt
+++ b/kmath-stat/src/jvmMain/kotlin/kscience/kmath/stat/distributions.kt
@@ -51,7 +51,7 @@ private fun normalSampler(method: NormalSamplerMethod, provider: UniformRandomPr
 
 public fun Distribution.Companion.normal(
     method: NormalSamplerMethod = NormalSamplerMethod.Ziggurat
-): Distribution<Double> = object : ContinuousSamplerDistribution() {
+): ContinuousSamplerDistribution = object : ContinuousSamplerDistribution() {
     override fun buildCMSampler(generator: RandomGenerator): ContinuousSampler {
         val provider = generator.asUniformRandomProvider()
         return normalSampler(method, provider)
@@ -60,6 +60,9 @@ public fun Distribution.Companion.normal(
     override fun probability(arg: Double): Double = exp(-arg.pow(2) / 2) / sqrt(PI * 2)
 }
 
+/**
+ * A univariate normal distribution with given [mean] and [sigma]. [method] defines commons-rng generation method
+ */
 public fun Distribution.Companion.normal(
     mean: Double,
     sigma: Double,