Merge pull request #2 from SciProgCentre/dev

Dev

.github/workflows/build.yml

@@ -7,26 +7,18 @@ on:
 
 jobs:
   build:
-    strategy:
-      matrix:
-        os: [ macOS-latest, windows-latest ]
-    runs-on: ${{matrix.os}}
-    timeout-minutes: 40
+    runs-on: windows-latest
+    timeout-minutes: 20
     steps:
-      - uses: actions/checkout@v3.0.0
-      - uses: actions/setup-java@v3.0.0
+      - uses: actions/checkout@v3
+      - uses: actions/setup-java@v3.5.1
         with:
-          java-version: 11
-          distribution: liberica
-      - name: Cache konan
-        uses: actions/cache@v3.0.1
-        with:
-          path: ~/.konan
-          key: ${{ runner.os }}-gradle-${{ hashFiles('*.gradle.kts') }}
-          restore-keys: |
-            ${{ runner.os }}-gradle-
+          java-version: '11'
+          distribution: 'liberica'
+          cache: 'gradle'
       - name: Gradle Wrapper Validation
         uses: gradle/wrapper-validation-action@v1.0.4
-      - uses: gradle/gradle-build-action@v2.1.5
+      - name: Gradle Build
+        uses: gradle/gradle-build-action@v2.3.2
         with:
-          arguments: build
+          arguments: test jvmTest

CHANGELOG.md

@@ -2,11 +2,15 @@
 
 ## [Unreleased]
 ### Added
+- Type-aliases for numbers like `Float64`
 - 2D optimal trajectory computation in a separate module `kmath-trajectory`
 - Autodiff for generic algebra elements in core!
 - Algebra now has an obligatory `bufferFactory` (#477).
 
 ### Changed
+- Tensor operations switched to prefix notation
+- Row-wise and column-wise ND shapes in the core
+- Shape is read-only
 - Major refactor of tensors (only minor API changes)
 - Kotlin 1.7.20
 - `LazyStructure` `deffered` -> `async` to comply with coroutines code style
@@ -16,6 +20,7 @@
 ### Deprecated
 
 ### Removed
+- Polynomials moved to https://github.com/SciProgCentre/kmath-polynomial
 
 ### Fixed
 

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

@@ -6,34 +6,75 @@
 package space.kscience.kmath.benchmarks
 
 import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
 import kotlinx.benchmark.Scope
 import kotlinx.benchmark.State
 import space.kscience.kmath.complex.Complex
+import space.kscience.kmath.complex.ComplexField
 import space.kscience.kmath.complex.complex
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.MutableBuffer
+import space.kscience.kmath.structures.getDouble
+import space.kscience.kmath.structures.permute
 
 @State(Scope.Benchmark)
 internal class BufferBenchmark {
-    @Benchmark
-    fun genericDoubleBufferReadWrite() {
-        val buffer = DoubleBuffer(size) { it.toDouble() }
 
+    @Benchmark
+    fun doubleArrayReadWrite(blackhole: Blackhole) {
+        val buffer = DoubleArray(size) { it.toDouble() }
+        var res = 0.0
         (0 until size).forEach {
-            buffer[it]
+            res += buffer[it]
         }
+        blackhole.consume(res)
     }
 
     @Benchmark
-    fun complexBufferReadWrite() {
-        val buffer = MutableBuffer.complex(size / 2) { Complex(it.toDouble(), -it.toDouble()) }
-
-        (0 until size / 2).forEach {
-            buffer[it]
+    fun doubleBufferReadWrite(blackhole: Blackhole) {
+        val buffer = DoubleBuffer(size) { it.toDouble() }
+        var res = 0.0
+        (0 until size).forEach {
+            res += buffer[it]
         }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun bufferViewReadWrite(blackhole: Blackhole) {
+        val buffer = DoubleBuffer(size) { it.toDouble() }.permute(reversedIndices)
+        var res = 0.0
+        (0 until size).forEach {
+            res += buffer[it]
+        }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun bufferViewReadWriteSpecialized(blackhole: Blackhole) {
+        val buffer = DoubleBuffer(size) { it.toDouble() }.permute(reversedIndices)
+        var res = 0.0
+        (0 until size).forEach {
+            res += buffer.getDouble(it)
+        }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun complexBufferReadWrite(blackhole: Blackhole) = ComplexField {
+        val buffer = Buffer.complex(size / 2) { Complex(it.toDouble(), -it.toDouble()) }
+
+        var res = zero
+        (0 until size / 2).forEach {
+            res += buffer[it]
+        }
+
+        blackhole.consume(res)
     }
 
     private companion object {
         private const val size = 100
+        private val reversedIndices = IntArray(size){it}.apply { reverse() }
     }
 }

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

@@ -13,10 +13,8 @@ import org.jetbrains.kotlinx.multik.api.Multik
 import org.jetbrains.kotlinx.multik.api.ones
 import org.jetbrains.kotlinx.multik.ndarray.data.DN
 import org.jetbrains.kotlinx.multik.ndarray.data.DataType
-import space.kscience.kmath.nd.BufferedFieldOpsND
-import space.kscience.kmath.nd.StructureND
-import space.kscience.kmath.nd.ndAlgebra
-import space.kscience.kmath.nd.one
+import space.kscience.kmath.misc.UnsafeKMathAPI
+import space.kscience.kmath.nd.*
 import space.kscience.kmath.nd4j.nd4j
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.tensors.core.DoubleTensor
@@ -69,9 +67,10 @@ internal class NDFieldBenchmark {
         blackhole.consume(res)
     }
 
+    @OptIn(UnsafeKMathAPI::class)
     @Benchmark
     fun multikInPlaceAdd(blackhole: Blackhole) = with(multikAlgebra) {
-        val res = Multik.ones<Double, DN>(shape, DataType.DoubleDataType).wrap()
+        val res = Multik.ones<Double, DN>(shape.asArray(), DataType.DoubleDataType).wrap()
         repeat(n) { res += 1.0 }
         blackhole.consume(res)
     }
@@ -86,7 +85,7 @@ internal class NDFieldBenchmark {
     private companion object {
         private const val dim = 1000
         private const val n = 100
-        private val shape = intArrayOf(dim, dim)
+        private val shape = ShapeND(dim, dim)
         private val specializedField = DoubleField.ndAlgebra
         private val genericField = BufferedFieldOpsND(DoubleField)
         private val nd4jField = DoubleField.nd4j

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

@@ -13,6 +13,8 @@ import space.kscience.kmath.linear.linearSpace
 import space.kscience.kmath.linear.matrix
 import space.kscience.kmath.linear.symmetric
 import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.tensors.core.symEigJacobi
+import space.kscience.kmath.tensors.core.symEigSvd
 import space.kscience.kmath.tensors.core.tensorAlgebra
 import kotlin.random.Random
 
@@ -27,11 +29,11 @@ internal class TensorAlgebraBenchmark {
 
     @Benchmark
     fun tensorSymEigSvd(blackhole: Blackhole) = with(Double.tensorAlgebra) {
-        blackhole.consume(matrix.symEigSvd(1e-10))
+        blackhole.consume(symEigSvd(matrix, 1e-10))
     }
 
     @Benchmark
     fun tensorSymEigJacobi(blackhole: Blackhole) = with(Double.tensorAlgebra) {
-        blackhole.consume(matrix.symEigJacobi(50, 1e-10))
+        blackhole.consume(symEigJacobi(matrix, 50, 1e-10))
     }
 }

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

@@ -10,7 +10,7 @@ import kotlinx.benchmark.Blackhole
 import kotlinx.benchmark.Scope
 import kotlinx.benchmark.State
 import org.jetbrains.bio.viktor.F64Array
-import space.kscience.kmath.nd.Shape
+import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.nd.ndAlgebra
 import space.kscience.kmath.nd.one
@@ -49,7 +49,7 @@ internal class ViktorBenchmark {
     private companion object {
         private const val dim = 1000
         private const val n = 100
-        private val shape = Shape(dim, dim)
+        private val shape = ShapeND(dim, dim)
 
         // automatically build context most suited for given type.
         private val doubleField = DoubleField.ndAlgebra

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

@@ -10,7 +10,7 @@ import kotlinx.benchmark.Blackhole
 import kotlinx.benchmark.Scope
 import kotlinx.benchmark.State
 import org.jetbrains.bio.viktor.F64Array
-import space.kscience.kmath.nd.Shape
+import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.nd.ndAlgebra
 import space.kscience.kmath.nd.one
 import space.kscience.kmath.operations.DoubleField
@@ -49,7 +49,7 @@ internal class ViktorLogBenchmark {
     private companion object {
         private const val dim = 1000
         private const val n = 100
-        private val shape = Shape(dim, dim)
+        private val shape = ShapeND(dim, dim)
 
         // automatically build context most suited for given type.
         private val doubleField = DoubleField.ndAlgebra

build.gradle.kts

@@ -15,7 +15,7 @@ allprojects {
     }
 
     group = "space.kscience"
-    version = "0.3.1-dev-4"
+    version = "0.3.1-dev-7"
 }
 
 subprojects {
@@ -79,9 +79,9 @@ ksciencePublish {
     github("kmath", "SciProgCentre")
     space(
         if (isInDevelopment) {
-            "https://maven.pkg.jetbrains.space/mipt-npm/p/sci/dev"
+            "https://maven.pkg.jetbrains.space/spc/p/sci/dev"
         } else {
-            "https://maven.pkg.jetbrains.space/mipt-npm/p/sci/release"
+            "https://maven.pkg.jetbrains.space/spc/p/sci/maven"
         }
     )
     sonatype()

docs/templates/ARTIFACT-TEMPLATE.md

@@ -3,10 +3,12 @@
 The Maven coordinates of this project are `${group}:${name}:${version}`.
 
 **Gradle:**
-```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 {
@@ -18,6 +20,8 @@ dependencies {
 repositories {
     maven("https://repo.kotlin.link")
     mavenCentral()
+    // development and snapshot versions
+    maven("https://maven.pkg.jetbrains.space/spc/p/sci/dev")
 }
 
 dependencies {

examples/build.gradle.kts

@@ -18,7 +18,6 @@ dependencies {
     implementation(project(":kmath-commons"))
     implementation(project(":kmath-complex"))
     implementation(project(":kmath-functions"))
-    implementation(project(":kmath-polynomial"))
     implementation(project(":kmath-optimization"))
     implementation(project(":kmath-stat"))
     implementation(project(":kmath-viktor"))

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

@@ -18,10 +18,10 @@ import space.kscience.kmath.optimization.FunctionOptimizationTarget
 import space.kscience.kmath.optimization.optimizeWith
 import space.kscience.kmath.optimization.resultPoint
 import space.kscience.kmath.optimization.resultValue
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.real.DoubleVector
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
-import space.kscience.kmath.stat.RandomGenerator
 import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
 import space.kscience.plotly.models.TraceValues

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

@@ -19,9 +19,9 @@ import space.kscience.kmath.optimization.QowOptimizer
 import space.kscience.kmath.optimization.chiSquaredOrNull
 import space.kscience.kmath.optimization.fitWith
 import space.kscience.kmath.optimization.resultPoint
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
-import space.kscience.kmath.stat.RandomGenerator
 import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
 import kotlin.math.abs

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

@@ -1,399 +0,0 @@
-/*
- * Copyright 2018-2022 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("LocalVariableName")
-
-package space.kscience.kmath.functions
-
-import space.kscience.kmath.expressions.Symbol
-import space.kscience.kmath.expressions.symbol
-import space.kscience.kmath.operations.algebra
-import space.kscience.kmath.operations.invoke
-
-
-/**
- * Shows [ListPolynomial]s' and [ListRationalFunction]s' capabilities.
- */
-fun listPolynomialsExample() {
-    // [ListPolynomial] is a representation of a univariate polynomial as a list of coefficients from the least term to
-    // the greatest term. For example,
-    val polynomial1: ListPolynomial<Int> = ListPolynomial(listOf(2, -3, 1))
-    // represents polynomial 2 + (-3) x + x^2
-
-    // There are also shortcut fabrics:
-    val polynomial2: ListPolynomial<Int> = ListPolynomial(2, -3, 1)
-    println(polynomial1 == polynomial2) // true
-    // and even
-    val polynomial3: ListPolynomial<Int> = 57.asListPolynomial()
-    val polynomial4: ListPolynomial<Int> = ListPolynomial(listOf(57))
-    println(polynomial3 == polynomial4) // true
-
-    val polynomial5: ListPolynomial<Int> = ListPolynomial(3, -1)
-    // For every ring there can be provided a polynomial ring:
-    Int.algebra.listPolynomialSpace {
-        println(-polynomial5 == ListPolynomial(-3, 1)) // true
-        println(polynomial1 + polynomial5 == ListPolynomial(5, -4, 1)) // true
-        println(polynomial1 - polynomial5 == ListPolynomial(-1, -2, 1)) // true
-        println(polynomial1 * polynomial5 == ListPolynomial(6, -11, 6, -1)) // true
-    }
-    // You can even write
-    val x: ListPolynomial<Double> = ListPolynomial(0.0, 1.0)
-    val polynomial6: ListPolynomial<Double> = ListPolynomial(2.0, -3.0, 1.0)
-    Double.algebra.listPolynomialSpace {
-        println(2 - 3 * x + x * x == polynomial6)
-        println(2.0 - 3.0 * x + x * x == polynomial6)
-    }
-
-    // Also there are some utilities for polynomials:
-    println(polynomial1.substitute(Int.algebra, 1) == 0) // true, because 2 + (-3) * 1 + 1^2 = 0
-    println(polynomial1.substitute(Int.algebra, polynomial5) == polynomial1) // true, because 2 + (-3) * (3-x) + (3-x)^2 = 2 - 3x + x^2
-    println(polynomial1.derivative(Int.algebra) == ListPolynomial(-3, 2)) // true, (2 - 3x + x^2)' = -3 + 2x
-    println(polynomial1.nthDerivative(Int.algebra, 2) == 2.asListPolynomial()) // true, (2 - 3x + x^2)'' = 2
-
-    // Lastly, there are rational functions and some other utilities:
-    Double.algebra.listRationalFunctionSpace {
-        val rationalFunction1: ListRationalFunction<Double> = ListRationalFunction(listOf(2.0, -3.0, 1.0), listOf(3.0, -1.0))
-        // It's just (2 - 3x + x^2)/(3 - x)
-
-        val rationalFunction2 : ListRationalFunction<Double> = ListRationalFunction(listOf(5.0, -4.0, 1.0), listOf(3.0, -1.0))
-        // It's just (5 - 4x + x^2)/(3 - x)
-
-        println(rationalFunction1 + 1 == rationalFunction2)
-    }
-}
-
-/**
- * Shows [NumberedPolynomial]s' and [NumberedRationalFunction]s' capabilities.
- */
-fun numberedPolynomialsExample() {
-    // Consider polynomial
-    // 3 + 5 x_2 - 7 x_1^2 x_3
-    // Consider, for example, its term -7 x_1^2 x_3. -7 is a coefficient of the term, whereas (2, 0, 1, 0, 0, ...) is
-    // description of degrees of variables x_1, x_2, ... in the term. Such description with removed leading zeros
-    // [2, 0, 1] is called "signature" of the term -7 x_1^2 x_3.
-
-    val polynomial1: NumberedPolynomial<Int>
-    with(Int.algebra) {
-        // [NumberedPolynomial] is a representation of a multivariate polynomial, that stores terms in a map with terms'
-        // signatures as the map's keys and terms' coefficients as corresponding values. For example,
-        polynomial1 = NumberedPolynomial(
-            mapOf(
-                listOf<UInt>() to 3,
-                listOf(0u, 1u) to 5,
-                listOf(2u, 0u, 1u) to -7,
-            )
-        )
-        // represents polynomial 3 + 5 x_2 - 7 x_1^2 x_3
-
-        // This `NumberedPolynomial` function needs context of either ring of constant (as `Int.algebra` in this example)
-        // or space of NumberedPolynomials over it. To understand why it is like this see documentations of functions
-        // NumberedPolynomial and NumberedPolynomialWithoutCheck
-
-        // There are also shortcut fabrics:
-        val polynomial2: NumberedPolynomial<Int> = NumberedPolynomial(
-            listOf<UInt>() to 3,
-            listOf(0u, 1u) to 5,
-            listOf(2u, 0u, 1u) to -7,
-        )
-        println(polynomial1 == polynomial2) // true
-        // and even
-        val polynomial3: NumberedPolynomial<Int> = 57.asNumberedPolynomial() // This one actually does not algebraic context!
-        val polynomial4: NumberedPolynomial<Int> = NumberedPolynomial(listOf<UInt>() to 57)
-        println(polynomial3 == polynomial4) // true
-
-        numberedPolynomialSpace {
-            // Also there is DSL for constructing NumberedPolynomials:
-            val polynomial5: NumberedPolynomial<Int> = NumberedPolynomialDSL1 {
-                3 {}
-                5 { 1 inPowerOf 1u }
-                -7 with { 0 pow 2u; 2 pow 1u }
-                // `pow` and `inPowerOf` are the same
-                // `with` is omittable
-            }
-            println(polynomial1 == polynomial5) // true
-
-            // Unfortunately the DSL does not work good in bare context of constants' ring, so for now it's disabled and
-            // works only in NumberedPolynomialSpace and NumberedRationalFunctionSpace
-        }
-    }
-
-    val polynomial6: NumberedPolynomial<Int> = Int.algebra {
-        NumberedPolynomial(
-            listOf<UInt>() to 7,
-            listOf(0u, 1u) to -5,
-            listOf(2u, 0u, 1u) to 0,
-            listOf(0u, 0u, 0u, 4u) to 4,
-        )
-    }
-    // For every ring there can be provided a polynomial ring:
-    Int.algebra.numberedPolynomialSpace {
-        println(
-            -polynomial6 == NumberedPolynomial(
-                listOf<UInt>() to -7,
-                listOf(0u, 1u) to 5,
-                listOf(2u, 0u, 1u) to 0,
-                listOf(0u, 0u, 0u, 4u) to (-4),
-            )
-        ) // true
-        println(
-            polynomial1 + polynomial6 == NumberedPolynomial(
-                listOf<UInt>() to 10,
-                listOf(0u, 1u) to 0,
-                listOf(2u, 0u, 1u) to -7,
-                listOf(0u, 0u, 0u, 4u) to 4,
-            )
-        ) // true
-        println(
-            polynomial1 - polynomial6 == NumberedPolynomial(
-                listOf<UInt>() to -4,
-                listOf(0u, 1u) to 10,
-                listOf(2u, 0u, 1u) to -7,
-                listOf(0u, 0u, 0u, 4u) to -4,
-            )
-        ) // true
-
-        polynomial1 * polynomial6 // Multiplication works too
-    }
-
-    Double.algebra.numberedPolynomialSpace {
-        // You can even write
-        val x_1: NumberedPolynomial<Double> = NumberedPolynomial(listOf(1u) to 1.0)
-        val x_2: NumberedPolynomial<Double> = NumberedPolynomial(listOf(0u, 1u) to 1.0)
-        val x_3: NumberedPolynomial<Double> = NumberedPolynomial(listOf(0u, 0u, 1u) to 1.0)
-        val polynomial7: NumberedPolynomial<Double> = NumberedPolynomial(
-            listOf<UInt>() to 3.0,
-            listOf(0u, 1u) to 5.0,
-            listOf(2u, 0u, 1u) to -7.0,
-        )
-        Double.algebra.listPolynomialSpace {
-            println(3 + 5 * x_2 - 7 * x_1 * x_1 * x_3 == polynomial7)
-            println(3.0 + 5.0 * x_2 - 7.0 * x_1 * x_1 * x_3 == polynomial7)
-        }
-    }
-
-    Int.algebra.numberedPolynomialSpace {
-        val x_4: NumberedPolynomial<Int> = NumberedPolynomial(listOf(0u, 0u, 0u, 4u) to 1)
-        // Also there are some utilities for polynomials:
-        println(polynomial1.substitute(mapOf(0 to 1, 1 to -2, 2 to -1)) == 0.asNumberedPolynomial()) // true,
-            // because it's substitution x_1 -> 1, x_2 -> -2, x_3 -> -1,
-            // so 3 + 5 x_2 - 7 x_1^2 x_3 = 3 + 5 * (-2) - 7 * 1^2 * (-1) = 3 - 10 + 7 = 0
-        println(
-            polynomial1.substitute(mapOf(1 to x_4)) == NumberedPolynomial(
-                listOf<UInt>() to 3,
-                listOf(0u, 1u) to 5,
-                listOf(2u, 0u, 1u) to -7,
-            )
-        ) // true, because it's substitution x_2 -> x_4, so result is 3 + 5 x_4 - 7 x_1^2 x_3
-        println(
-            polynomial1.derivativeWithRespectTo(Int.algebra, 1) ==
-                    NumberedPolynomial(listOf<UInt>() to 5)
-        ) // true, d/dx_2 (3 + 5 x_2 - 7 x_1^2 x_3) = 5
-    }
-
-    // Lastly, there are rational functions and some other utilities:
-    Double.algebra.numberedRationalFunctionSpace {
-        val rationalFunction1: NumberedRationalFunction<Double> = NumberedRationalFunction(
-            NumberedPolynomial(
-                listOf<UInt>() to 2.0,
-                listOf(1u) to -3.0,
-                listOf(2u) to 1.0,
-            ),
-            NumberedPolynomial(
-                listOf<UInt>() to 3.0,
-                listOf(1u) to -1.0,
-            )
-        )
-        // It's just (2 - 3x + x^2)/(3 - x) where x = x_1
-
-        val rationalFunction2: NumberedRationalFunction<Double> = NumberedRationalFunction(
-            NumberedPolynomial(
-                listOf<UInt>() to 5.0,
-                listOf(1u) to -4.0,
-                listOf(2u) to 1.0,
-            ),
-            NumberedPolynomial(
-                listOf<UInt>() to 3.0,
-                listOf(1u) to -1.0,
-            )
-        )
-        // It's just (5 - 4x + x^2)/(3 - x) where x = x_1
-
-        println(rationalFunction1 + 1 == rationalFunction2)
-    }
-}
-
-/**
- * Shows [LabeledPolynomial]s' and [LabeledRationalFunction]s' capabilities.
- */
-fun labeledPolynomialsExample() {
-    val x by symbol
-    val y by symbol
-    val z by symbol
-    val t by symbol
-
-    // Consider polynomial
-    // 3 + 5 y - 7 x^2 z
-    // Consider, for example, its term -7 x^2 z. -7 is a coefficient of the term, whereas matching (x -> 2, z -> 3) is
-    // description of degrees of variables x_1, x_2, ... in the term. Such description is called "signature" of the
-    // term -7 x_1^2 x_3.
-
-    val polynomial1: LabeledPolynomial<Int>
-    with(Int.algebra) {
-        // [LabeledPolynomial] is a representation of a multivariate polynomial, that stores terms in a map with terms'
-        // signatures as the map's keys and terms' coefficients as corresponding values. For example,
-        polynomial1 = LabeledPolynomial(
-            mapOf(
-                mapOf<Symbol, UInt>() to 3,
-                mapOf(y to 1u) to 5,
-                mapOf(x to 2u, z to 1u) to -7,
-            )
-        )
-        // represents polynomial 3 + 5 y - 7 x^2 z
-
-        // This `LabeledPolynomial` function needs context of either ring of constant (as `Int.algebra` in this example)
-        // or space of LabeledPolynomials over it. To understand why it is like this see documentations of functions
-        // LabeledPolynomial and LabeledPolynomialWithoutCheck
-
-        // There are also shortcut fabrics:
-        val polynomial2: LabeledPolynomial<Int> = LabeledPolynomial(
-            mapOf<Symbol, UInt>() to 3,
-            mapOf(y to 1u) to 5,
-            mapOf(x to 2u, z to 1u) to -7,
-        )
-        println(polynomial1 == polynomial2) // true
-        // and even
-        val polynomial3: LabeledPolynomial<Int> = 57.asLabeledPolynomial() // This one actually does not algebraic context!
-        val polynomial4: LabeledPolynomial<Int> = LabeledPolynomial(mapOf<Symbol, UInt>() to 57)
-        println(polynomial3 == polynomial4) // true
-
-        labeledPolynomialSpace {
-            // Also there is DSL for constructing NumberedPolynomials:
-            val polynomial5: LabeledPolynomial<Int> = LabeledPolynomialDSL1 {
-                3 {}
-                5 { y inPowerOf 1u }
-                -7 with { x pow 2u; z pow 1u }
-                // `pow` and `inPowerOf` are the same
-                // `with` is omittable
-            }
-            println(polynomial1 == polynomial5) // true
-
-            // Unfortunately the DSL does not work good in bare context of constants' ring, so for now it's disabled and
-            // works only in NumberedPolynomialSpace and NumberedRationalFunctionSpace
-        }
-    }
-
-    val polynomial6: LabeledPolynomial<Int> = Int.algebra {
-        LabeledPolynomial(
-            mapOf<Symbol, UInt>() to 7,
-            mapOf(y to 1u) to -5,
-            mapOf(x to 2u, z to 1u) to 0,
-            mapOf(t to 4u) to 4,
-        )
-    }
-    // For every ring there can be provided a polynomial ring:
-    Int.algebra.labeledPolynomialSpace {
-        println(
-            -polynomial6 == LabeledPolynomial(
-                mapOf<Symbol, UInt>() to -7,
-                mapOf(y to 1u) to 5,
-                mapOf(x to 2u, z to 1u) to 0,
-                mapOf(t to 4u) to -4,
-            )
-        ) // true
-        println(
-            polynomial1 + polynomial6 == LabeledPolynomial(
-                mapOf<Symbol, UInt>() to 10,
-                mapOf(y to 1u) to 0,
-                mapOf(x to 2u, z to 1u) to -7,
-                mapOf(t to 4u) to 4,
-            )
-        ) // true
-        println(
-            polynomial1 - polynomial6 == LabeledPolynomial(
-                mapOf<Symbol, UInt>() to -4,
-                mapOf(y to 1u) to 10,
-                mapOf(x to 2u, z to 1u) to -7,
-                mapOf(t to 4u) to -4,
-            )
-        ) // true
-
-        polynomial1 * polynomial6 // Multiplication works too
-    }
-
-    Double.algebra.labeledPolynomialSpace {
-        // You can even write
-        val polynomial7: LabeledPolynomial<Double> = LabeledPolynomial(
-            mapOf<Symbol, UInt>() to 3.0,
-            mapOf(y to 1u) to 5.0,
-            mapOf(x to 2u, z to 1u) to -7.0,
-        )
-        Double.algebra.listPolynomialSpace {
-            println(3 + 5 * y - 7 * x * x * z == polynomial7)
-            println(3.0 + 5.0 * y - 7.0 * x * x * z == polynomial7)
-        }
-    }
-
-    Int.algebra.labeledPolynomialSpace {
-        // Also there are some utilities for polynomials:
-        println(polynomial1.substitute(mapOf(x to 1, y to -2, z to -1)) == 0.asLabeledPolynomial()) // true,
-        // because it's substitution x -> 1, y -> -2, z -> -1,
-        // so 3 + 5 y - 7 x^2 z = 3 + 5 * (-2) - 7 * 1^2 * (-1) = 3 - 10 + 7 = 0
-        println(
-            polynomial1.substitute(mapOf(y to t.asPolynomial())) == LabeledPolynomial(
-                mapOf<Symbol, UInt>() to 3,
-                mapOf(t to 1u) to 5,
-                mapOf(x to 2u, z to 1u) to -7,
-                )
-        ) // true, because it's substitution y -> t, so result is 3 + 5 t - 7 x^2 z
-        println(
-            polynomial1.derivativeWithRespectTo(Int.algebra, y) == LabeledPolynomial(mapOf<Symbol, UInt>() to 5)
-        ) // true, d/dy (3 + 5 y - 7 x^2 z) = 5
-    }
-
-    // Lastly, there are rational functions and some other utilities:
-    Double.algebra.labeledRationalFunctionSpace {
-        val rationalFunction1: LabeledRationalFunction<Double> = LabeledRationalFunction(
-            LabeledPolynomial(
-                mapOf<Symbol, UInt>() to 2.0,
-                mapOf(x to 1u) to -3.0,
-                mapOf(x to 2u) to 1.0,
-            ),
-            LabeledPolynomial(
-                mapOf<Symbol, UInt>() to 3.0,
-                mapOf(x to 1u) to -1.0,
-            )
-        )
-        // It's just (2 - 3x + x^2)/(3 - x)
-
-        val rationalFunction2: LabeledRationalFunction<Double> = LabeledRationalFunction(
-            LabeledPolynomial(
-                mapOf<Symbol, UInt>() to 5.0,
-                mapOf(x to 1u) to -4.0,
-                mapOf(x to 2u) to 1.0,
-            ),
-            LabeledPolynomial(
-                mapOf<Symbol, UInt>() to 3.0,
-                mapOf(x to 1u) to -1.0,
-            )
-        )
-        // It's just (5 - 4x + x^2)/(3 - x)
-
-        println(rationalFunction1 + 1 == rationalFunction2)
-    }
-}
-
-fun main() {
-    println("ListPolynomials:")
-    listPolynomialsExample()
-    println()
-
-    println("NumberedPolynomials:")
-    numberedPolynomialsExample()
-    println()
-
-    println("ListPolynomials:")
-    labeledPolynomialsExample()
-    println()
-}

examples/src/main/kotlin/space/kscience/kmath/operations/mixedNDOperations.kt

@@ -8,14 +8,14 @@ package space.kscience.kmath.operations
 import space.kscience.kmath.commons.linear.CMLinearSpace
 import space.kscience.kmath.linear.matrix
 import space.kscience.kmath.nd.DoubleBufferND
-import space.kscience.kmath.nd.Shape
+import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.nd.Structure2D
 import space.kscience.kmath.nd.ndAlgebra
 import space.kscience.kmath.viktor.ViktorStructureND
 import space.kscience.kmath.viktor.viktorAlgebra
 
 fun main() {
-    val viktorStructure: ViktorStructureND = DoubleField.viktorAlgebra.structureND(Shape(2, 2)) { (i, j) ->
+    val viktorStructure: ViktorStructureND = DoubleField.viktorAlgebra.structureND(ShapeND(2, 2)) { (i, j) ->
         if (i == j) 2.0 else 0.0
     }
 

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

@@ -0,0 +1,49 @@
+package space.kscience.kmath.series
+
+
+import kotlinx.html.FlowContent
+import kotlinx.html.h1
+import space.kscience.kmath.operations.DoubleBufferOps
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
+import space.kscience.kmath.operations.toList
+import space.kscience.kmath.stat.KMComparisonResult
+import space.kscience.kmath.stat.ksComparisonStatistic
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.slice
+import space.kscience.plotly.*
+import kotlin.math.PI
+
+fun main() = with(Double.algebra.bufferAlgebra.seriesAlgebra()) {
+    fun FlowContent.plotSeries(buffer: Buffer<Double>) {
+        val ls = buffer.labels
+        plot {
+            scatter {
+                x.numbers = ls
+                y.numbers = buffer.toList()
+            }
+            layout {
+                xaxis {
+                    range(0.0..100.0)
+                }
+            }
+        }
+    }
+
+
+    val s1 = series(100) { sin(2 * PI * it / 100) + 1.0 }
+    val s2 = s1.slice(20..50).moveTo(40)
+
+    val s3: Buffer<Double> = s1.zip(s2) { l, r -> l + r } //s1 + s2
+    val s4 = DoubleBufferOps.ln(s3)
+
+    val kmTest: KMComparisonResult<Double> = ksComparisonStatistic(s1, s2)
+
+    Plotly.page {
+        h1 { +"This is my plot" }
+        plotSeries(s1)
+        plotSeries(s2)
+        plotSeries(s4)
+    }.makeFile()
+
+}

examples/src/main/kotlin/space/kscience/kmath/stat/DistributionBenchmark.kt

@@ -10,6 +10,7 @@ import kotlinx.coroutines.async
 import kotlinx.coroutines.runBlocking
 import org.apache.commons.rng.sampling.distribution.BoxMullerNormalizedGaussianSampler
 import org.apache.commons.rng.simple.RandomSource
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.samplers.GaussianSampler
 import java.time.Duration
 import java.time.Instant

examples/src/main/kotlin/space/kscience/kmath/stat/DistributionDemo.kt

@@ -9,6 +9,7 @@ import kotlinx.coroutines.runBlocking
 import space.kscience.kmath.chains.Chain
 import space.kscience.kmath.chains.combineWithState
 import space.kscience.kmath.distributions.NormalDistribution
+import space.kscience.kmath.random.RandomGenerator
 
 private data class AveragingChainState(var num: Int = 0, var value: Double = 0.0)
 

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

@@ -29,7 +29,7 @@ fun main() {
     Nd4j.zeros(0)
     val dim = 1000
     val n = 1000
-    val shape = Shape(dim, dim)
+    val shape = ShapeND(dim, dim)
 
 
     // specialized nd-field for Double. It works as generic Double field as well.

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

@@ -17,11 +17,11 @@ import java.util.stream.IntStream
  * A demonstration implementation of NDField over Real using Java [java.util.stream.DoubleStream] for parallel
  * execution.
  */
-class StreamDoubleFieldND(override val shape: IntArray) : FieldND<Double, DoubleField>,
+class StreamDoubleFieldND(override val shape: ShapeND) : FieldND<Double, DoubleField>,
     NumbersAddOps<StructureND<Double>>,
     ExtendedField<StructureND<Double>> {
 
-    private val strides = DefaultStrides(shape)
+    private val strides = ColumnStrides(shape)
     override val elementAlgebra: DoubleField get() = DoubleField
     override val zero: BufferND<Double> by lazy { structureND(shape) { zero } }
     override val one: BufferND<Double> by lazy { structureND(shape) { one } }
@@ -31,17 +31,19 @@ class StreamDoubleFieldND(override val shape: IntArray) : FieldND<Double, Double
         return structureND(shape) { d }
     }
 
+    @OptIn(PerformancePitfall::class)
     private val StructureND<Double>.buffer: DoubleBuffer
         get() = when {
             !shape.contentEquals(this@StreamDoubleFieldND.shape) -> throw ShapeMismatchException(
                 this@StreamDoubleFieldND.shape,
                 shape
             )
-            this is BufferND && this.indices == this@StreamDoubleFieldND.strides -> this.buffer as DoubleBuffer
+
+            this is BufferND && indices == this@StreamDoubleFieldND.strides -> this.buffer as DoubleBuffer
             else -> DoubleBuffer(strides.linearSize) { offset -> get(strides.index(offset)) }
         }
 
-    override fun structureND(shape: Shape, initializer: DoubleField.(IntArray) -> Double): BufferND<Double> {
+    override fun structureND(shape: ShapeND, initializer: DoubleField.(IntArray) -> Double): BufferND<Double> {
         val array = IntStream.range(0, strides.linearSize).parallel().mapToDouble { offset ->
             val index = strides.index(offset)
             DoubleField.initializer(index)
@@ -109,4 +111,4 @@ class StreamDoubleFieldND(override val shape: IntArray) : FieldND<Double, Double
     override fun atanh(arg: StructureND<Double>): BufferND<Double> = arg.map { atanh(it) }
 }
 
-fun DoubleField.ndStreaming(vararg shape: Int): StreamDoubleFieldND = StreamDoubleFieldND(shape)
+fun DoubleField.ndStreaming(vararg shape: Int): StreamDoubleFieldND = StreamDoubleFieldND(ShapeND(shape))

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

@@ -5,16 +5,19 @@
 
 package space.kscience.kmath.structures
 
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.nd.BufferND
-import space.kscience.kmath.nd.DefaultStrides
+import space.kscience.kmath.nd.ColumnStrides
+import space.kscience.kmath.nd.ShapeND
 import kotlin.system.measureTimeMillis
 
 @Suppress("ASSIGNED_BUT_NEVER_ACCESSED_VARIABLE")
+@OptIn(PerformancePitfall::class)
 fun main() {
     val n = 6000
     val array = DoubleArray(n * n) { 1.0 }
     val buffer = DoubleBuffer(array)
-    val strides = DefaultStrides(intArrayOf(n, n))
+    val strides = ColumnStrides(ShapeND(n, n))
     val structure = BufferND(strides, buffer)
 
     measureTimeMillis {

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

@@ -5,16 +5,23 @@
 
 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.nd.mapToBuffer
+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(),
+    crossinline block: (T) -> R,
+): BufferND<R> = BufferND(indices, buffer.mapToBuffer(bufferFactory, block))
+
 @Suppress("UNUSED_VARIABLE")
 fun main() {
     val n = 6000
-    val structure = StructureND.buffered(intArrayOf(n, n), Buffer.Companion::auto) { 1.0 }
-    structure.mapToBuffer { it + 1 } // warm-up
-    val time1 = measureTimeMillis { val res = structure.mapToBuffer { it + 1 } }
+    val structure = StructureND.buffered(ShapeND(n, n), Buffer.Companion::auto) { 1.0 }
+    structure.mapToBufferND { it + 1 } // warm-up
+    val time1 = measureTimeMillis { val res = structure.mapToBufferND { it + 1 } }
     println("Structure mapping finished in $time1 millis")
     val array = DoubleArray(n * n) { 1.0 }
 

examples/src/main/kotlin/space/kscience/kmath/tensors/OLSWithSVD.kt

@@ -5,16 +5,17 @@
 
 package space.kscience.kmath.tensors
 
-import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.contentEquals
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.tensors.core.DoubleTensor
 import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
-
+import space.kscience.kmath.tensors.core.randomNormal
+import space.kscience.kmath.tensors.core.randomNormalLike
 import kotlin.math.abs
 
 // OLS estimator using SVD
 
-@OptIn(PerformancePitfall::class)
 fun main() {
     //seed for random
     val randSeed = 100500L
@@ -23,10 +24,10 @@ fun main() {
     DoubleTensorAlgebra {
         // take coefficient vector from normal distribution
         val alpha = randomNormal(
-            intArrayOf(5),
+            ShapeND(5),
             randSeed
         ) + fromArray(
-            intArrayOf(5),
+            ShapeND(5),
             doubleArrayOf(1.0, 2.5, 3.4, 5.0, 10.1)
         )
 
@@ -34,27 +35,29 @@ fun main() {
 
         // also take sample of size 20 from normal distribution for x
         val x = randomNormal(
-            intArrayOf(20, 5),
+            ShapeND(20, 5),
             randSeed
         )
 
         // calculate y and add gaussian noise (N(0, 0.05))
         val y = x dot alpha
-        y += y.randomNormalLike(randSeed) * 0.05
+        y += randomNormalLike(y, randSeed) * 0.05
 
         // now restore the coefficient vector with OSL estimator with SVD
-        val (u, singValues, v) = x.svd()
+        val (u, singValues, v) = svd(x)
 
         // we have to make sure the singular values of the matrix are not close to zero
         println("Singular values:\n$singValues")
 
 
         // inverse Sigma matrix can be restored from singular values with diagonalEmbedding function
-        val sigma = diagonalEmbedding(singValues.map{ if (abs(it) < 1e-3) 0.0 else 1.0/it })
+        val sigma = diagonalEmbedding(singValues.map { if (abs(it) < 1e-3) 0.0 else 1.0 / it })
 
         val alphaOLS = v dot sigma dot u.transposed() dot y
-        println("Estimated alpha:\n" +
-                "$alphaOLS")
+        println(
+            "Estimated alpha:\n" +
+                    "$alphaOLS"
+        )
 
         // figure out MSE of approximation
         fun mse(yTrue: DoubleTensor, yPred: DoubleTensor): Double {
@@ -62,7 +65,7 @@ fun main() {
             require(yTrue.shape contentEquals yPred.shape)
 
             val diff = yTrue - yPred
-            return diff.dot(diff).sqrt().value()
+            return sqrt(diff.dot(diff)).value()
         }
 
         println("MSE: ${mse(alpha, alphaOLS)}")

examples/src/main/kotlin/space/kscience/kmath/tensors/PCA.kt

@@ -5,8 +5,8 @@
 
 package space.kscience.kmath.tensors
 
-import space.kscience.kmath.tensors.core.tensorAlgebra
-import space.kscience.kmath.tensors.core.withBroadcast
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.tensors.core.*
 
 
 // simple PCA
@@ -16,12 +16,12 @@ fun main(): Unit = Double.tensorAlgebra.withBroadcast {  // work in context with
 
     // assume x is range from 0 until 10
     val x = fromArray(
-        intArrayOf(10),
+        ShapeND(10),
         DoubleArray(10) { it.toDouble() }
     )
 
     // take y dependent on x with noise
-    val y = 2.0 * x + (3.0 + x.randomNormalLike(seed) * 1.5)
+    val y = 2.0 * x + (3.0 + randomNormalLike(x, seed) * 1.5)
 
     println("x:\n$x")
     println("y:\n$y")
@@ -30,34 +30,34 @@ fun main(): Unit = Double.tensorAlgebra.withBroadcast {  // work in context with
     val dataset = stack(listOf(x, y)).transposed()
 
     // normalize both x and y
-    val xMean = x.mean()
-    val yMean = y.mean()
+    val xMean = mean(x)
+    val yMean = mean(y)
 
-    val xStd = x.std()
-    val yStd = y.std()
+    val xStd = std(x)
+    val yStd = std(y)
 
-    val xScaled = (x - xMean) / xStd
-    val yScaled = (y - yMean) / yStd
+    val xScaled: DoubleTensor = (x - xMean) / xStd
+    val yScaled: DoubleTensor = (y - yMean) / yStd
 
     // save means ans standard deviations for further recovery
     val mean = fromArray(
-        intArrayOf(2),
+        ShapeND(2),
         doubleArrayOf(xMean, yMean)
     )
     println("Means:\n$mean")
 
     val std = fromArray(
-        intArrayOf(2),
+        ShapeND(2),
         doubleArrayOf(xStd, yStd)
     )
     println("Standard deviations:\n$std")
 
     // calculate the covariance matrix of scaled x and y
-    val covMatrix = cov(listOf(xScaled, yScaled))
+    val covMatrix = covariance(listOf(xScaled.asDoubleTensor1D(), yScaled.asDoubleTensor1D()))
     println("Covariance matrix:\n$covMatrix")
 
     // and find out eigenvector of it
-    val (_, evecs) = covMatrix.symEig()
+    val (_, evecs) = symEig(covMatrix)
     val v = evecs.getTensor(0)
     println("Eigenvector:\n$v")
 
@@ -68,7 +68,7 @@ fun main(): Unit = Double.tensorAlgebra.withBroadcast {  // work in context with
     // we can restore original data from reduced data;
     // for example, find 7th element of dataset.
     val n = 7
-    val restored = (datasetReduced.getTensor(n) dot v.view(intArrayOf(1, 2))) * std + mean
+    val restored = (datasetReduced.getTensor(n) dot v.view(ShapeND(1, 2))) * std + mean
     println("Original value:\n${dataset.getTensor(n)}")
     println("Restored value:\n$restored")
 }

examples/src/main/kotlin/space/kscience/kmath/tensors/dataSetNormalization.kt

@@ -5,6 +5,8 @@
 
 package space.kscience.kmath.tensors
 
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.tensors.core.randomNormal
 import space.kscience.kmath.tensors.core.tensorAlgebra
 import space.kscience.kmath.tensors.core.withBroadcast
 
@@ -13,17 +15,17 @@ import space.kscience.kmath.tensors.core.withBroadcast
 
 fun main() = Double.tensorAlgebra.withBroadcast {  // work in context with broadcast methods
     // take dataset of 5-element vectors from normal distribution
-    val dataset = randomNormal(intArrayOf(100, 5)) * 1.5 // all elements from N(0, 1.5)
+    val dataset = randomNormal(ShapeND(100, 5)) * 1.5 // all elements from N(0, 1.5)
 
     dataset += fromArray(
-        intArrayOf(5),
+        ShapeND(5),
         doubleArrayOf(0.0, 1.0, 1.5, 3.0, 5.0) // row means
     )
 
 
     // find out mean and standard deviation of each column
-    val mean = dataset.mean(0, false)
-    val std = dataset.std(0, false)
+    val mean = mean(dataset, 0, false)
+    val std = std(dataset, 0, false)
 
     println("Mean:\n$mean")
     println("Standard deviation:\n$std")
@@ -35,8 +37,8 @@ fun main() = Double.tensorAlgebra.withBroadcast {  // work in context with broad
     // now we can scale dataset with mean normalization
     val datasetScaled = (dataset - mean) / std
 
-    // find out mean and std of scaled dataset
+    // find out mean and standardDiviation of scaled dataset
 
-    println("Mean of scaled:\n${datasetScaled.mean(0, false)}")
-    println("Mean of scaled:\n${datasetScaled.std(0, false)}")
+    println("Mean of scaled:\n${mean(datasetScaled, 0, false)}")
+    println("Mean of scaled:\n${std(datasetScaled, 0, false)}")
 }

examples/src/main/kotlin/space/kscience/kmath/tensors/linearSystemSolvingWithLUP.kt

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.tensors
 
+import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.tensors.core.DoubleTensor
 import space.kscience.kmath.tensors.core.tensorAlgebra
 import space.kscience.kmath.tensors.core.withBroadcast
@@ -15,13 +16,13 @@ fun main() = Double.tensorAlgebra.withBroadcast {// work in context with linear
 
     // set true value of x
     val trueX = fromArray(
-        intArrayOf(4),
+        ShapeND(4),
         doubleArrayOf(-2.0, 1.5, 6.8, -2.4)
     )
 
     // and A matrix
     val a = fromArray(
-        intArrayOf(4, 4),
+        ShapeND(4, 4),
         doubleArrayOf(
             0.5, 10.5, 4.5, 1.0,
             8.5, 0.9, 12.8, 0.1,
@@ -40,7 +41,7 @@ fun main() = Double.tensorAlgebra.withBroadcast {// work in context with linear
     // solve `Ax = b` system using LUP decomposition
 
     // get P, L, U such that PA = LU
-    val (p, l, u) = a.lu()
+    val (p, l, u) = lu(a)
 
     // check P is permutation matrix
     println("P:\n$p")
@@ -64,7 +65,7 @@ fun main() = Double.tensorAlgebra.withBroadcast {// work in context with linear
     // this function returns solution x of a system lx = b, l should be lower triangular
     fun solveLT(l: DoubleTensor, b: DoubleTensor): DoubleTensor {
         val n = l.shape[0]
-        val x = zeros(intArrayOf(n))
+        val x = zeros(ShapeND(n))
         for (i in 0 until n) {
             x[intArrayOf(i)] = (b[intArrayOf(i)] - l.getTensor(i).dot(x).value()) / l[intArrayOf(i, i)]
         }

examples/src/main/kotlin/space/kscience/kmath/tensors/neuralNetwork.kt

@@ -5,12 +5,11 @@
 
 package space.kscience.kmath.tensors
 
+import space.kscience.kmath.nd.ShapeND
+import space.kscience.kmath.nd.contentEquals
 import space.kscience.kmath.operations.asIterable
 import space.kscience.kmath.operations.invoke
-import space.kscience.kmath.tensors.core.BroadcastDoubleTensorAlgebra
-import space.kscience.kmath.tensors.core.DoubleTensor
-import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
-import space.kscience.kmath.tensors.core.toDoubleTensor
+import space.kscience.kmath.tensors.core.*
 import kotlin.math.sqrt
 
 const val seed = 100500L
@@ -49,7 +48,7 @@ fun reluDer(x: DoubleTensor): DoubleTensor = DoubleTensorAlgebra {
 class ReLU : Activation(::relu, ::reluDer)
 
 fun sigmoid(x: DoubleTensor): DoubleTensor = DoubleTensorAlgebra {
-    1.0 / (1.0 + (-x).exp())
+    1.0 / (1.0 + exp((-x)))
 }
 
 fun sigmoidDer(x: DoubleTensor): DoubleTensor = DoubleTensorAlgebra {
@@ -68,12 +67,12 @@ class Dense(
 
     private val weights: DoubleTensor = DoubleTensorAlgebra {
         randomNormal(
-            intArrayOf(inputUnits, outputUnits),
+            ShapeND(inputUnits, outputUnits),
             seed
         ) * sqrt(2.0 / (inputUnits + outputUnits))
     }
 
-    private val bias: DoubleTensor = DoubleTensorAlgebra { zeros(intArrayOf(outputUnits)) }
+    private val bias: DoubleTensor = DoubleTensorAlgebra { zeros(ShapeND(outputUnits)) }
 
     override fun forward(input: DoubleTensor): DoubleTensor = BroadcastDoubleTensorAlgebra {
         (input dot weights) + bias
@@ -83,7 +82,7 @@ class Dense(
         val gradInput = outputError dot weights.transposed()
 
         val gradW = input.transposed() dot outputError
-        val gradBias = outputError.mean(dim = 0, keepDim = false) * input.shape[0].toDouble()
+        val gradBias = mean(structureND = outputError, dim = 0, keepDim = false) * input.shape[0].toDouble()
 
         weights -= learningRate * gradW
         bias -= learningRate * gradBias
@@ -116,7 +115,7 @@ class NeuralNetwork(private val layers: List<Layer>) {
             onesForAnswers[intArrayOf(index, label)] = 1.0
         }
 
-        val softmaxValue = yPred.exp() / yPred.exp().sum(dim = 1, keepDim = true)
+        val softmaxValue = exp(yPred) / exp(yPred).sum(dim = 1, keepDim = true)
 
         (-onesForAnswers + softmaxValue) / (yPred.shape[0].toDouble())
     }
@@ -174,7 +173,6 @@ class NeuralNetwork(private val layers: List<Layer>) {
 }
 
 
-@OptIn(ExperimentalStdlibApi::class)
 fun main() = BroadcastDoubleTensorAlgebra {
     val features = 5
     val sampleSize = 250
@@ -182,17 +180,17 @@ fun main() = BroadcastDoubleTensorAlgebra {
     //val testSize = sampleSize - trainSize
 
     // take sample of features from normal distribution
-    val x = randomNormal(intArrayOf(sampleSize, features), seed) * 2.5
+    val x = randomNormal(ShapeND(sampleSize, features), seed) * 2.5
 
     x += fromArray(
-        intArrayOf(5),
+        ShapeND(5),
         doubleArrayOf(0.0, -1.0, -2.5, -3.0, 5.5) // row means
     )
 
 
     // define class like '1' if the sum of features > 0 and '0' otherwise
     val y = fromArray(
-        intArrayOf(sampleSize, 1),
+        ShapeND(sampleSize, 1),
         DoubleArray(sampleSize) { i ->
             if (x.getTensor(i).sum() > 0.0) {
                 1.0

gradle.properties

@@ -3,13 +3,15 @@
 # Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
 #
 kotlin.code.style=official
-kotlin.jupyter.add.scanner=false
 kotlin.mpp.stability.nowarn=true
 kotlin.native.ignoreDisabledTargets=true
 kotlin.incremental.js.ir=true
 
 org.gradle.configureondemand=true
-org.gradle.parallel=true
 org.gradle.jvmargs=-Xmx4096m
 
-toolsVersion=0.13.0-kotlin-1.7.20-Beta
+toolsVersion=0.13.1-kotlin-1.7.20
+
+
+org.gradle.parallel=true
+org.gradle.workers.max=4

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/random/CMRandomGeneratorWrapper.kt

@@ -6,12 +6,12 @@
 package space.kscience.kmath.commons.random
 
 import kotlinx.coroutines.runBlocking
-import space.kscience.kmath.misc.PerformancePitfall
-import space.kscience.kmath.samplers.GaussianSampler
 import space.kscience.kmath.misc.toIntExact
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
+import space.kscience.kmath.samplers.GaussianSampler
 import space.kscience.kmath.stat.next
 
+
 public class CMRandomGeneratorWrapper(
     public val factory: (IntArray) -> RandomGenerator,
 ) : org.apache.commons.math3.random.RandomGenerator {

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/transform/Transformations.kt

@@ -10,28 +10,18 @@ import kotlinx.coroutines.flow.Flow
 import kotlinx.coroutines.flow.map
 import org.apache.commons.math3.transform.*
 import space.kscience.kmath.complex.Complex
-import space.kscience.kmath.operations.SuspendBufferTransform
+import space.kscience.kmath.operations.BufferTransform
 import space.kscience.kmath.streaming.chunked
 import space.kscience.kmath.streaming.spread
-import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.VirtualBuffer
-import space.kscience.kmath.structures.asBuffer
-
+import space.kscience.kmath.structures.*
 
 /**
- * Streaming and buffer transformations
+ * Streaming and buffer transformations with Commons-math algorithms
  */
 public object Transformations {
-    private fun Buffer<Complex>.toArray(): Array<org.apache.commons.math3.complex.Complex> =
+    private fun Buffer<Complex>.toCmComplexArray(): Array<org.apache.commons.math3.complex.Complex> =
         Array(size) { org.apache.commons.math3.complex.Complex(get(it).re, get(it).im) }
 
-    private fun Buffer<Double>.asArray() = if (this is DoubleBuffer) {
-        array
-    } else {
-        DoubleArray(size) { i -> get(i) }
-    }
-
     /**
      * Create a virtual buffer on top of array
      */
@@ -43,70 +33,67 @@ public object Transformations {
     public fun fourier(
         normalization: DftNormalization = DftNormalization.STANDARD,
         direction: TransformType = TransformType.FORWARD,
-    ): SuspendBufferTransform<Complex, Complex> = {
-        FastFourierTransformer(normalization).transform(it.toArray(), direction).asBuffer()
+    ): BufferTransform<Complex, Complex> = BufferTransform {
+        FastFourierTransformer(normalization).transform(it.toCmComplexArray(), direction).asBuffer()
     }
 
     public fun realFourier(
         normalization: DftNormalization = DftNormalization.STANDARD,
         direction: TransformType = TransformType.FORWARD,
-    ): SuspendBufferTransform<Double, Complex> = {
-        FastFourierTransformer(normalization).transform(it.asArray(), direction).asBuffer()
+    ): BufferTransform<Double, Complex> = BufferTransform {
+        FastFourierTransformer(normalization).transform(it.toDoubleArray(), direction).asBuffer()
     }
 
     public fun sine(
         normalization: DstNormalization = DstNormalization.STANDARD_DST_I,
         direction: TransformType = TransformType.FORWARD,
-    ): SuspendBufferTransform<Double, Double> = {
-        FastSineTransformer(normalization).transform(it.asArray(), direction).asBuffer()
+    ): BufferTransform<Double, Double> = DoubleBufferTransform {
+        FastSineTransformer(normalization).transform(it.array, direction).asBuffer()
     }
 
     public fun cosine(
         normalization: DctNormalization = DctNormalization.STANDARD_DCT_I,
         direction: TransformType = TransformType.FORWARD,
-    ): SuspendBufferTransform<Double, Double> = {
-        FastCosineTransformer(normalization).transform(it.asArray(), direction).asBuffer()
+    ): BufferTransform<Double, Double> = BufferTransform {
+        FastCosineTransformer(normalization).transform(it.toDoubleArray(), direction).asBuffer()
     }
 
     public fun hadamard(
         direction: TransformType = TransformType.FORWARD,
-    ): SuspendBufferTransform<Double, Double> = {
-        FastHadamardTransformer().transform(it.asArray(), direction).asBuffer()
+    ): BufferTransform<Double, Double> = DoubleBufferTransform {
+        FastHadamardTransformer().transform(it.array, direction).asBuffer()
     }
 }
 
 /**
  * Process given [Flow] with commons-math fft transformation
  */
-@FlowPreview
-public fun Flow<Buffer<Complex>>.FFT(
+public fun Flow<Buffer<Complex>>.fft(
     normalization: DftNormalization = DftNormalization.STANDARD,
     direction: TransformType = TransformType.FORWARD,
 ): Flow<Buffer<Complex>> {
     val transform = Transformations.fourier(normalization, direction)
-    return map { transform(it) }
+    return map(transform::transform)
 }
 
-@FlowPreview
 @JvmName("realFFT")
-public fun Flow<Buffer<Double>>.FFT(
+public fun Flow<Buffer<Double>>.fft(
     normalization: DftNormalization = DftNormalization.STANDARD,
     direction: TransformType = TransformType.FORWARD,
 ): Flow<Buffer<Complex>> {
     val transform = Transformations.realFourier(normalization, direction)
-    return map(transform)
+    return map(transform::transform)
 }
 
 /**
  * Process a continuous flow of real numbers in FFT splitting it in chunks of [bufferSize].
  */
-@FlowPreview
 @JvmName("realFFT")
-public fun Flow<Double>.FFT(
+public fun Flow<Double>.fft(
     bufferSize: Int = Int.MAX_VALUE,
     normalization: DftNormalization = DftNormalization.STANDARD,
     direction: TransformType = TransformType.FORWARD,
-): Flow<Complex> = chunked(bufferSize).FFT(normalization, direction).spread()
+): Flow<Complex> = chunked(bufferSize).fft(normalization, direction).spread()
 
 /**
  * Map a complex flow into real flow by taking real part of each number

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

@@ -13,12 +13,11 @@ import space.kscience.kmath.expressions.Symbol.Companion.x
 import space.kscience.kmath.expressions.Symbol.Companion.y
 import space.kscience.kmath.expressions.chiSquaredExpression
 import space.kscience.kmath.expressions.symbol
-import space.kscience.kmath.operations.map
+import space.kscience.kmath.operations.DoubleBufferOps.Companion.map
 import space.kscience.kmath.optimization.*
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.structures.DoubleBuffer
 import space.kscience.kmath.structures.asBuffer
-import kotlin.math.pow
 import kotlin.test.Test
 
 internal class OptimizeTest {

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

@@ -193,7 +193,6 @@ public object ComplexField :
  * @property re The real part.
  * @property im The imaginary part.
  */
-@OptIn(UnstableKMathAPI::class)
 public data class Complex(val re: Double, val im: Double) {
     public constructor(re: Number, im: Number) : this(re.toDouble(), im.toDouble())
     public constructor(re: Number) : this(re.toDouble(), 0.0)

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

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.complex
 
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.*
@@ -20,6 +21,7 @@ import kotlin.contracts.contract
 public sealed class ComplexFieldOpsND : BufferedFieldOpsND<Complex, ComplexField>(ComplexField.bufferAlgebra),
     ScaleOperations<StructureND<Complex>>, ExtendedFieldOps<StructureND<Complex>>, PowerOperations<StructureND<Complex>> {
 
+    @OptIn(PerformancePitfall::class)
     override fun StructureND<Complex>.toBufferND(): BufferND<Complex> = when (this) {
         is BufferND -> this
         else -> {
@@ -57,7 +59,7 @@ public sealed class ComplexFieldOpsND : BufferedFieldOpsND<Complex, ComplexField
 }
 
 @OptIn(UnstableKMathAPI::class)
-public class ComplexFieldND(override val shape: Shape) :
+public class ComplexFieldND(override val shape: ShapeND) :
     ComplexFieldOpsND(), FieldND<Complex, ComplexField>,
     NumbersAddOps<StructureND<Complex>> {
 
@@ -69,12 +71,12 @@ public class ComplexFieldND(override val shape: Shape) :
 
 public val ComplexField.ndAlgebra: ComplexFieldOpsND get() = ComplexFieldOpsND
 
-public fun ComplexField.ndAlgebra(vararg shape: Int): ComplexFieldND = ComplexFieldND(shape)
+public fun ComplexField.ndAlgebra(vararg shape: Int): ComplexFieldND = ComplexFieldND(ShapeND(shape))
 
 /**
  * Produce a context for n-dimensional operations inside this real field
  */
 public inline fun <R> ComplexField.withNdAlgebra(vararg shape: Int, action: ComplexFieldND.() -> R): R {
     contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
-    return ComplexFieldND(shape).action()
+    return ComplexFieldND(ShapeND(shape)).action()
 }

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

@@ -268,7 +268,7 @@ public open class DSRing<T, A>(
 
     protected fun DS<T, A>.mapData(block: A.(T) -> T): DS<T, A> {
         require(derivativeAlgebra == this@DSRing) { "All derivative operations should be done in the same algebra" }
-        val newData: Buffer<T> = data.map(valueBufferFactory) {
+        val newData: Buffer<T> = data.mapToBuffer(valueBufferFactory) {
             algebra.block(it)
         }
         return DS(newData)
@@ -276,7 +276,7 @@ public open class DSRing<T, A>(
 
     protected fun DS<T, A>.mapDataIndexed(block: (Int, T) -> T): DS<T, A> {
         require(derivativeAlgebra == this@DSRing) { "All derivative operations should be done in the same algebra" }
-        val newData: Buffer<T> = data.mapIndexed(valueBufferFactory, block)
+        val newData: Buffer<T> = data.mapIndexedToBuffer(valueBufferFactory, block)
         return DS(newData)
     }
 

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

@@ -11,6 +11,8 @@ import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.indices
 import kotlin.jvm.JvmName
 
+//TODO move to stat
+
 /**
  * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic
  * differentiation.

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

@@ -6,9 +6,7 @@
 package space.kscience.kmath.linear
 
 import space.kscience.kmath.misc.PerformancePitfall
-import space.kscience.kmath.nd.BufferedRingOpsND
-import space.kscience.kmath.nd.as2D
-import space.kscience.kmath.nd.asND
+import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.*
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.VirtualBuffer
@@ -23,7 +21,7 @@ public class BufferedLinearSpace<T, out A : Ring<T>>(
     private val ndAlgebra = BufferedRingOpsND(bufferAlgebra)
 
     override fun buildMatrix(rows: Int, columns: Int, initializer: A.(i: Int, j: Int) -> T): Matrix<T> =
-        ndAlgebra.structureND(intArrayOf(rows, columns)) { (i, j) -> elementAlgebra.initializer(i, j) }.as2D()
+        ndAlgebra.structureND(ShapeND(rows, columns)) { (i, j) -> elementAlgebra.initializer(i, j) }.as2D()
 
     override fun buildVector(size: Int, initializer: A.(Int) -> T): Point<T> =
         bufferAlgebra.buffer(size) { elementAlgebra.initializer(it) }

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

@@ -6,9 +6,7 @@
 package space.kscience.kmath.linear
 
 import space.kscience.kmath.misc.PerformancePitfall
-import space.kscience.kmath.nd.DoubleFieldOpsND
-import space.kscience.kmath.nd.as2D
-import space.kscience.kmath.nd.asND
+import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.DoubleBufferOps
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.invoke
@@ -23,7 +21,7 @@ public object DoubleLinearSpace : LinearSpace<Double, DoubleField> {
         rows: Int,
         columns: Int,
         initializer: DoubleField.(i: Int, j: Int) -> Double
-    ): Matrix<Double> = DoubleFieldOpsND.structureND(intArrayOf(rows, columns)) { (i, j) ->
+    ): Matrix<Double> = DoubleFieldOpsND.structureND(ShapeND(rows, columns)) { (i, j) ->
         DoubleField.initializer(i, j)
     }.as2D()
 

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

@@ -5,6 +5,9 @@
 
 package space.kscience.kmath.linear
 
+import space.kscience.kmath.nd.ShapeND
+
+
 /**
  * The matrix where each element is evaluated each time when is being accessed.
  *
@@ -16,7 +19,7 @@ public class VirtualMatrix<out T : Any>(
     public val generator: (i: Int, j: Int) -> T,
 ) : Matrix<T> {
 
-    override val shape: IntArray get() = intArrayOf(rowNum, colNum)
+    override val shape: ShapeND get() = ShapeND(rowNum, colNum)
 
     override operator fun get(i: Int, j: Int): T = generator(i, j)
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/annotations.kt

@@ -29,3 +29,16 @@ public annotation class UnstableKMathAPI
 public annotation class PerformancePitfall(
     val message: String = "Potential performance problem",
 )
+
+/**
+ * Marks API that is public, but should not be used without clear understanding what it does.
+ */
+@MustBeDocumented
+@Retention(value = AnnotationRetention.BINARY)
+@RequiresOptIn(
+    "This API is unsafe and should be used carefully",
+    RequiresOptIn.Level.ERROR,
+)
+public annotation class UnsafeKMathAPI(
+    val message: String = "Unsafe API",
+)

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

@@ -25,7 +25,7 @@ public interface AlgebraND<T, out C : Algebra<T>>: Algebra<StructureND<T>> {
     /**
      * Produces a new [StructureND] using given initializer function.
      */
-    public fun structureND(shape: Shape, initializer: C.(IntArray) -> T): StructureND<T>
+    public fun structureND(shape: ShapeND, initializer: C.(IntArray) -> T): StructureND<T>
 
     /**
      * Maps elements from one structure to another one by applying [transform] to them.

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

@@ -12,11 +12,11 @@ import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.*
 
 public interface BufferAlgebraND<T, out A : Algebra<T>> : AlgebraND<T, A> {
-    public val indexerBuilder: (IntArray) -> ShapeIndexer
+    public val indexerBuilder: (ShapeND) -> ShapeIndexer
     public val bufferAlgebra: BufferAlgebra<T, A>
     override val elementAlgebra: A get() = bufferAlgebra.elementAlgebra
 
-    override fun structureND(shape: Shape, initializer: A.(IntArray) -> T): BufferND<T> {
+    override fun structureND(shape: ShapeND, initializer: A.(IntArray) -> T): BufferND<T> {
         val indexer = indexerBuilder(shape)
         return BufferND(
             indexer,
@@ -26,6 +26,7 @@ public interface BufferAlgebraND<T, out A : Algebra<T>> : AlgebraND<T, A> {
         )
     }
 
+    @OptIn(PerformancePitfall::class)
     public fun StructureND<T>.toBufferND(): BufferND<T> = when (this) {
         is BufferND -> this
         else -> {
@@ -46,7 +47,7 @@ public interface BufferAlgebraND<T, out A : Algebra<T>> : AlgebraND<T, A> {
         zipInline(left.toBufferND(), right.toBufferND(), transform)
 
     public companion object {
-        public val defaultIndexerBuilder: (IntArray) -> ShapeIndexer = ::Strides
+        public val defaultIndexerBuilder: (ShapeND) -> ShapeIndexer = ::Strides
     }
 }
 
@@ -98,24 +99,24 @@ internal inline fun <T, A : Algebra<T>> BufferAlgebraND<T, A>.zipInline(
 @OptIn(PerformancePitfall::class)
 public open class BufferedGroupNDOps<T, out A : Group<T>>(
     override val bufferAlgebra: BufferAlgebra<T, A>,
-    override val indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
+    override val indexerBuilder: (ShapeND) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
 ) : GroupOpsND<T, A>, BufferAlgebraND<T, A> {
     override fun StructureND<T>.unaryMinus(): StructureND<T> = map { -it }
 }
 
 public open class BufferedRingOpsND<T, out A : Ring<T>>(
     bufferAlgebra: BufferAlgebra<T, A>,
-    indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
+    indexerBuilder: (ShapeND) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
 ) : BufferedGroupNDOps<T, A>(bufferAlgebra, indexerBuilder), RingOpsND<T, A>
 
 public open class BufferedFieldOpsND<T, out A : Field<T>>(
     bufferAlgebra: BufferAlgebra<T, A>,
-    indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
+    indexerBuilder: (ShapeND) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
 ) : BufferedRingOpsND<T, A>(bufferAlgebra, indexerBuilder), FieldOpsND<T, A> {
 
     public constructor(
         elementAlgebra: A,
-        indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
+        indexerBuilder: (ShapeND) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
     ) : this(BufferFieldOps(elementAlgebra), indexerBuilder)
 
     @OptIn(PerformancePitfall::class)
@@ -130,7 +131,7 @@ public val <T, A : Field<T>> BufferAlgebra<T, A>.nd: BufferedFieldOpsND<T, A> ge
 public fun <T, A : Algebra<T>> BufferAlgebraND<T, A>.structureND(
     vararg shape: Int,
     initializer: A.(IntArray) -> T,
-): BufferND<T> = structureND(shape, initializer)
+): BufferND<T> = structureND(ShapeND(shape), initializer)
 
 public fun <T, EA : Algebra<T>, A> A.structureND(
     initializer: EA.(IntArray) -> T,

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

@@ -5,10 +5,9 @@
 
 package space.kscience.kmath.nd
 
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.BufferFactory
 import space.kscience.kmath.structures.MutableBuffer
-import space.kscience.kmath.structures.MutableBufferFactory
 
 /**
  * Represents [StructureND] over [Buffer].
@@ -22,32 +21,33 @@ public open class BufferND<out T>(
     public open val buffer: Buffer<T>,
 ) : StructureND<T> {
 
+    @PerformancePitfall
     override operator fun get(index: IntArray): T = buffer[indices.offset(index)]
 
-    override val shape: IntArray get() = indices.shape
+    override val shape: ShapeND get() = indices.shape
 
     override fun toString(): String = StructureND.toString(this)
 }
 
-/**
- * Transform structure to a new structure using provided [BufferFactory] and optimizing if argument is [BufferND]
- */
-public inline fun <T, R : Any> StructureND<T>.mapToBuffer(
-    factory: BufferFactory<R>,
-    crossinline transform: (T) -> R,
-): BufferND<R> = if (this is BufferND<T>)
-    BufferND(this.indices, factory.invoke(indices.linearSize) { transform(buffer[it]) })
-else {
-    val strides = DefaultStrides(shape)
-    BufferND(strides, factory.invoke(strides.linearSize) { transform(get(strides.index(it))) })
-}
-
-/**
- * Transform structure to a new structure using inferred [BufferFactory]
- */
-public inline fun <T, reified R : Any> StructureND<T>.mapToBuffer(
-    crossinline transform: (T) -> R,
-): BufferND<R> = mapToBuffer(Buffer.Companion::auto, transform)
+///**
+// * Transform structure to a new structure using provided [BufferFactory] and optimizing if argument is [BufferND]
+// */
+//public inline fun <T, R : Any> StructureND<T>.mapToBuffer(
+//    factory: BufferFactory<R>,
+//    crossinline transform: (T) -> R,
+//): BufferND<R> = if (this is BufferND<T>)
+//    BufferND(this.indices, factory.invoke(indices.linearSize) { transform(buffer[it]) })
+//else {
+//    val strides = ColumnStrides(shape)
+//    BufferND(strides, factory.invoke(strides.linearSize) { transform(get(strides.index(it))) })
+//}
+//
+///**
+// * Transform structure to a new structure using inferred [BufferFactory]
+// */
+//public inline fun <T, reified R : Any> StructureND<T>.mapToBuffer(
+//    crossinline transform: (T) -> R,
+//): BufferND<R> = mapToBuffer(Buffer.Companion::auto, transform)
 
 /**
  * Represents [MutableStructureND] over [MutableBuffer].
@@ -60,22 +60,24 @@ public open class MutableBufferND<T>(
     strides: ShapeIndexer,
     override val buffer: MutableBuffer<T>,
 ) : MutableStructureND<T>, BufferND<T>(strides, buffer) {
+
+    @PerformancePitfall
     override fun set(index: IntArray, value: T) {
         buffer[indices.offset(index)] = value
     }
 }
 
-/**
- * Transform structure to a new structure using provided [MutableBufferFactory] and optimizing if argument is [MutableBufferND]
- */
-public inline fun <T, reified R : Any> MutableStructureND<T>.mapToMutableBuffer(
-    factory: MutableBufferFactory<R> = MutableBufferFactory(MutableBuffer.Companion::auto),
-    crossinline transform: (T) -> R,
-): MutableBufferND<R> {
-    return if (this is MutableBufferND<T>)
-        MutableBufferND(this.indices, factory.invoke(indices.linearSize) { transform(buffer[it]) })
-    else {
-        val strides = DefaultStrides(shape)
-        MutableBufferND(strides, factory.invoke(strides.linearSize) { transform(get(strides.index(it))) })
-    }
-}
+///**
+// * Transform structure to a new structure using provided [MutableBufferFactory] and optimizing if argument is [MutableBufferND]
+// */
+//public inline fun <T, reified R : Any> MutableStructureND<T>.mapToMutableBuffer(
+//    factory: MutableBufferFactory<R> = MutableBufferFactory(MutableBuffer.Companion::auto),
+//    crossinline transform: (T) -> R,
+//): MutableBufferND<R> {
+//    return if (this is MutableBufferND<T>)
+//        MutableBufferND(this.indices, factory.invoke(indices.linearSize) { transform(buffer[it]) })
+//    else {
+//        val strides = ColumnStrides(shape)
+//        MutableBufferND(strides, factory.invoke(strides.linearSize) { transform(get(strides.index(it))) })
+//    }
+//}

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

@@ -14,15 +14,25 @@ import kotlin.contracts.contract
 import kotlin.math.pow
 import kotlin.math.pow as kpow
 
+/**
+ * A simple mutable [StructureND] of doubles
+ */
 public class DoubleBufferND(
     indexes: ShapeIndexer,
     override val buffer: DoubleBuffer,
-) : MutableBufferND<Double>(indexes, buffer)
+) : MutableBufferND<Double>(indexes, buffer), MutableStructureNDOfDouble{
+    override fun getDouble(index: IntArray): Double = buffer[indices.offset(index)]
+
+    override fun setDouble(index: IntArray, value: Double) {
+        buffer[indices.offset(index)] = value
+    }
+}
 
 
 public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(DoubleField.bufferAlgebra),
     ScaleOperations<StructureND<Double>>, ExtendedFieldOps<StructureND<Double>> {
 
+    @OptIn(PerformancePitfall::class)
     override fun StructureND<Double>.toBufferND(): DoubleBufferND = when (this) {
         is DoubleBufferND -> this
         else -> {
@@ -64,7 +74,7 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
         transform: DoubleField.(Double, Double) -> Double,
     ): BufferND<Double> = zipInline(left.toBufferND(), right.toBufferND()) { l, r -> DoubleField.transform(l, r) }
 
-    override fun structureND(shape: Shape, initializer: DoubleField.(IntArray) -> Double): DoubleBufferND {
+    override fun structureND(shape: ShapeND, initializer: DoubleField.(IntArray) -> Double): DoubleBufferND {
         val indexer = indexerBuilder(shape)
         return DoubleBufferND(
             indexer,
@@ -179,7 +189,7 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
 }
 
 @OptIn(UnstableKMathAPI::class)
-public class DoubleFieldND(override val shape: Shape) :
+public class DoubleFieldND(override val shape: ShapeND) :
     DoubleFieldOpsND(), FieldND<Double, DoubleField>, NumbersAddOps<StructureND<Double>>,
     ExtendedField<StructureND<Double>> {
 
@@ -221,7 +231,8 @@ public class DoubleFieldND(override val shape: Shape) :
 
 public val DoubleField.ndAlgebra: DoubleFieldOpsND get() = DoubleFieldOpsND
 
-public fun DoubleField.ndAlgebra(vararg shape: Int): DoubleFieldND = DoubleFieldND(shape)
+public fun DoubleField.ndAlgebra(vararg shape: Int): DoubleFieldND = DoubleFieldND(ShapeND(shape))
+public fun DoubleField.ndAlgebra(shape: ShapeND): DoubleFieldND = DoubleFieldND(shape)
 
 /**
  * Produce a context for n-dimensional operations inside this real field
@@ -229,5 +240,5 @@ public fun DoubleField.ndAlgebra(vararg shape: Int): DoubleFieldND = DoubleField
 @UnstableKMathAPI
 public inline fun <R> DoubleField.withNdAlgebra(vararg shape: Int, action: DoubleFieldND.() -> R): R {
     contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
-    return DoubleFieldND(shape).run(action)
+    return DoubleFieldND(ShapeND(shape)).run(action)
 }

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

@@ -20,7 +20,7 @@ public class IntBufferND(
 
 public sealed class IntRingOpsND : BufferedRingOpsND<Int, IntRing>(IntRing.bufferAlgebra) {
 
-    override fun structureND(shape: Shape, initializer: IntRing.(IntArray) -> Int): IntBufferND {
+    override fun structureND(shape: ShapeND, initializer: IntRing.(IntArray) -> Int): IntBufferND {
         val indexer = indexerBuilder(shape)
         return IntBufferND(
             indexer,
@@ -35,7 +35,7 @@ public sealed class IntRingOpsND : BufferedRingOpsND<Int, IntRing>(IntRing.buffe
 
 @OptIn(UnstableKMathAPI::class)
 public class IntRingND(
-    override val shape: Shape
+    override val shape: ShapeND
 ) : IntRingOpsND(), RingND<Int, IntRing>, NumbersAddOps<StructureND<Int>> {
 
     override fun number(value: Number): BufferND<Int> {
@@ -46,5 +46,5 @@ public class IntRingND(
 
 public inline fun <R> IntRing.withNdAlgebra(vararg shape: Int, action: IntRingND.() -> R): R {
     contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
-    return IntRingND(shape).run(action)
+    return IntRingND(ShapeND(shape)).run(action)
 }

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

@@ -0,0 +1,50 @@
+/*
+ * Copyright 2018-2022 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.nd
+
+import space.kscience.kmath.misc.PerformancePitfall
+
+
+public class PermutedStructureND<T>(
+    public val origin: StructureND<T>,
+    public val permutation: (IntArray) -> IntArray,
+) : StructureND<T> {
+
+    override val shape: ShapeND
+        get() = origin.shape
+
+    @OptIn(PerformancePitfall::class)
+    override fun get(index: IntArray): T {
+        return origin[permutation(index)]
+    }
+}
+
+public fun <T> StructureND<T>.permute(
+    permutation: (IntArray) -> IntArray,
+): PermutedStructureND<T> = PermutedStructureND(this, permutation)
+
+public class PermutedMutableStructureND<T>(
+    public val origin: MutableStructureND<T>,
+    override val shape: ShapeND = origin.shape,
+    public val permutation: (IntArray) -> IntArray,
+) : MutableStructureND<T> {
+
+
+    @OptIn(PerformancePitfall::class)
+    override fun set(index: IntArray, value: T) {
+        origin[permutation(index)] = value
+    }
+
+    @OptIn(PerformancePitfall::class)
+    override fun get(index: IntArray): T {
+        return origin[permutation(index)]
+    }
+}
+
+public fun <T> MutableStructureND<T>.permute(
+    newShape: ShapeND = shape,
+    permutation: (IntArray) -> IntArray,
+): PermutedMutableStructureND<T> = PermutedMutableStructureND(this, newShape, permutation)

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

@@ -1,35 +0,0 @@
-/*
- * Copyright 2018-2022 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.nd
-
-/**
- * An exception is thrown when the expected and actual shape of NDArray differ.
- *
- * @property expected the expected shape.
- * @property actual the actual shape.
- */
-public class ShapeMismatchException(public val expected: IntArray, public val actual: IntArray) :
-    RuntimeException("Shape ${actual.contentToString()} doesn't fit in expected shape ${expected.contentToString()}.")
-
-public class IndexOutOfShapeException(public val shape: Shape, public val index: IntArray) :
-    RuntimeException("Index ${index.contentToString()} is out of shape ${shape.contentToString()}")
-
-public typealias Shape = IntArray
-
-public fun Shape(shapeFirst: Int, vararg shapeRest: Int): Shape = intArrayOf(shapeFirst, *shapeRest)
-
-public interface WithShape {
-    public val shape: Shape
-
-    public val indices: ShapeIndexer get() = DefaultStrides(shape)
-}
-
-internal fun requireIndexInShape(index: IntArray, shape: Shape) {
-    if (index.size != shape.size) throw IndexOutOfShapeException(index, shape)
-    shape.forEachIndexed { axis, axisShape ->
-        if (index[axis] !in 0 until axisShape) throw IndexOutOfShapeException(index, shape)
-    }
-}

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

@@ -1,127 +0,0 @@
-/*
- * Copyright 2018-2022 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.nd
-
-import kotlin.native.concurrent.ThreadLocal
-
-/**
- * A converter from linear index to multivariate index
- */
-public interface ShapeIndexer : Iterable<IntArray> {
-    public val shape: Shape
-
-    /**
-     * Get linear index from multidimensional index
-     */
-    public fun offset(index: IntArray): Int
-
-    /**
-     * Get multidimensional from linear
-     */
-    public fun index(offset: Int): IntArray
-
-    /**
-     * The size of linear buffer to accommodate all elements of ND-structure corresponding to strides
-     */
-    public val linearSize: Int
-
-    // TODO introduce a fast way to calculate index of the next element?
-
-    /**
-     * Iterate over ND indices in a natural order
-     */
-    public fun asSequence(): Sequence<IntArray>
-
-    override fun iterator(): Iterator<IntArray> = asSequence().iterator()
-
-    override fun equals(other: Any?): Boolean
-    override fun hashCode(): Int
-}
-
-/**
- * Linear transformation of indexes
- */
-public abstract class Strides : ShapeIndexer {
-    /**
-     * Array strides
-     */
-    public abstract val strides: IntArray
-
-    public override fun offset(index: IntArray): Int = index.mapIndexed { i, value ->
-        if (value < 0 || value >= shape[i]) throw IndexOutOfBoundsException("Index $value out of shape bounds: (0,${this.shape[i]})")
-        value * strides[i]
-    }.sum()
-
-    // TODO introduce a fast way to calculate index of the next element?
-
-    /**
-     * Iterate over ND indices in a natural order
-     */
-    public override fun asSequence(): Sequence<IntArray> = (0 until linearSize).asSequence().map(::index)
-}
-
-/**
- * Simple implementation of [Strides].
- */
-public class DefaultStrides(override val shape: IntArray) : Strides() {
-    override val linearSize: Int get() = strides[shape.size]
-
-    /**
-     * Strides for memory access
-     */
-    override val strides: IntArray by lazy {
-        sequence {
-            var current = 1
-            yield(1)
-
-            shape.forEach {
-                current *= it
-                yield(current)
-            }
-        }.toList().toIntArray()
-    }
-
-    override fun index(offset: Int): IntArray {
-        val res = IntArray(shape.size)
-        var current = offset
-        var strideIndex = strides.size - 2
-
-        while (strideIndex >= 0) {
-            res[strideIndex] = (current / strides[strideIndex])
-            current %= strides[strideIndex]
-            strideIndex--
-        }
-
-        return res
-    }
-
-    override fun equals(other: Any?): Boolean {
-        if (this === other) return true
-        if (other !is DefaultStrides) return false
-        if (!shape.contentEquals(other.shape)) return false
-        return true
-    }
-
-    override fun hashCode(): Int = shape.contentHashCode()
-
-
-    public companion object {
-        /**
-         * Cached builder for default strides
-         */
-        @Deprecated("Replace by Strides(shape)")
-        public operator fun invoke(shape: IntArray): Strides =
-            defaultStridesCache.getOrPut(shape) { DefaultStrides(shape) }
-    }
-}
-
-@ThreadLocal
-private val defaultStridesCache = HashMap<IntArray, Strides>()
-
-/**
- * Cached builder for default strides
- */
-public fun Strides(shape: IntArray): Strides = defaultStridesCache.getOrPut(shape) { DefaultStrides(shape) }

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

@@ -0,0 +1,174 @@
+/*
+ * Copyright 2018-2022 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.nd
+
+import kotlin.math.max
+import kotlin.native.concurrent.ThreadLocal
+
+/**
+ * A converter from linear index to multivariate index
+ */
+public interface ShapeIndexer : Iterable<IntArray> {
+    public val shape: ShapeND
+
+    /**
+     * Get linear index from multidimensional index
+     */
+    public fun offset(index: IntArray): Int
+
+    /**
+     * Get multidimensional from linear
+     */
+    public fun index(offset: Int): IntArray
+
+    /**
+     * The size of linear buffer to accommodate all elements of ND-structure corresponding to strides
+     */
+    public val linearSize: Int
+
+    // TODO introduce a fast way to calculate index of the next element?
+
+    /**
+     * Iterate over ND indices in a natural order
+     */
+    public fun asSequence(): Sequence<IntArray>
+
+    override fun iterator(): Iterator<IntArray> = asSequence().iterator()
+
+    override fun equals(other: Any?): Boolean
+    override fun hashCode(): Int
+}
+
+/**
+ * Linear transformation of indexes
+ */
+public abstract class Strides : ShapeIndexer {
+    /**
+     * Array strides
+     */
+    internal abstract val strides: IntArray
+
+    public override fun offset(index: IntArray): Int {
+        var res = 0
+        index.forEachIndexed { i, value ->
+            if (value !in 0 until shape[i]) throw IndexOutOfBoundsException("Index $value out of shape bounds: (0, ${this.shape[i]})")
+            res += value * strides[i]
+
+        }
+        return res
+    }
+
+    // TODO introduce a fast way to calculate index of the next element?
+
+    /**
+     * Iterate over ND indices in a natural order
+     */
+    public override fun asSequence(): Sequence<IntArray> = (0 until linearSize).asSequence().map(::index)
+
+}
+
+/**
+ * Column-first [Strides]. Columns are represented as continuous arrays
+ */
+public class ColumnStrides(override val shape: ShapeND) : Strides() {
+    override val linearSize: Int get() = strides[shape.size]
+
+    /**
+     * Strides for memory access
+     */
+    override val strides: IntArray = sequence {
+        var current = 1
+        yield(1)
+
+        shape.forEach {
+            current *= it
+            yield(current)
+        }
+    }.toList().toIntArray()
+
+    override fun index(offset: Int): IntArray {
+        val res = IntArray(shape.size)
+        var current = offset
+        var strideIndex = strides.size - 2
+
+        while (strideIndex >= 0) {
+            res[strideIndex] = (current / strides[strideIndex])
+            current %= strides[strideIndex]
+            strideIndex--
+        }
+
+        return res
+    }
+
+    override fun equals(other: Any?): Boolean {
+        if (this === other) return true
+        if (other !is ColumnStrides) return false
+        return shape.contentEquals(other.shape)
+    }
+
+    override fun hashCode(): Int = shape.contentHashCode()
+
+
+    public companion object
+}
+
+/**
+ * This [Strides] implementation follows the last dimension first convention
+ * For more information: https://numpy.org/doc/stable/reference/generated/numpy.ndarray.strides.html
+ *
+ * @param shape the shape of the tensor.
+ */
+public class RowStrides(override val shape: ShapeND) : Strides() {
+
+    override val strides: IntArray = run {
+        val nDim = shape.size
+        val res = IntArray(nDim)
+        if (nDim == 0) return@run res
+
+        var current = nDim - 1
+        res[current] = 1
+
+        while (current > 0) {
+            res[current - 1] = max(1, shape[current]) * res[current]
+            current--
+        }
+        res
+    }
+
+    override fun index(offset: Int): IntArray {
+        val res = IntArray(shape.size)
+        var current = offset
+        var strideIndex = 0
+
+        while (strideIndex < shape.size) {
+            res[strideIndex] = (current / strides[strideIndex])
+            current %= strides[strideIndex]
+            strideIndex++
+        }
+        return res
+    }
+
+    override val linearSize: Int get() = shape.linearSize
+
+    override fun equals(other: Any?): Boolean {
+        if (this === other) return true
+        if (other !is RowStrides) return false
+        return shape.contentEquals(other.shape)
+    }
+
+    override fun hashCode(): Int = shape.contentHashCode()
+
+    public companion object
+
+}
+
+@ThreadLocal
+private val defaultStridesCache = HashMap<ShapeND, Strides>()
+
+/**
+ * Cached builder for default strides
+ */
+public fun Strides(shape: ShapeND): Strides = defaultStridesCache.getOrPut(shape) { RowStrides(shape) }

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

@@ -0,0 +1,102 @@
+/*
+ * Copyright 2018-2022 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.nd
+
+import space.kscience.kmath.misc.UnsafeKMathAPI
+import kotlin.jvm.JvmInline
+
+/**
+ * A read-only ND shape
+ */
+@JvmInline
+public value class ShapeND(@PublishedApi internal val array: IntArray) {
+    public val size: Int get() = array.size
+    public operator fun get(index: Int): Int = array[index]
+    override fun toString(): String = array.contentToString()
+}
+
+public inline fun ShapeND.forEach(block: (value: Int) -> Unit): Unit = array.forEach(block)
+
+public inline fun ShapeND.forEachIndexed(block: (index: Int, value: Int) -> Unit): Unit = array.forEachIndexed(block)
+
+public infix fun ShapeND.contentEquals(other: ShapeND): Boolean = array.contentEquals(other.array)
+
+public fun ShapeND.contentHashCode(): Int = array.contentHashCode()
+
+public val ShapeND.indices: IntRange get() = array.indices
+public val ShapeND.linearSize: Int get() = array.reduce(Int::times)
+
+public fun ShapeND.slice(range: IntRange): ShapeND = ShapeND(array.sliceArray(range))
+
+public fun ShapeND.last(): Int = array.last()
+
+/**
+ * A shape including last [n] dimensions of this shape
+ */
+public fun ShapeND.last(n: Int): ShapeND = ShapeND(array.copyOfRange(size - n, size))
+
+public fun ShapeND.first(): Int = array.first()
+
+/**
+ * A shape including first [n] dimensions of this shape
+ */
+public fun ShapeND.first(n: Int): ShapeND = ShapeND(array.copyOfRange(0, n))
+
+public operator fun ShapeND.plus(add: IntArray): ShapeND = ShapeND(array + add)
+
+public operator fun ShapeND.plus(add: ShapeND): ShapeND = ShapeND(array + add.array)
+
+public fun ShapeND.isEmpty(): Boolean = size == 0
+public fun ShapeND.isNotEmpty(): Boolean = size > 0
+
+public fun ShapeND.transposed(i: Int, j: Int): ShapeND = ShapeND(array.copyOf().apply {
+    val ith = get(i)
+    val jth = get(j)
+    set(i, jth)
+    set(j, ith)
+})
+
+public operator fun ShapeND.component1(): Int = get(0)
+public operator fun ShapeND.component2(): Int = get(1)
+public operator fun ShapeND.component3(): Int = get(2)
+
+/**
+ * Convert to array with protective copy
+ */
+public fun ShapeND.toArray(): IntArray = array.copyOf()
+
+@UnsafeKMathAPI
+public fun ShapeND.asArray(): IntArray = array
+
+public fun ShapeND.asList(): List<Int> = array.asList()
+
+
+/**
+ * An exception is thrown when the expected and actual shape of NDArray differ.
+ *
+ * @property expected the expected shape.
+ * @property actual the actual shape.
+ */
+public class ShapeMismatchException(public val expected: ShapeND, public val actual: ShapeND) :
+    RuntimeException("Shape $actual doesn't fit in expected shape ${expected}.")
+
+public class IndexOutOfShapeException(public val shape: ShapeND, public val index: IntArray) :
+    RuntimeException("Index ${index.contentToString()} is out of shape ${shape}")
+
+public fun ShapeND(shapeFirst: Int, vararg shapeRest: Int): ShapeND = ShapeND(intArrayOf(shapeFirst, *shapeRest))
+
+public interface WithShape {
+    public val shape: ShapeND
+
+    public val indices: ShapeIndexer get() = ColumnStrides(shape)
+}
+
+internal fun requireIndexInShape(index: IntArray, shape: ShapeND) {
+    if (index.size != shape.size) throw IndexOutOfShapeException(shape, index)
+    shape.forEachIndexed { axis, axisShape ->
+        if (index[axis] !in 0 until axisShape) throw IndexOutOfShapeException(shape, index)
+    }
+}

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

@@ -18,7 +18,7 @@ public sealed class ShortRingOpsND : BufferedRingOpsND<Short, ShortRing>(ShortRi
 
 @OptIn(UnstableKMathAPI::class)
 public class ShortRingND(
-    override val shape: Shape
+    override val shape: ShapeND
 ) : ShortRingOpsND(), RingND<Short, ShortRing>, NumbersAddOps<StructureND<Short>> {
 
     override fun number(value: Number): BufferND<Short> {
@@ -30,5 +30,5 @@ public class ShortRingND(
 
 public inline fun <R> ShortRing.withNdAlgebra(vararg shape: Int, action: ShortRingND.() -> R): R {
     contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
-    return ShortRingND(shape).run(action)
+    return ShortRingND(ShapeND(shape)).run(action)
 }

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

@@ -18,6 +18,7 @@ import kotlin.jvm.JvmInline
 public interface Structure1D<out T> : StructureND<T>, Buffer<T> {
     override val dimension: Int get() = 1
 
+    @PerformancePitfall
     override operator fun get(index: IntArray): T {
         require(index.size == 1) { "Index dimension mismatch. Expected 1 but found ${index.size}" }
         return get(index[0])
@@ -32,6 +33,8 @@ public interface Structure1D<out T> : StructureND<T>, Buffer<T> {
  * A mutable structure that is guaranteed to be one-dimensional
  */
 public interface MutableStructure1D<T> : Structure1D<T>, MutableStructureND<T>, MutableBuffer<T> {
+
+    @PerformancePitfall
     override operator fun set(index: IntArray, value: T) {
         require(index.size == 1) { "Index dimension mismatch. Expected 1 but found ${index.size}" }
         set(index[0], value)
@@ -43,9 +46,10 @@ public interface MutableStructure1D<T> : Structure1D<T>, MutableStructureND<T>,
  */
 @JvmInline
 private value class Structure1DWrapper<out T>(val structure: StructureND<T>) : Structure1D<T> {
-    override val shape: IntArray get() = structure.shape
+    override val shape: ShapeND get() = structure.shape
     override val size: Int get() = structure.shape[0]
 
+    @PerformancePitfall
     override operator fun get(index: Int): T = structure[index]
 
     @PerformancePitfall
@@ -56,13 +60,16 @@ private value class Structure1DWrapper<out T>(val structure: StructureND<T>) : S
  * A 1D wrapper for a mutable nd-structure
  */
 private class MutableStructure1DWrapper<T>(val structure: MutableStructureND<T>) : MutableStructure1D<T> {
-    override val shape: IntArray get() = structure.shape
+    override val shape: ShapeND get() = structure.shape
     override val size: Int get() = structure.shape[0]
 
     @PerformancePitfall
     override fun elements(): Sequence<Pair<IntArray, T>> = structure.elements()
 
+    @PerformancePitfall
     override fun get(index: Int): T = structure[index]
+
+    @PerformancePitfall
     override fun set(index: Int, value: T) {
         structure[intArrayOf(index)] = value
     }
@@ -83,7 +90,7 @@ private class MutableStructure1DWrapper<T>(val structure: MutableStructureND<T>)
  */
 @JvmInline
 private value class Buffer1DWrapper<out T>(val buffer: Buffer<T>) : Structure1D<T> {
-    override val shape: IntArray get() = intArrayOf(buffer.size)
+    override val shape: ShapeND get() = ShapeND(buffer.size)
     override val size: Int get() = buffer.size
 
     @PerformancePitfall
@@ -95,7 +102,7 @@ private value class Buffer1DWrapper<out T>(val buffer: Buffer<T>) : Structure1D<
 }
 
 internal class MutableBuffer1DWrapper<T>(val buffer: MutableBuffer<T>) : MutableStructure1D<T> {
-    override val shape: IntArray get() = intArrayOf(buffer.size)
+    override val shape: ShapeND get() = ShapeND(buffer.size)
     override val size: Int get() = buffer.size
 
     @PerformancePitfall

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

@@ -7,6 +7,7 @@ package space.kscience.kmath.nd
 
 import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.MutableBuffer
 import space.kscience.kmath.structures.MutableListBuffer
 import space.kscience.kmath.structures.VirtualBuffer
 import kotlin.jvm.JvmInline
@@ -28,7 +29,7 @@ public interface Structure2D<out T> : StructureND<T> {
      */
     public val colNum: Int
 
-    override val shape: IntArray get() = intArrayOf(rowNum, colNum)
+    override val shape: ShapeND get() = ShapeND(rowNum, colNum)
 
     /**
      * The buffer of rows of this structure. It gets elements from the structure dynamically.
@@ -53,6 +54,7 @@ public interface Structure2D<out T> : StructureND<T> {
      */
     public operator fun get(i: Int, j: Int): T
 
+    @PerformancePitfall
     override operator fun get(index: IntArray): T {
         require(index.size == 2) { "Index dimension mismatch. Expected 2 but found ${index.size}" }
         return get(index[0], index[1])
@@ -84,15 +86,15 @@ public interface MutableStructure2D<T> : Structure2D<T>, MutableStructureND<T> {
      * The buffer of rows of this structure. It gets elements from the structure dynamically.
      */
     @PerformancePitfall
-    override val rows: List<MutableStructure1D<T>>
-        get() = List(rowNum) { i -> MutableBuffer1DWrapper(MutableListBuffer(colNum) { j -> get(i, j) }) }
+    override val rows: List<MutableBuffer<T>>
+        get() = List(rowNum) { i -> MutableListBuffer(colNum) { j -> get(i, j) } }
 
     /**
      * The buffer of columns of this structure. It gets elements from the structure dynamically.
      */
     @PerformancePitfall
-    override val columns: List<MutableStructure1D<T>>
-        get() = List(colNum) { j -> MutableBuffer1DWrapper(MutableListBuffer(rowNum) { i -> get(i, j) }) }
+    override val columns: List<MutableBuffer<T>>
+        get() = List(colNum) { j -> MutableListBuffer(rowNum) { i -> get(i, j) } }
 }
 
 /**
@@ -100,11 +102,12 @@ public interface MutableStructure2D<T> : Structure2D<T>, MutableStructureND<T> {
  */
 @JvmInline
 private value class Structure2DWrapper<out T>(val structure: StructureND<T>) : Structure2D<T> {
-    override val shape: Shape get() = structure.shape
+    override val shape: ShapeND get() = structure.shape
 
     override val rowNum: Int get() = shape[0]
     override val colNum: Int get() = shape[1]
 
+    @PerformancePitfall
     override operator fun get(i: Int, j: Int): T = structure[i, j]
 
     override fun <F : StructureFeature> getFeature(type: KClass<out F>): F? = structure.getFeature(type)
@@ -117,17 +120,20 @@ private value class Structure2DWrapper<out T>(val structure: StructureND<T>) : S
  * A 2D wrapper for a mutable nd-structure
  */
 private class MutableStructure2DWrapper<T>(val structure: MutableStructureND<T>) : MutableStructure2D<T> {
-    override val shape: Shape get() = structure.shape
+    override val shape: ShapeND get() = structure.shape
 
     override val rowNum: Int get() = shape[0]
     override val colNum: Int get() = shape[1]
 
+    @PerformancePitfall
     override operator fun get(i: Int, j: Int): T = structure[i, j]
 
+    @PerformancePitfall
     override fun set(index: IntArray, value: T) {
         structure[index] = value
     }
 
+    @PerformancePitfall
     override operator fun set(i: Int, j: Int, value: T) {
         structure[intArrayOf(i, j)] = value
     }

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

@@ -33,7 +33,7 @@ public interface StructureND<out T> : Featured<StructureFeature>, WithShape {
      * The shape of structure i.e., non-empty sequence of non-negative integers that specify sizes of dimensions of
      * this structure.
      */
-    override val shape: Shape
+    override val shape: ShapeND
 
     /**
      * The count of dimensions in this structure. It should be equal to size of [shape].
@@ -46,6 +46,7 @@ public interface StructureND<out T> : Featured<StructureFeature>, WithShape {
      * @param index the indices.
      * @return the value.
      */
+    @PerformancePitfall
     public operator fun get(index: IntArray): T
 
     /**
@@ -97,6 +98,7 @@ public interface StructureND<out T> : Featured<StructureFeature>, WithShape {
         /**
          * Debug output to string
          */
+        @OptIn(PerformancePitfall::class)
         public fun toString(structure: StructureND<*>): String {
             val bufferRepr: String = when (structure.shape.size) {
                 1 -> (0 until structure.shape[0]).map { structure[it] }
@@ -116,7 +118,7 @@ public interface StructureND<out T> : Featured<StructureFeature>, WithShape {
             }
             val className = structure::class.simpleName ?: "StructureND"
 
-            return "$className(shape=${structure.shape.contentToString()}, buffer=$bufferRepr)"
+            return "$className(shape=${structure.shape}, buffer=$bufferRepr)"
         }
 
         /**
@@ -145,28 +147,28 @@ public interface StructureND<out T> : Featured<StructureFeature>, WithShape {
         ): BufferND<T> = BufferND(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) })
 
         public fun <T> buffered(
-            shape: IntArray,
+            shape: ShapeND,
             bufferFactory: BufferFactory<T> = BufferFactory.boxing(),
             initializer: (IntArray) -> T,
-        ): BufferND<T> = buffered(DefaultStrides(shape), bufferFactory, initializer)
+        ): BufferND<T> = buffered(ColumnStrides(shape), bufferFactory, initializer)
 
         public inline fun <reified T : Any> auto(
-            shape: IntArray,
+            shape: ShapeND,
             crossinline initializer: (IntArray) -> T,
-        ): BufferND<T> = auto(DefaultStrides(shape), initializer)
+        ): BufferND<T> = auto(ColumnStrides(shape), initializer)
 
         @JvmName("autoVarArg")
         public inline fun <reified T : Any> auto(
             vararg shape: Int,
             crossinline initializer: (IntArray) -> T,
         ): BufferND<T> =
-            auto(DefaultStrides(shape), initializer)
+            auto(ColumnStrides(ShapeND(shape)), initializer)
 
         public inline fun <T : Any> auto(
             type: KClass<T>,
             vararg shape: Int,
             crossinline initializer: (IntArray) -> T,
-        ): BufferND<T> = auto(type, DefaultStrides(shape), initializer)
+        ): BufferND<T> = auto(type, ColumnStrides(ShapeND(shape)), initializer)
     }
 }
 
@@ -214,8 +216,13 @@ public fun <T : Comparable<T>> LinearSpace<T, Ring<T>>.contentEquals(
  * @param index the indices.
  * @return the value.
  */
+@PerformancePitfall
 public operator fun <T> StructureND<T>.get(vararg index: Int): T = get(index)
 
+public operator fun StructureND<Double>.get(vararg index: Int): Double = getDouble(index)
+
+public operator fun StructureND<Int>.get(vararg index: Int): Int = getInt(index)
+
 //@UnstableKMathAPI
 //public inline fun <reified T : StructureFeature> StructureND<*>.getFeature(): T? = getFeature(T::class)
 
@@ -229,12 +236,14 @@ public interface MutableStructureND<T> : StructureND<T> {
      * @param index the indices.
      * @param value the value.
      */
+    @PerformancePitfall
     public operator fun set(index: IntArray, value: T)
 }
 
 /**
  * Set value at specified indices
  */
+@PerformancePitfall
 public operator fun <T> MutableStructureND<T>.set(vararg index: Int, value: T) {
     set(index, value)
 }

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

@@ -5,12 +5,15 @@
 
 package space.kscience.kmath.nd
 
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 
 public open class VirtualStructureND<T>(
-    override val shape: Shape,
+    override val shape: ShapeND,
     public val producer: (IntArray) -> T,
 ) : StructureND<T> {
+
+    @PerformancePitfall
     override fun get(index: IntArray): T {
         requireIndexInShape(index, shape)
         return producer(index)
@@ -19,12 +22,12 @@ public open class VirtualStructureND<T>(
 
 @UnstableKMathAPI
 public class VirtualDoubleStructureND(
-    shape: Shape,
+    shape: ShapeND,
     producer: (IntArray) -> Double,
 ) : VirtualStructureND<Double>(shape, producer)
 
 @UnstableKMathAPI
 public class VirtualIntStructureND(
-    shape: Shape,
+    shape: ShapeND,
     producer: (IntArray) -> Int,
 ) : VirtualStructureND<Int>(shape, producer)

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

@@ -15,9 +15,9 @@ public fun <T, A : Algebra<T>> AlgebraND<T, A>.structureND(
     shapeFirst: Int,
     vararg shapeRest: Int,
     initializer: A.(IntArray) -> T
-): StructureND<T> = structureND(Shape(shapeFirst, *shapeRest), initializer)
+): StructureND<T> = structureND(ShapeND(shapeFirst, *shapeRest), initializer)
 
-public fun <T, A : Group<T>> AlgebraND<T, A>.zero(shape: Shape): StructureND<T> = structureND(shape) { zero }
+public fun <T, A : Group<T>> AlgebraND<T, A>.zero(shape: ShapeND): StructureND<T> = structureND(shape) { zero }
 
 @JvmName("zeroVarArg")
 public fun <T, A : Group<T>> AlgebraND<T, A>.zero(
@@ -25,7 +25,7 @@ public fun <T, A : Group<T>> AlgebraND<T, A>.zero(
     vararg shapeRest: Int,
 ): StructureND<T> = structureND(shapeFirst, *shapeRest) { zero }
 
-public fun <T, A : Ring<T>> AlgebraND<T, A>.one(shape: Shape): StructureND<T> = structureND(shape) { one }
+public fun <T, A : Ring<T>> AlgebraND<T, A>.one(shape: ShapeND): StructureND<T> = structureND(shape) { one }
 
 @JvmName("oneVarArg")
 public fun <T, A : Ring<T>> AlgebraND<T, A>.one(

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

@@ -5,6 +5,9 @@
 
 package space.kscience.kmath.nd
 
+import space.kscience.kmath.misc.PerformancePitfall
+
+@OptIn(PerformancePitfall::class)
 public fun <T> StructureND<T>.roll(axis: Int, step: Int = 1): StructureND<T> {
     require(axis in shape.indices) { "Axis $axis is outside of shape dimensions: [0, ${shape.size})" }
     return VirtualStructureND(shape) { index ->
@@ -19,6 +22,7 @@ public fun <T> StructureND<T>.roll(axis: Int, step: Int = 1): StructureND<T> {
     }
 }
 
+@OptIn(PerformancePitfall::class)
 public fun <T> StructureND<T>.roll(pair: Pair<Int, Int>, vararg others: Pair<Int, Int>): StructureND<T> {
     val axisMap: Map<Int, Int> = mapOf(pair, *others)
     require(axisMap.keys.all { it in shape.indices }) { "Some of axes ${axisMap.keys} is outside of shape dimensions: [0, ${shape.size})" }

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

@@ -0,0 +1,45 @@
+/*
+ * Copyright 2018-2022 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.nd
+
+import space.kscience.kmath.misc.PerformancePitfall
+
+public interface StructureNDOfDouble : StructureND<Double> {
+    /**
+     * Guaranteed non-blocking access to content
+     */
+    public fun getDouble(index: IntArray): Double
+}
+
+/**
+ * Optimized method to access primitive without boxing if possible
+ */
+@OptIn(PerformancePitfall::class)
+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
+     */
+    public fun setDouble(index: IntArray, value: Double)
+}
+
+@OptIn(PerformancePitfall::class)
+public fun MutableStructureND<Double>.getDouble(index: IntArray): Double =
+    if (this is StructureNDOfDouble) getDouble(index) else get(index)
+
+
+public interface StructureNDOfInt : StructureND<Int> {
+    /**
+     * Guaranteed non-blocking access to content
+     */
+    public fun getInt(index: IntArray): Int
+}
+
+@OptIn(PerformancePitfall::class)
+public fun StructureND<Int>.getInt(index: IntArray): Int =
+    if (this is StructureNDOfInt) getInt(index) else get(index)

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

@@ -5,8 +5,6 @@
 
 package space.kscience.kmath.operations
 
-import space.kscience.kmath.misc.UnstableKMathAPI
-import space.kscience.kmath.operations.DoubleField.pow
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
 
@@ -32,9 +30,9 @@ public class DoubleBufferField(public val size: Int) : ExtendedField<Buffer<Doub
     override fun atanh(arg: Buffer<Double>): DoubleBuffer = super<DoubleBufferOps>.atanh(arg)
 
     override fun power(arg: Buffer<Double>, pow: Number): DoubleBuffer = if (pow.isInteger()) {
-        arg.mapInline { it.pow(pow.toInt()) }
+        arg.map { it.pow(pow.toInt()) }
     } else {
-        arg.mapInline {
+        arg.map {
             if(it<0) throw IllegalArgumentException("Negative argument $it could not be raised to the fractional power")
             it.pow(pow.toDouble())
         }
@@ -42,103 +40,4 @@ public class DoubleBufferField(public val size: Int) : ExtendedField<Buffer<Doub
 
     override fun unaryOperationFunction(operation: String): (arg: Buffer<Double>) -> Buffer<Double> =
         super<ExtendedField>.unaryOperationFunction(operation)
-
-    //    override fun number(value: Number): Buffer<Double> = DoubleBuffer(size) { value.toDouble() }
-//
-//    override fun Buffer<Double>.unaryMinus(): Buffer<Double> = DoubleBufferOperations.run {
-//        -this@unaryMinus
-//    }
-//
-//    override fun add(a: Buffer<Double>, b: Buffer<Double>): DoubleBuffer {
-//        require(a.size == size) { "The buffer size ${a.size} does not match context size $size" }
-//        return DoubleBufferOperations.add(a, b)
-//    }
-//
-
-//
-//    override fun multiply(a: Buffer<Double>, b: Buffer<Double>): DoubleBuffer {
-//        require(a.size == size) { "The buffer size ${a.size} does not match context size $size" }
-//        return DoubleBufferOperations.multiply(a, b)
-//    }
-//
-//    override fun divide(a: Buffer<Double>, b: Buffer<Double>): DoubleBuffer {
-//        require(a.size == size) { "The buffer size ${a.size} does not match context size $size" }
-//        return DoubleBufferOperations.divide(a, b)
-//    }
-//
-//    override fun sin(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.sin(arg)
-//    }
-//
-//    override fun cos(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.cos(arg)
-//    }
-//
-//    override fun tan(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.tan(arg)
-//    }
-//
-//    override fun asin(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.asin(arg)
-//    }
-//
-//    override fun acos(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.acos(arg)
-//    }
-//
-//    override fun atan(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.atan(arg)
-//    }
-//
-//    override fun sinh(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.sinh(arg)
-//    }
-//
-//    override fun cosh(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.cosh(arg)
-//    }
-//
-//    override fun tanh(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.tanh(arg)
-//    }
-//
-//    override fun asinh(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.asinh(arg)
-//    }
-//
-//    override fun acosh(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.acosh(arg)
-//    }
-//
-//    override fun atanh(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.atanh(arg)
-//    }
-//
-//    override fun power(arg: Buffer<Double>, pow: Number): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.power(arg, pow)
-//    }
-//
-//    override fun exp(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.exp(arg)
-//    }
-//
-//    override fun ln(arg: Buffer<Double>): DoubleBuffer {
-//        require(arg.size == size) { "The buffer size ${arg.size} does not match context size $size" }
-//        return DoubleBufferOperations.ln(arg)
-//    }
-
 }

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

@@ -6,10 +6,8 @@
 package space.kscience.kmath.operations
 
 import space.kscience.kmath.linear.Point
-import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.MutableBufferFactory
-import space.kscience.kmath.structures.asBuffer
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.structures.*
 import kotlin.math.*
 
 /**
@@ -19,10 +17,29 @@ public abstract class DoubleBufferOps : BufferAlgebra<Double, DoubleField>, Exte
     Norm<Buffer<Double>, Double> {
 
     override val elementAlgebra: DoubleField get() = DoubleField
+
     override val elementBufferFactory: MutableBufferFactory<Double> get() = elementAlgebra.bufferFactory
 
-    override fun Buffer<Double>.map(block: DoubleField.(Double) -> Double): DoubleBuffer =
-        mapInline { DoubleField.block(it) }
+    @Suppress("OVERRIDE_BY_INLINE")
+    @OptIn(UnstableKMathAPI::class)
+    final override inline fun Buffer<Double>.map(block: DoubleField.(Double) -> Double): DoubleBuffer =
+        DoubleArray(size) { DoubleField.block(getDouble(it)) }.asBuffer()
+
+
+    @OptIn(UnstableKMathAPI::class)
+    @Suppress("OVERRIDE_BY_INLINE")
+    final override inline fun Buffer<Double>.mapIndexed(block: DoubleField.(index: Int, arg: Double) -> Double): DoubleBuffer =
+        DoubleBuffer(size) { DoubleField.block(it, getDouble(it)) }
+
+    @OptIn(UnstableKMathAPI::class)
+    @Suppress("OVERRIDE_BY_INLINE")
+    final override inline fun Buffer<Double>.zip(
+        other: Buffer<Double>,
+        block: DoubleField.(left: Double, right: Double) -> Double,
+    ): DoubleBuffer {
+        require(size == other.size) { "Incompatible buffer sizes. left: ${size}, right: ${other.size}" }
+        return DoubleBuffer(size) { DoubleField.block(getDouble(it), other.getDouble(it)) }
+    }
 
     override fun unaryOperationFunction(operation: String): (arg: Buffer<Double>) -> Buffer<Double> =
         super<ExtendedFieldOps>.unaryOperationFunction(operation)
@@ -30,7 +47,7 @@ public abstract class DoubleBufferOps : BufferAlgebra<Double, DoubleField>, Exte
     override fun binaryOperationFunction(operation: String): (left: Buffer<Double>, right: Buffer<Double>) -> Buffer<Double> =
         super<ExtendedFieldOps>.binaryOperationFunction(operation)
 
-    override fun Buffer<Double>.unaryMinus(): DoubleBuffer = mapInline { -it }
+    override fun Buffer<Double>.unaryMinus(): DoubleBuffer = map { -it }
 
     override fun add(left: Buffer<Double>, right: Buffer<Double>): DoubleBuffer {
         require(right.size == left.size) {
@@ -77,6 +94,7 @@ public abstract class DoubleBufferOps : BufferAlgebra<Double, DoubleField>, Exte
 //        } else RealBuffer(DoubleArray(a.size) { a[it] / kValue })
 //    }
 
+    @UnstableKMathAPI
     override fun multiply(left: Buffer<Double>, right: Buffer<Double>): DoubleBuffer {
         require(right.size == left.size) {
             "The size of the first buffer ${left.size} should be the same as for second one: ${right.size} "
@@ -101,55 +119,83 @@ public abstract class DoubleBufferOps : BufferAlgebra<Double, DoubleField>, Exte
         } else DoubleBuffer(DoubleArray(left.size) { left[it] / right[it] })
     }
 
-    override fun sin(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::sin)
+    override fun sin(arg: Buffer<Double>): DoubleBuffer = arg.map { sin(it) }
 
-    override fun cos(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::cos)
+    override fun cos(arg: Buffer<Double>): DoubleBuffer = arg.map { cos(it) }
 
-    override fun tan(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::tan)
+    override fun tan(arg: Buffer<Double>): DoubleBuffer = arg.map { tan(it) }
 
-    override fun asin(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::asin)
+    override fun asin(arg: Buffer<Double>): DoubleBuffer = arg.map { asin(it) }
 
-    override fun acos(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::acos)
+    override fun acos(arg: Buffer<Double>): DoubleBuffer = arg.map { acos(it) }
 
-    override fun atan(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::atan)
+    override fun atan(arg: Buffer<Double>): DoubleBuffer = arg.map { atan(it) }
 
-    override fun sinh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::sinh)
+    override fun sinh(arg: Buffer<Double>): DoubleBuffer = arg.map { sinh(it) }
 
-    override fun cosh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::cosh)
+    override fun cosh(arg: Buffer<Double>): DoubleBuffer = arg.map { cosh(it) }
 
-    override fun tanh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::tanh)
+    override fun tanh(arg: Buffer<Double>): DoubleBuffer = arg.map { tanh(it) }
 
-    override fun asinh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::asinh)
+    override fun asinh(arg: Buffer<Double>): DoubleBuffer = arg.map { asinh(it) }
 
-    override fun acosh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::acosh)
+    override fun acosh(arg: Buffer<Double>): DoubleBuffer = arg.map { acosh(it) }
 
-    override fun atanh(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::atanh)
+    override fun atanh(arg: Buffer<Double>): DoubleBuffer = arg.map { atanh(it) }
 
-    override fun exp(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::exp)
+    override fun exp(arg: Buffer<Double>): DoubleBuffer = arg.map { exp(it) }
 
-    override fun ln(arg: Buffer<Double>): DoubleBuffer = arg.mapInline(::ln)
+    override fun ln(arg: Buffer<Double>): DoubleBuffer = arg.map { ln(it) }
 
     override fun norm(arg: Buffer<Double>): Double = DoubleL2Norm.norm(arg)
 
-    override fun scale(a: Buffer<Double>, value: Double): DoubleBuffer = a.mapInline { it * value }
+    override fun scale(a: Buffer<Double>, value: Double): DoubleBuffer = a.map { it * value }
 
     override fun power(arg: Buffer<Double>, pow: Number): Buffer<Double> = if (pow is Int) {
-        arg.mapInline { it.pow(pow) }
+        arg.map { it.pow(pow) }
     } else {
-        arg.mapInline { it.pow(pow.toDouble()) }
+        arg.map { it.pow(pow.toDouble()) }
     }
 
-    public companion object : DoubleBufferOps() {
-        public inline fun Buffer<Double>.mapInline(block: (Double) -> Double): DoubleBuffer =
-            if (this is DoubleBuffer) {
-                DoubleArray(size) { block(array[it]) }.asBuffer()
-            } else {
-                DoubleArray(size) { block(get(it)) }.asBuffer()
-            }
-    }
+    public companion object : DoubleBufferOps()
 }
 
 public object DoubleL2Norm : Norm<Point<Double>, Double> {
     override fun norm(arg: Point<Double>): Double = sqrt(arg.fold(0.0) { acc: Double, d: Double -> acc + d.pow(2) })
 }
 
+public fun DoubleBufferOps.sum(buffer: Buffer<Double>): Double = buffer.reduce(Double::plus)
+
+/**
+ * Sum of elements using given [conversion]
+ */
+public inline fun <T> DoubleBufferOps.sumOf(buffer: Buffer<T>, conversion: (T) -> Double): Double =
+    buffer.fold(0.0) { acc, value -> acc + conversion(value) }
+
+public fun DoubleBufferOps.average(buffer: Buffer<Double>): Double = sum(buffer) / buffer.size
+
+/**
+ * Average of elements using given [conversion]
+ */
+public inline fun <T> DoubleBufferOps.averageOf(buffer: Buffer<T>, conversion: (T) -> Double): Double =
+    sumOf(buffer, conversion) / buffer.size
+
+public fun DoubleBufferOps.dispersion(buffer: Buffer<Double>): Double {
+    val av = average(buffer)
+    return buffer.fold(0.0) { acc, value -> acc + (value - av).pow(2) } / buffer.size
+}
+
+public fun DoubleBufferOps.std(buffer: Buffer<Double>): Double = sqrt(dispersion(buffer))
+
+public fun DoubleBufferOps.covariance(x: Buffer<Double>, y: Buffer<Double>): Double {
+    require(x.size == y.size) { "Expected buffers of the same size, but x.size == ${x.size} and y.size == ${y.size}" }
+    val xMean = average(x)
+    val yMean = average(y)
+    var sum = 0.0
+    x.indices.forEach {
+        sum += (x[it] - xMean) * (y[it] - yMean)
+    }
+    return sum / (x.size - 1)
+}
+
+

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

@@ -5,6 +5,21 @@
 
 package space.kscience.kmath.operations
 
+import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.structures.Buffer
+
+/**
+ * Returns the sum of all elements in the iterable in this [Group].
+ *
+ * @receiver the algebra that provides addition.
+ * @param data the iterable to sum up.
+ * @return the sum.
+ */
+@PerformancePitfall("Potential boxing access to buffer elements")
+public fun <T> Group<T>.sum(data: Buffer<T>): T = data.fold(zero) { left, right ->
+    add(left, right)
+}
+
 /**
  * Returns the sum of all elements in the iterable in this [Group].
  *
@@ -29,6 +44,18 @@ public fun <T> Group<T>.sum(data: Sequence<T>): T = data.fold(zero) { left, righ
     add(left, right)
 }
 
+/**
+ * Returns an average value of elements in the iterable in this [Group].
+ *
+ * @receiver the algebra that provides addition and division.
+ * @param data the iterable to find average.
+ * @return the average value.
+ * @author Iaroslav Postovalov
+ */
+@PerformancePitfall("Potential boxing access to buffer elements")
+public fun <T, S> S.average(data: Buffer<T>): T where S : Group<T>, S : ScaleOperations<T> =
+    sum(data) / data.size
+
 /**
  * Returns an average value of elements in the iterable in this [Group].
  *
@@ -95,4 +122,3 @@ public fun <T, S> Iterable<T>.averageWith(space: S): T where S : Group<T>, S : S
  */
 public fun <T, S> Sequence<T>.averageWith(space: S): T where S : Group<T>, S : ScaleOperations<T> =
     space.average(this)
-

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

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.structures
 
+import space.kscience.kmath.operations.WithSize
 import space.kscience.kmath.operations.asSequence
 import kotlin.jvm.JvmInline
 import kotlin.reflect.KClass
@@ -50,11 +51,11 @@ public fun interface MutableBufferFactory<T> : BufferFactory<T> {
  *
  * @param T the type of elements contained in the buffer.
  */
-public interface Buffer<out T> {
+public interface Buffer<out T> : WithSize {
     /**
      * The size of this buffer.
      */
-    public val size: Int
+    override val size: Int
 
     /**
      * Gets element at given index.
@@ -64,7 +65,7 @@ public interface Buffer<out T> {
     /**
      * Iterates over all elements.
      */
-    public operator fun iterator(): Iterator<T>
+    public operator fun iterator(): Iterator<T> = indices.asSequence().map(::get).iterator()
 
     override fun toString(): String
 
@@ -122,7 +123,14 @@ public interface Buffer<out T> {
 /**
  * Returns an [IntRange] of the valid indices for this [Buffer].
  */
-public val Buffer<*>.indices: IntRange get() = 0 until size
+public val <T> Buffer<T>.indices: IntRange get() = 0 until size
+
+public operator fun <T> Buffer<T>.get(index: UInt): T = get(index.toInt())
+
+/**
+ * if index is in range of buffer, return the value. Otherwise, return null.
+ */
+public fun <T> Buffer<T>.getOrNull(index: Int): T? = if (index in indices) get(index) else null
 
 public fun <T> Buffer<T>.first(): T {
     require(size > 0) { "Can't get the first element of empty buffer" }

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

@@ -5,10 +5,7 @@
 
 package space.kscience.kmath.structures
 
-import space.kscience.kmath.nd.DefaultStrides
-import space.kscience.kmath.nd.Structure2D
-import space.kscience.kmath.nd.StructureND
-import space.kscience.kmath.nd.as2D
+import space.kscience.kmath.nd.*
 
 /**
  * A context that allows to operate on a [MutableBuffer] as on 2d array
@@ -31,7 +28,7 @@ internal class BufferAccessor2D<T>(
 
     //TODO optimize wrapper
     fun MutableBuffer<T>.collect(): Structure2D<T> = StructureND.buffered(
-        DefaultStrides(intArrayOf(rowNum, colNum)),
+        ColumnStrides(ShapeND(rowNum, colNum)),
         factory
     ) { (i, j) ->
         get(i, j)

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

@@ -0,0 +1,192 @@
+package space.kscience.kmath.structures
+
+import space.kscience.kmath.misc.UnstableKMathAPI
+
+/**
+ * A buffer that wraps an original buffer
+ */
+public interface BufferView<T> : Buffer<T> {
+    public val origin: Buffer<T>
+
+    /**
+     * Get the index in [origin] buffer from index in this buffer.
+     * Return -1 if element not present in the original buffer
+     * This method should be used internally to optimize non-boxing access.
+     */
+    @UnstableKMathAPI
+    public fun originIndex(index: Int): Int
+
+    @OptIn(UnstableKMathAPI::class)
+    override fun get(index: Int): T = origin[originIndex(index)]
+}
+
+/**
+ * A zero-copy buffer that "sees" only part of original buffer. Slice can't go beyond original buffer borders.
+ */
+public class BufferSlice<T>(
+    override val origin: Buffer<T>,
+    public val offset: Int = 0,
+    override val size: Int,
+) : BufferView<T> {
+
+    init {
+        require(size > 0) { "Size must be positive" }
+        require(offset + size <= origin.size) {
+            "End of buffer ${offset + size} is beyond the end of origin buffer size ${origin.size}"
+        }
+    }
+
+    override fun get(index: Int): T = if (index >= size) {
+        throw IndexOutOfBoundsException("$index is out of ${0 until size} rage")
+    } else {
+        origin[index + offset]
+    }
+
+    override fun iterator(): Iterator<T> =
+        (offset until (offset + size)).asSequence().map { origin[it] }.iterator()
+
+    @UnstableKMathAPI
+    override fun originIndex(index: Int): Int = if (index >= size) -1 else index - offset
+
+    override fun toString(): String = "$origin[$offset..${offset + size}"
+}
+
+/**
+ * An expanded buffer that could include the whole initial buffer or its part and fills all space beyond it borders with [defaultValue].
+ *
+ * The [offset] parameter shows the shift of expanded buffer start relative to origin start and could be both positive and negative.
+ */
+public class BufferExpanded<T>(
+    override val origin: Buffer<T>,
+    private val defaultValue: T,
+    public val offset: Int = 0,
+    override val size: Int = origin.size,
+) : BufferView<T> {
+
+    init {
+        require(size > 0) { "Size must be positive" }
+    }
+
+    override fun get(index: Int): T = when (index) {
+        !in 0 until size -> throw IndexOutOfBoundsException("Index $index is not in $indices")
+        in -offset until origin.size - offset -> origin[index + offset]
+        else -> defaultValue
+    }
+
+    @UnstableKMathAPI
+    override fun originIndex(index: Int): Int = if (index in -offset until origin.size - offset) index + offset else -1
+
+    override fun toString(): String = "$origin[$offset..${offset + size}]"
+}
+
+/**
+ * Zero-copy select a slice inside the original buffer
+ */
+public fun <T> Buffer<T>.slice(range: IntRange): BufferView<T> = if (this is BufferSlice) {
+    BufferSlice(
+        origin,
+        this.offset + range.first,
+        (range.last - range.first) + 1
+    )
+} else {
+    BufferSlice(
+        this,
+        range.first,
+        (range.last - range.first) + 1
+    )
+}
+
+/**
+ * Resize original buffer to a given range using given [range], filling additional segments with [defaultValue].
+ * Range left border could be negative to designate adding new blank segment to the beginning of the buffer
+ */
+public fun <T> Buffer<T>.expand(
+    range: IntRange,
+    defaultValue: T,
+): BufferView<T> = if (range.first >= 0 && range.last < size) {
+    BufferSlice(
+        this,
+        range.first,
+        (range.last - range.first) + 1
+    )
+} else {
+    BufferExpanded(
+        this,
+        defaultValue,
+        range.first,
+        (range.last - range.first) + 1
+    )
+}
+
+/**
+ * A [BufferView] that overrides indexing of the original buffer
+ */
+public class PermutedBuffer<T>(
+    override val origin: Buffer<T>,
+    private val permutations: IntArray,
+) : BufferView<T> {
+    init {
+        permutations.forEach { index ->
+            if (index !in origin.indices) {
+                throw IndexOutOfBoundsException("Index $index is not in ${origin.indices}")
+            }
+        }
+    }
+
+    override val size: Int get() = permutations.size
+
+    override fun get(index: Int): T = origin[permutations[index]]
+
+    override fun iterator(): Iterator<T> = permutations.asSequence().map { origin[it] }.iterator()
+
+    @UnstableKMathAPI
+    override fun originIndex(index: Int): Int = if (index in permutations.indices) permutations[index] else -1
+
+    override fun toString(): String = Buffer.toString(this)
+}
+
+/**
+ * Created a permuted view of given buffer using provided [indices]
+ */
+public fun <T> Buffer<T>.permute(indices: IntArray): PermutedBuffer<T> =
+    PermutedBuffer(this, indices)
+
+/**
+ * A [BufferView] that overrides indexing of the original buffer
+ */
+public class PermutedMutableBuffer<T>(
+    override val origin: MutableBuffer<T>,
+    private val permutations: IntArray,
+) : BufferView<T>, MutableBuffer<T> {
+    init {
+        permutations.forEach { index ->
+            if (index !in origin.indices) {
+                throw IndexOutOfBoundsException("Index $index is not in ${origin.indices}")
+            }
+        }
+    }
+
+    override val size: Int get() = permutations.size
+
+    override fun get(index: Int): T = origin[permutations[index]]
+
+    override fun set(index: Int, value: T) {
+        origin[permutations[index]] = value
+    }
+
+    override fun copy(): MutableBuffer<T> = PermutedMutableBuffer(origin.copy(), permutations)
+    //TODO Probably could be optimized
+
+    override fun iterator(): Iterator<T> = permutations.asSequence().map { origin[it] }.iterator()
+
+    @UnstableKMathAPI
+    override fun originIndex(index: Int): Int = if (index in permutations.indices) permutations[index] else -1
+
+    override fun toString(): String = Buffer.toString(this)
+}
+
+/**
+ * Created a permuted mutable view of given buffer using provided [indices]
+ */
+public fun <T> MutableBuffer<T>.permute(indices: IntArray): PermutedMutableBuffer<T> =
+    PermutedMutableBuffer(this, indices)

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

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.structures
 
+import space.kscience.kmath.operations.BufferTransform
 import kotlin.jvm.JvmInline
 
 /**
@@ -13,7 +14,7 @@ import kotlin.jvm.JvmInline
  * @property array the underlying array.
  */
 @JvmInline
-public value class DoubleBuffer(public val array: DoubleArray) : MutableBuffer<Double> {
+public value class DoubleBuffer(public val array: DoubleArray) : PrimitiveBuffer<Double> {
     override val size: Int get() = array.size
 
     override operator fun get(index: Int): Double = array[index]
@@ -28,7 +29,7 @@ public value class DoubleBuffer(public val array: DoubleArray) : MutableBuffer<D
 
     override fun toString(): String = Buffer.toString(this)
 
-    public companion object{
+    public companion object {
         public fun zero(size: Int): DoubleBuffer = DoubleArray(size).asBuffer()
     }
 }
@@ -40,7 +41,8 @@ public value class DoubleBuffer(public val array: DoubleArray) : MutableBuffer<D
  * The function [init] is called for each array element sequentially starting from the first one.
  * It should return the value for a buffer element given its index.
  */
-public inline fun DoubleBuffer(size: Int, init: (Int) -> Double): DoubleBuffer = DoubleBuffer(DoubleArray(size) { init(it) })
+public inline fun DoubleBuffer(size: Int, init: (Int) -> Double): DoubleBuffer =
+    DoubleBuffer(DoubleArray(size) { init(it) })
 
 /**
  * Returns a new [DoubleBuffer] of given elements.
@@ -51,10 +53,18 @@ public fun DoubleBuffer(vararg doubles: Double): DoubleBuffer = DoubleBuffer(dou
  * Returns a new [DoubleArray] containing all the elements of this [Buffer].
  */
 public fun Buffer<Double>.toDoubleArray(): DoubleArray = when (this) {
-    is DoubleBuffer -> array.copyOf()
+    is DoubleBuffer -> array
     else -> DoubleArray(size, ::get)
 }
 
+/**
+ * Represent this buffer as [DoubleBuffer]. Does not guarantee that changes in the original buffer are reflected on this buffer.
+ */
+public fun Buffer<Double>.toDoubleBuffer(): DoubleBuffer = when (this) {
+    is DoubleBuffer -> this
+    else -> DoubleArray(size, ::get).asBuffer()
+}
+
 /**
  * Returns [DoubleBuffer] over this array.
  *
@@ -62,3 +72,10 @@ public fun Buffer<Double>.toDoubleArray(): DoubleArray = when (this) {
  * @return the new buffer.
  */
 public fun DoubleArray.asBuffer(): DoubleBuffer = DoubleBuffer(this)
+
+
+public fun interface DoubleBufferTransform : BufferTransform<Double, Double> {
+    public fun transform(arg: DoubleBuffer): DoubleBuffer
+
+    override fun transform(arg: Buffer<Double>): DoubleBuffer = arg.toDoubleBuffer()
+}

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

@@ -14,7 +14,7 @@ import kotlin.jvm.JvmInline
  * @author Iaroslav Postovalov
  */
 @JvmInline
-public value class FloatBuffer(public val array: FloatArray) : MutableBuffer<Float> {
+public value class FloatBuffer(public val array: FloatArray) : PrimitiveBuffer<Float> {
     override val size: Int get() = array.size
 
     override operator fun get(index: Int): Float = array[index]

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

@@ -13,7 +13,7 @@ import kotlin.jvm.JvmInline
  * @property array the underlying array.
  */
 @JvmInline
-public value class IntBuffer(public val array: IntArray) : MutableBuffer<Int> {
+public value class IntBuffer(public val array: IntArray) : PrimitiveBuffer<Int> {
     override val size: Int get() = array.size
 
     override operator fun get(index: Int): Int = array[index]

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

@@ -13,7 +13,7 @@ import kotlin.jvm.JvmInline
  * @property array the underlying array.
  */
 @JvmInline
-public value class LongBuffer(public val array: LongArray) : MutableBuffer<Long> {
+public value class LongBuffer(public val array: LongArray) : PrimitiveBuffer<Long> {
     override val size: Int get() = array.size
 
     override operator fun get(index: Int): Long = array[index]

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

@@ -94,4 +94,7 @@ public interface MutableBuffer<T> : Buffer<T> {
         public inline fun <reified T : Any> auto(size: Int, initializer: (Int) -> T): MutableBuffer<T> =
             auto(T::class, size, initializer)
     }
-}
+}
+
+
+public sealed interface PrimitiveBuffer<T>: MutableBuffer<T>

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

@@ -9,14 +9,18 @@ import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.structures.*
 
 /**
- * Typealias for buffer transformations.
+ * Type alias for buffer transformations.
  */
-public typealias BufferTransform<T, R> = (Buffer<T>) -> Buffer<R>
+public fun interface BufferTransform<T, R> {
+    public fun transform(arg: Buffer<T>): Buffer<R>
+}
 
-/**
- * Typealias for buffer transformations with suspend function.
- */
-public typealias SuspendBufferTransform<T, R> = suspend (Buffer<T>) -> Buffer<R>
+///**
+// * Type alias for buffer transformations with suspend function.
+// */
+//public fun interface SuspendBufferTransform<T, R>{
+//    public suspend fun transform(arg: Buffer<T>): Buffer<R>
+//}
 
 
 /**
@@ -57,18 +61,18 @@ public fun <T> Buffer<T>.toMutableList(): MutableList<T> = when (this) {
  */
 @UnstableKMathAPI
 public inline fun <reified T> Buffer<T>.toTypedArray(): Array<T> = Array(size, ::get)
-
-/**
- * Create a new buffer from this one with the given mapping function and using [Buffer.Companion.auto] buffer factory.
- */
-public inline fun <T, reified R : Any> Buffer<T>.map(block: (T) -> R): Buffer<R> =
-    Buffer.auto(size) { block(get(it)) }
+//
+///**
+// * Create a new buffer from this one with the given mapping function and using [Buffer.Companion.auto] buffer factory.
+// */
+//public inline fun <T, reified R : Any> Buffer<T>.map(block: (T) -> R): Buffer<R> =
+//    Buffer.auto(size) { block(get(it)) }
 
 /**
  * Create a new buffer from this one with the given mapping function.
  * Provided [bufferFactory] is used to construct the new buffer.
  */
-public inline fun <T, R> Buffer<T>.map(
+public inline fun <T, R> Buffer<T>.mapToBuffer(
     bufferFactory: BufferFactory<R>,
     crossinline block: (T) -> R,
 ): Buffer<R> = bufferFactory(size) { block(get(it)) }
@@ -77,23 +81,24 @@ public inline fun <T, R> Buffer<T>.map(
  * Create a new buffer from this one with the given mapping (indexed) function.
  * Provided [bufferFactory] is used to construct the new buffer.
  */
-public inline fun <T, R> Buffer<T>.mapIndexed(
+public inline fun <T, R> Buffer<T>.mapIndexedToBuffer(
     bufferFactory: BufferFactory<R>,
     crossinline block: (index: Int, value: T) -> R,
 ): Buffer<R> = bufferFactory(size) { block(it, get(it)) }
-
-/**
- * Create a new buffer from this one with the given indexed mapping function.
- * Provided [BufferFactory] is used to construct the new buffer.
- */
-public inline fun <T, reified R : Any> Buffer<T>.mapIndexed(
-    crossinline block: (index: Int, value: T) -> R,
-): Buffer<R> = Buffer.auto(size) { block(it, get(it)) }
+//
+///**
+// * Create a new buffer from this one with the given indexed mapping function.
+// * Provided [BufferFactory] is used to construct the new buffer.
+// */
+//public inline fun <T, reified R : Any> Buffer<T>.mapIndexed(
+//    crossinline block: (index: Int, value: T) -> R,
+//): Buffer<R> = Buffer.auto(size) { block(it, get(it)) }
 
 /**
  * Fold given buffer according to [operation]
  */
 public inline fun <T, R> Buffer<T>.fold(initial: R, operation: (acc: R, T) -> R): R {
+    if (size == 0) return initial
     var accumulator = initial
     for (index in this.indices) accumulator = operation(accumulator, get(index))
     return accumulator
@@ -103,18 +108,31 @@ public inline fun <T, R> Buffer<T>.fold(initial: R, operation: (acc: R, T) -> R)
  * Fold given buffer according to indexed [operation]
  */
 public inline fun <T : Any, R> Buffer<T>.foldIndexed(initial: R, operation: (index: Int, acc: R, T) -> R): R {
+    if (size == 0) return initial
     var accumulator = initial
     for (index in this.indices) accumulator = operation(index, accumulator, get(index))
     return accumulator
 }
 
+/**
+ * Reduce a buffer from left to right according to [operation]
+ */
+public inline fun <T> Buffer<T>.reduce(operation: (left: T, value: T) -> T): T {
+    require(size > 0) { "Buffer must have elements" }
+    var current = get(0)
+    for (i in 1 until size) {
+        current = operation(current, get(i))
+    }
+    return current
+}
+
 /**
  * Zip two buffers using given [transform].
  */
 @UnstableKMathAPI
-public inline fun <T1, T2 : Any, reified R : Any> Buffer<T1>.zip(
+public inline fun <T1, T2, R> Buffer<T1>.combineToBuffer(
     other: Buffer<T2>,
-    bufferFactory: BufferFactory<R> = BufferFactory.auto(),
+    bufferFactory: BufferFactory<R>,
     crossinline transform: (T1, T2) -> R,
 ): Buffer<R> {
     require(size == other.size) { "Buffer size mismatch in zip: expected $size but found ${other.size}" }

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

@@ -0,0 +1,38 @@
+package space.kscience.kmath.structures
+
+import space.kscience.kmath.misc.UnstableKMathAPI
+
+/**
+ * Non-boxing access to primitive [Double]
+ */
+
+@UnstableKMathAPI
+public fun Buffer<Double>.getDouble(index: Int): Double = if (this is BufferView) {
+    val originIndex = originIndex(index)
+    if (originIndex >= 0) {
+        origin.getDouble(originIndex)
+    } else {
+        get(index)
+    }
+} else if (this is DoubleBuffer) {
+    array[index]
+} else {
+    get(index)
+}
+
+/**
+ * Non-boxing access to primitive [Int]
+ */
+@UnstableKMathAPI
+public fun Buffer<Int>.getInt(index: Int): Int = if (this is BufferView) {
+    val originIndex = originIndex(index)
+    if (originIndex >= 0) {
+        origin.getInt(originIndex)
+    } else {
+        get(index)
+    }
+} else if (this is IntBuffer) {
+    array[index]
+} else {
+    get(index)
+}

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

@@ -0,0 +1,20 @@
+/*
+ * Copyright 2018-2022 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
+
+
+public typealias Float32 = Float
+public typealias Float64 = Double
+
+public typealias Int8 = Byte
+public typealias Int16 = Short
+public typealias Int32 = Int
+public typealias Int64 = Long
+
+public typealias UInt8 = UByte
+public typealias UInt16 = UShort
+public typealias UInt32 = UInt
+public typealias UInt64 = ULong

kmath-core/src/commonTest/kotlin/space/kscience/kmath/nd/StridesTest.kt

@@ -0,0 +1,38 @@
+/*
+ * Copyright 2018-2022 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.nd
+
+import kotlin.test.Test
+
+class StridesTest {
+    @Test
+    fun checkRowBasedStrides() {
+        val strides = RowStrides(ShapeND(3, 3))
+        var counter = 0
+        for(i in 0..2){
+            for(j in 0..2){
+//                print(strides.offset(intArrayOf(i,j)).toString() + "\t")
+                require(strides.offset(intArrayOf(i,j)) == counter)
+                counter++
+            }
+            println()
+        }
+    }
+
+    @Test
+    fun checkColumnBasedStrides() {
+        val strides = ColumnStrides(ShapeND(3, 3))
+        var counter = 0
+        for(i in 0..2){
+            for(j in 0..2){
+//                print(strides.offset(intArrayOf(i,j)).toString() + "\t")
+                require(strides.offset(intArrayOf(j,i)) == counter)
+                counter++
+            }
+            println()
+        }
+    }
+}

kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/BufferExpandedTest.kt

@@ -0,0 +1,27 @@
+package space.kscience.kmath.structures
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertFails
+
+internal class BufferExpandedTest {
+    private val buffer = (0..100).toList().asBuffer()
+
+    @Test
+    fun shrink(){
+        val view = buffer.slice(20..30)
+        assertEquals(20, view[0])
+        assertEquals(30, view[10])
+        assertFails { view[11] }
+    }
+
+    @Test
+    fun expandNegative(){
+        val view: BufferView<Int> = buffer.expand(-20..113,0)
+        assertEquals(0,view[4])
+        assertEquals(0,view[123])
+        assertEquals(100, view[120])
+        assertFails { view[-2] }
+        assertFails { view[134] }
+    }
+}

kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt

@@ -88,7 +88,7 @@ class NumberNDFieldTest {
     @Test
     fun testInternalContext() {
         algebra {
-            (DoubleField.ndAlgebra(*array1.shape)) { with(L2Norm) { 1 + norm(array1) + exp(array2) } }
+            (DoubleField.ndAlgebra(array1.shape)) { with(L2Norm) { 1 + norm(array1) + exp(array2) } }
         }
     }
 }

kmath-core/src/jvmTest/kotlin/space/kscience/kmath/misc/JBigTest.kt

@@ -0,0 +1,20 @@
+/*
+ * Copyright 2018-2022 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.misc
+
+import org.junit.jupiter.api.Test
+import space.kscience.kmath.operations.JBigDecimalField
+import kotlin.test.assertEquals
+import kotlin.test.assertNotEquals
+
+class JBigTest {
+
+    @Test
+    fun testExact() = with(JBigDecimalField) {
+        assertNotEquals(0.3, 0.1 + 0.2)
+        assertEquals(one * 0.3, one * 0.1 + one * 0.2)
+    }
+}

kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/coroutines/coroutinesExtra.kt

@@ -81,9 +81,7 @@ public suspend fun <T> AsyncFlow<T>.collect(concurrency: Int, collector: FlowCol
 public suspend inline fun <T> AsyncFlow<T>.collect(
     concurrency: Int,
     crossinline action: suspend (value: T) -> Unit,
-): Unit = collect(concurrency, object : FlowCollector<T> {
-    override suspend fun emit(value: T): Unit = action(value)
-})
+): Unit = collect(concurrency, FlowCollector<T> { value -> action(value) })
 
 public inline fun <T, R> Flow<T>.mapParallel(
     dispatcher: CoroutineDispatcher = Dispatchers.Default,

kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/BufferFlow.kt

@@ -5,8 +5,10 @@
 
 package space.kscience.kmath.streaming
 
-import kotlinx.coroutines.FlowPreview
-import kotlinx.coroutines.flow.*
+import kotlinx.coroutines.flow.Flow
+import kotlinx.coroutines.flow.asFlow
+import kotlinx.coroutines.flow.flatMapConcat
+import kotlinx.coroutines.flow.flow
 import space.kscience.kmath.chains.BlockingDoubleChain
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.BufferFactory
@@ -20,7 +22,6 @@ public fun <T> Buffer<T>.asFlow(): Flow<T> = iterator().asFlow()
 /**
  * Flat map a [Flow] of [Buffer] into continuous [Flow] of elements
  */
-@FlowPreview
 public fun <T> Flow<Buffer<T>>.spread(): Flow<T> = flatMapConcat { it.asFlow() }
 
 /**

kmath-coroutines/src/jvmMain/kotlin/space/kscience/kmath/structures/LazyStructureND.kt

@@ -8,12 +8,13 @@ package space.kscience.kmath.structures
 import kotlinx.coroutines.*
 import space.kscience.kmath.coroutines.Math
 import space.kscience.kmath.misc.PerformancePitfall
-import space.kscience.kmath.nd.DefaultStrides
+import space.kscience.kmath.nd.ColumnStrides
+import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.nd.StructureND
 
 public class LazyStructureND<out T>(
     public val scope: CoroutineScope,
-    override val shape: IntArray,
+    override val shape: ShapeND,
     public val function: suspend (IntArray) -> T,
 ) : StructureND<T> {
     private val cache: MutableMap<IntArray, Deferred<T>> = HashMap()
@@ -23,30 +24,35 @@ public class LazyStructureND<out T>(
     }
 
     public suspend fun await(index: IntArray): T = async(index).await()
+    @PerformancePitfall
     override operator fun get(index: IntArray): T = runBlocking { async(index).await() }
 
     @OptIn(PerformancePitfall::class)
     override fun elements(): Sequence<Pair<IntArray, T>> {
-        val strides = DefaultStrides(shape)
+        val strides = ColumnStrides(shape)
         val res = runBlocking { strides.asSequence().toList().map { index -> index to await(index) } }
         return res.asSequence()
     }
 }
 
+@OptIn(PerformancePitfall::class)
 public fun <T> StructureND<T>.async(index: IntArray): Deferred<T> =
     if (this is LazyStructureND<T>) this@async.async(index) else CompletableDeferred(get(index))
 
+@OptIn(PerformancePitfall::class)
 public suspend fun <T> StructureND<T>.await(index: IntArray): T =
     if (this is LazyStructureND<T>) await(index) else get(index)
 
 /**
  * PENDING would benefit from KEEP-176
  */
+@OptIn(PerformancePitfall::class)
 public inline fun <T, R> StructureND<T>.mapAsyncIndexed(
     scope: CoroutineScope,
     crossinline function: suspend (T, index: IntArray) -> R,
 ): LazyStructureND<R> = LazyStructureND(scope, shape) { index -> function(get(index), index) }
 
+@OptIn(PerformancePitfall::class)
 public inline fun <T, R> StructureND<T>.mapAsync(
     scope: CoroutineScope,
     crossinline function: suspend (T) -> R,

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

@@ -6,6 +6,7 @@
 package space.kscience.kmath.dimensions
 
 import space.kscience.kmath.linear.*
+import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.nd.Structure2D
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.Ring
@@ -47,7 +48,7 @@ public interface DMatrix<out T, R : Dimension, C : Dimension> : Structure2D<T> {
 public value class DMatrixWrapper<out T, R : Dimension, C : Dimension>(
     private val structure: Structure2D<T>,
 ) : DMatrix<T, R, C> {
-    override val shape: IntArray get() = structure.shape
+    override val shape: ShapeND get() = structure.shape
     override val rowNum: Int get() = shape[0]
     override val colNum: Int get() = shape[1]
     override operator fun get(i: Int, j: Int): T = structure[i, j]

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

@@ -15,6 +15,7 @@ import space.kscience.kmath.linear.*
 import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.nd.toArray
 import space.kscience.kmath.operations.algebra
 import kotlin.random.Random
 import kotlin.random.asJavaRandom
@@ -52,7 +53,7 @@ internal class EjmlMatrixTest {
     fun shape() {
         val m = randomMatrix
         val w = EjmlDoubleMatrix(m)
-        assertContentEquals(intArrayOf(m.numRows, m.numCols), w.shape)
+        assertContentEquals(intArrayOf(m.numRows, m.numCols), w.shape.toArray())
     }
 
     @OptIn(UnstableKMathAPI::class)

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

@@ -131,7 +131,7 @@ public fun RealMatrix.extractColumn(columnIndex: Int): RealMatrix =
     extractColumns(columnIndex..columnIndex)
 
 public fun RealMatrix.sumByColumn(): DoubleBuffer = DoubleBuffer(colNum) { j ->
-    columns[j].asIterable().sum()
+    columns[j].sum()
 }
 
 public fun RealMatrix.minByColumn(): DoubleBuffer = DoubleBuffer(colNum) { j ->

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

@@ -5,7 +5,7 @@
 
 package space.kscience.kmath.integration
 
-import space.kscience.kmath.operations.map
+import space.kscience.kmath.operations.mapToBuffer
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
 import space.kscience.kmath.structures.asBuffer
@@ -30,14 +30,14 @@ public fun GaussIntegratorRuleFactory.build(
     numPoints: Int,
     range: ClosedRange<Double>,
 ): Pair<Buffer<Double>, Buffer<Double>> {
-    val normalized = build(numPoints)
+    val normalized: Pair<Buffer<Double>, Buffer<Double>> = build(numPoints)
     val length = range.endInclusive - range.start
 
-    val points = normalized.first.map(::DoubleBuffer) {
+    val points = normalized.first.mapToBuffer(::DoubleBuffer) {
         range.start + length / 2 + length / 2 * it
     }
 
-    val weights = normalized.second.map(::DoubleBuffer) {
+    val weights = normalized.second.mapToBuffer(::DoubleBuffer) {
         it * length / 2
     }
 

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

@@ -65,9 +65,9 @@ public class SplineIntegrator<T : Comparable<T>>(
             DoubleBuffer(numPoints) { i -> range.start + i * step }
         }
 
-        val values = nodes.map(bufferFactory) { integrand.function(it) }
+        val values = nodes.mapToBuffer(bufferFactory) { integrand.function(it) }
         val polynomials = interpolator.interpolatePolynomials(
-            nodes.map(bufferFactory) { number(it) },
+            nodes.mapToBuffer(bufferFactory) { number(it) },
             values
         )
         val res = polynomials.integrate(algebra, number(range.start)..number(range.endInclusive))
@@ -93,7 +93,7 @@ public object DoubleSplineIntegrator : UnivariateIntegrator<Double> {
             DoubleBuffer(numPoints) { i -> range.start + i * step }
         }
 
-        val values = nodes.map { integrand.function(it) }
+        val values = nodes.mapToBuffer(::DoubleBuffer) { integrand.function(it) }
         val polynomials = interpolator.interpolatePolynomials(nodes, values)
         val res = polynomials.integrate(DoubleField, range)
         return integrand + IntegrandValue(res) + IntegrandCallsPerformed(integrand.calls + nodes.size)

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

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.geometry
 
+
 import space.kscience.kmath.operations.toList
 import kotlin.test.Test
 import kotlin.test.assertEquals

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

@@ -7,9 +7,10 @@ package space.kscience.kmath.histogram
 
 import space.kscience.kmath.domains.Domain
 import space.kscience.kmath.linear.Point
-import space.kscience.kmath.nd.DefaultStrides
+import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.nd.ColumnStrides
 import space.kscience.kmath.nd.FieldOpsND
-import space.kscience.kmath.nd.Shape
+import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.operations.Group
 import space.kscience.kmath.operations.ScaleOperations
@@ -24,6 +25,7 @@ public class HistogramND<T : Comparable<T>, D : Domain<T>, V : Any>(
     internal val values: StructureND<V>,
 ) : Histogram<T, V, DomainBin<T, D, V>> {
 
+    @OptIn(PerformancePitfall::class)
     override fun get(point: Point<T>): DomainBin<T, D, V>? {
         val index = group.getIndexOrNull(point) ?: return null
         return group.produceBin(index, values[index])
@@ -31,8 +33,9 @@ public class HistogramND<T : Comparable<T>, D : Domain<T>, V : Any>(
 
     override val dimension: Int get() = group.shape.size
 
+    @OptIn(PerformancePitfall::class)
     override val bins: Iterable<DomainBin<T, D, V>>
-        get() = DefaultStrides(group.shape).asSequence().map {
+        get() = ColumnStrides(group.shape).asSequence().map {
             group.produceBin(it, values[it])
         }.asIterable()
 }
@@ -42,7 +45,7 @@ public class HistogramND<T : Comparable<T>, D : Domain<T>, V : Any>(
  */
 public interface HistogramGroupND<T : Comparable<T>, D : Domain<T>, V : Any> :
     Group<HistogramND<T, D, V>>, ScaleOperations<HistogramND<T, D, V>> {
-    public val shape: Shape
+    public val shape: ShapeND
     public val valueAlgebraND: FieldOpsND<V, *> //= NDAlgebra.space(valueSpace, Buffer.Companion::boxing, *shape),
 
     /**

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

@@ -9,11 +9,10 @@ package space.kscience.kmath.histogram
 
 import space.kscience.kmath.domains.HyperSquareDomain
 import space.kscience.kmath.linear.Point
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.*
-import space.kscience.kmath.operations.DoubleField
-import space.kscience.kmath.operations.Field
-import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.operations.*
 import space.kscience.kmath.structures.*
 import kotlin.math.floor
 
@@ -40,7 +39,7 @@ public class UniformHistogramGroupND<V : Any, A : Field<V>>(
 
     public val dimension: Int get() = lower.size
 
-    override val shape: IntArray = IntArray(binNums.size) { binNums[it] + 2 }
+    override val shape: ShapeND = ShapeND(IntArray(binNums.size) { binNums[it] + 2 })
 
     private val binSize = DoubleBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] }
 
@@ -83,8 +82,12 @@ public class UniformHistogramGroupND<V : Any, A : Field<V>>(
     }
 
 
-    override fun produce(builder: HistogramBuilder<Double, V>.() -> Unit): HistogramND<Double, HyperSquareDomain, V> {
-        val ndCounter = StructureND.buffered(shape) { Counter.of(valueAlgebraND.elementAlgebra) }
+    @OptIn(PerformancePitfall::class)
+    override fun produce(
+        builder: HistogramBuilder<Double, V>.() -> Unit,
+    ): HistogramND<Double, HyperSquareDomain, V> {
+        val ndCounter: BufferND<ObjectCounter<V>> =
+            StructureND.buffered(shape) { Counter.of(valueAlgebraND.elementAlgebra) }
         val hBuilder = object : HistogramBuilder<Double, V> {
             override val defaultValue: V get() = valueAlgebraND.elementAlgebra.one
 
@@ -94,7 +97,8 @@ public class UniformHistogramGroupND<V : Any, A : Field<V>>(
             }
         }
         hBuilder.apply(builder)
-        val values: BufferND<V> = ndCounter.mapToBuffer(valueBufferFactory) { it.value }
+        val values: BufferND<V> = BufferND(ndCounter.indices, ndCounter.buffer.mapToBuffer(valueBufferFactory) { it.value })
+
         return HistogramND(this, values)
     }
 

kmath-histograms/src/commonTest/kotlin/space/kscience/kmath/histogram/MultivariateHistogramTest.kt

@@ -7,8 +7,9 @@
 
 package space.kscience.kmath.histogram
 
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
-import space.kscience.kmath.nd.DefaultStrides
+import space.kscience.kmath.nd.ColumnStrides
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.real.DoubleVector
 import kotlin.random.Random
@@ -50,6 +51,7 @@ internal class MultivariateHistogramTest {
         assertEquals(n, histogram.bins.sumOf { it.binValue.toInt() })
     }
 
+    @OptIn(PerformancePitfall::class)
     @Test
     fun testHistogramAlgebra() {
         Histogram.uniformDoubleNDFromRanges(
@@ -73,7 +75,7 @@ internal class MultivariateHistogramTest {
             }
             val res = histogram1 - histogram2
             assertTrue {
-                DefaultStrides(shape).asSequence().all { index ->
+                ColumnStrides(shape).asSequence().all { index ->
                     res.values[index] <= histogram1.values[index]
                 }
             }

kmath-histograms/src/commonTest/kotlin/space/kscience/kmath/histogram/UniformHistogram1DTest.kt

@@ -10,7 +10,7 @@ import kotlinx.coroutines.test.runTest
 import space.kscience.kmath.distributions.NormalDistribution
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.DoubleField
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.stat.nextBuffer
 import kotlin.native.concurrent.ThreadLocal
 import kotlin.test.Test
@@ -36,7 +36,7 @@ internal class UniformHistogram1DTest {
     @Test
     fun rebinDown() = runTest {
         val h1 = Histogram.uniform1D(DoubleField, 0.01).produce(generator.nextDoubleBuffer(10000))
-        val h2 = Histogram.uniform1D(DoubleField,0.03).produceFrom(h1)
+        val h2 = Histogram.uniform1D(DoubleField, 0.03).produceFrom(h1)
 
         assertEquals(10000, h2.bins.sumOf { it.binValue }.toInt())
     }
@@ -44,13 +44,13 @@ internal class UniformHistogram1DTest {
     @Test
     fun rebinUp() = runTest {
         val h1 = Histogram.uniform1D(DoubleField, 0.03).produce(generator.nextDoubleBuffer(10000))
-        val h2 = Histogram.uniform1D(DoubleField,0.01).produceFrom(h1)
+        val h2 = Histogram.uniform1D(DoubleField, 0.01).produceFrom(h1)
 
         assertEquals(10000, h2.bins.sumOf { it.binValue }.toInt())
     }
 
     @ThreadLocal
-    companion object{
+    companion object {
         private val generator = RandomGenerator.default(123)
     }
 }

kmath-jupyter/src/main/kotlin/space/kscience/kmath/jupyter/KMathJupyter.kt

@@ -47,11 +47,13 @@ internal class KMathJupyter : JupyterIntegration() {
                 syntaxRender.renderPart(mathRender.render(MST.Numeric(it)), s)
                 +s.toString()
             }
+
             is MST -> {
                 val s = StringBuilder()
                 syntaxRender.renderPart(mathRender.render(it), s)
                 +s.toString()
             }
+
             else -> {
                 +"<ms>"
                 +it.toString()

kmath-multik/src/commonMain/kotlin/space/kscience/kmath/multik/MultikTensor.kt

@@ -7,20 +7,22 @@ package space.kscience.kmath.multik
 
 import org.jetbrains.kotlinx.multik.ndarray.data.*
 import space.kscience.kmath.misc.PerformancePitfall
-import space.kscience.kmath.nd.Shape
+import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.tensors.api.Tensor
 import kotlin.jvm.JvmInline
 
 @JvmInline
 public value class MultikTensor<T>(public val array: MutableMultiArray<T, DN>) : Tensor<T> {
-    override val shape: Shape get() = array.shape
+    override val shape: ShapeND get() = ShapeND(array.shape)
 
+    @PerformancePitfall
     override fun get(index: IntArray): T = array[index]
 
     @PerformancePitfall
     override fun elements(): Sequence<Pair<IntArray, T>> =
         array.multiIndices.iterator().asSequence().map { it to get(it) }
 
+    @PerformancePitfall
     override fun set(index: IntArray, value: T) {
         array[index] = value
     }

kmath-multik/src/commonMain/kotlin/space/kscience/kmath/multik/MultikTensorAlgebra.kt

@@ -14,6 +14,7 @@ import org.jetbrains.kotlinx.multik.api.stat.Statistics
 import org.jetbrains.kotlinx.multik.ndarray.data.*
 import org.jetbrains.kotlinx.multik.ndarray.operations.*
 import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.misc.UnsafeKMathAPI
 import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.*
 import space.kscience.kmath.tensors.api.Tensor
@@ -30,21 +31,22 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>>(
     protected val multikLinAl: LinAlg = multikEngine.getLinAlg()
     protected val multikStat: Statistics = multikEngine.getStatistics()
 
-    override fun structureND(shape: Shape, initializer: A.(IntArray) -> T): MultikTensor<T> {
-        val strides = DefaultStrides(shape)
+    @OptIn(UnsafeKMathAPI::class)
+    override fun structureND(shape: ShapeND, initializer: A.(IntArray) -> T): MultikTensor<T> {
+        val strides = ColumnStrides(shape)
         val memoryView = initMemoryView<T>(strides.linearSize, type)
         strides.asSequence().forEachIndexed { linearIndex, tensorIndex ->
             memoryView[linearIndex] = elementAlgebra.initializer(tensorIndex)
         }
-        return MultikTensor(NDArray(memoryView, shape = shape, dim = DN(shape.size)))
+        return MultikTensor(NDArray(memoryView, shape = shape.asArray(), dim = DN(shape.size)))
     }
 
-    @OptIn(PerformancePitfall::class)
+    @OptIn(PerformancePitfall::class, UnsafeKMathAPI::class)
     override fun StructureND<T>.map(transform: A.(T) -> T): MultikTensor<T> = if (this is MultikTensor) {
         val data = initMemoryView<T>(array.size, type)
         var count = 0
         for (el in array) data[count++] = elementAlgebra.transform(el)
-        NDArray(data, shape = shape, dim = array.dim).wrap()
+        NDArray(data, shape = shape.asArray(), dim = array.dim).wrap()
     } else {
         structureND(shape) { index ->
             transform(get(index))
@@ -75,6 +77,7 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>>(
     /**
      * Transform a structure element-by element in place.
      */
+    @OptIn(PerformancePitfall::class)
     public inline fun <T> MutableStructureND<T>.mapIndexedInPlace(operation: (index: IntArray, t: T) -> T): Unit {
         if (this is MultikTensor) {
             array.multiIndices.iterator().forEach {
@@ -106,10 +109,11 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>>(
      * Convert a tensor to [MultikTensor] if necessary. If tensor is converted, changes on the resulting tensor
      * are not reflected back onto the source
      */
+    @OptIn(UnsafeKMathAPI::class, PerformancePitfall::class)
     public fun StructureND<T>.asMultik(): MultikTensor<T> = if (this is MultikTensor) {
         this
     } else {
-        val res = mk.zeros<T, DN>(shape, type).asDNArray()
+        val res = mk.zeros<T, DN>(shape.asArray(), type).asDNArray()
         for (index in res.multiIndices) {
             res[index] = this[index]
         }
@@ -118,7 +122,8 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>>(
 
     public fun MutableMultiArray<T, *>.wrap(): MultikTensor<T> = MultikTensor(this.asDNArray())
 
-    override fun StructureND<T>.valueOrNull(): T? = if (shape contentEquals intArrayOf(1)) {
+    @OptIn(PerformancePitfall::class)
+    override fun StructureND<T>.valueOrNull(): T? = if (shape contentEquals ShapeND(1)) {
         get(intArrayOf(0))
     } else null
 
@@ -139,6 +144,7 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>>(
         }
     }
 
+    @OptIn(PerformancePitfall::class)
     override fun Tensor<T>.plusAssign(arg: StructureND<T>) {
         if (this is MultikTensor) {
             array.plusAssign(arg.asMultik().array)
@@ -163,6 +169,7 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>>(
         }
     }
 
+    @OptIn(PerformancePitfall::class)
     override fun Tensor<T>.minusAssign(arg: StructureND<T>) {
         if (this is MultikTensor) {
             array.minusAssign(arg.asMultik().array)
@@ -188,6 +195,7 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>>(
         }
     }
 
+    @OptIn(PerformancePitfall::class)
     override fun Tensor<T>.timesAssign(arg: StructureND<T>) {
         if (this is MultikTensor) {
             array.timesAssign(arg.asMultik().array)
@@ -201,13 +209,13 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>>(
 
     override fun Tensor<T>.getTensor(i: Int): MultikTensor<T> = asMultik().array.mutableView(i).wrap()
 
-    override fun Tensor<T>.transposed(i: Int, j: Int): MultikTensor<T> = asMultik().array.transpose(i, j).wrap()
+    override fun StructureND<T>.transposed(i: Int, j: Int): MultikTensor<T> = asMultik().array.transpose(i, j).wrap()
 
-    override fun Tensor<T>.view(shape: IntArray): MultikTensor<T> {
-        require(shape.all { it > 0 })
-        require(shape.fold(1, Int::times) == this.shape.size) {
+    override fun Tensor<T>.view(shape: ShapeND): MultikTensor<T> {
+        require(shape.asList().all { it > 0 })
+        require(shape.linearSize == this.shape.size) {
             "Cannot reshape array of size ${this.shape.size} into a new shape ${
-                shape.joinToString(
+                shape.asList().joinToString(
                     prefix = "(",
                     postfix = ")"
                 )
@@ -215,10 +223,11 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>>(
         }
 
         val mt = asMultik().array
-        return if (mt.shape.contentEquals(shape)) {
+        return if (ShapeND(mt.shape).contentEquals(shape)) {
             mt
         } else {
-            NDArray(mt.data, mt.offset, shape, dim = DN(shape.size), base = mt.base ?: mt)
+            @OptIn(UnsafeKMathAPI::class)
+            NDArray(mt.data, mt.offset, shape.asArray(), dim = DN(shape.size), base = mt.base ?: mt)
         }.wrap()
     }
 
@@ -241,7 +250,7 @@ public abstract class MultikTensorAlgebra<T, A : Ring<T>>(
             TODO("Not implemented for broadcasting")
         }
 
-    override fun diagonalEmbedding(diagonalEntries: Tensor<T>, offset: Int, dim1: Int, dim2: Int): MultikTensor<T> {
+    override fun diagonalEmbedding(diagonalEntries: StructureND<T>, offset: Int, dim1: Int, dim2: Int): MultikTensor<T> {
 
         TODO("Diagonal embedding not implemented")
     }
@@ -284,8 +293,9 @@ public abstract class MultikDivisionTensorAlgebra<T, A : Field<T>>(
     multikEngine: Engine,
 ) : MultikTensorAlgebra<T, A>(multikEngine), TensorPartialDivisionAlgebra<T, A> where T : Number, T : Comparable<T> {
 
+    @OptIn(UnsafeKMathAPI::class)
     override fun T.div(arg: StructureND<T>): MultikTensor<T> =
-        Multik.ones<T, DN>(arg.shape, type).apply { divAssign(arg.asMultik().array) }.wrap()
+        Multik.ones<T, DN>(arg.shape.asArray(), type).apply { divAssign(arg.asMultik().array) }.wrap()
 
     override fun StructureND<T>.div(arg: T): MultikTensor<T> =
         asMultik().array.div(arg).wrap()
@@ -301,6 +311,7 @@ public abstract class MultikDivisionTensorAlgebra<T, A : Field<T>>(
         }
     }
 
+    @OptIn(PerformancePitfall::class)
     override fun Tensor<T>.divAssign(arg: StructureND<T>) {
         if (this is MultikTensor) {
             array.divAssign(arg.asMultik().array)

kmath-multik/src/commonTest/kotlin/space/kscience/kmath/multik/MultikNDTest.kt

@@ -7,10 +7,12 @@ package space.kscience.kmath.multik
 
 import org.jetbrains.kotlinx.multik.default.DefaultEngine
 import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.nd.ShapeND
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.nd.one
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
+import space.kscience.kmath.tensors.core.randomNormal
 import space.kscience.kmath.tensors.core.tensorAlgebra
 import kotlin.test.Test
 import kotlin.test.assertTrue
@@ -28,8 +30,8 @@ internal class MultikNDTest {
     fun dotResult() {
         val dim = 100
 
-        val tensor1 = DoubleTensorAlgebra.randomNormal(shape = intArrayOf(dim, dim), 12224)
-        val tensor2 = DoubleTensorAlgebra.randomNormal(shape = intArrayOf(dim, dim), 12225)
+        val tensor1 = DoubleTensorAlgebra.randomNormal(shape = ShapeND(dim, dim), 12224)
+        val tensor2 = DoubleTensorAlgebra.randomNormal(shape = ShapeND(dim, dim), 12225)
 
         val multikResult = with(multikAlgebra) {
             tensor1 dot tensor2

kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jArrayAlgebra.kt

@@ -11,6 +11,7 @@ import org.nd4j.linalg.api.ops.impl.transforms.strict.ASinh
 import org.nd4j.linalg.factory.Nd4j
 import org.nd4j.linalg.ops.transforms.Transforms
 import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.misc.UnsafeKMathAPI
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.*
@@ -32,8 +33,10 @@ public sealed interface Nd4jArrayAlgebra<T, out C : Algebra<T>> : AlgebraND<T, C
      */
     public val StructureND<T>.ndArray: INDArray
 
-    override fun structureND(shape: Shape, initializer: C.(IntArray) -> T): Nd4jArrayStructure<T> {
-        val struct = Nd4j.create(*shape)!!.wrap()
+    @OptIn(PerformancePitfall::class)
+    override fun structureND(shape: ShapeND, initializer: C.(IntArray) -> T): Nd4jArrayStructure<T> {
+        @OptIn(UnsafeKMathAPI::class)
+        val struct: Nd4jArrayStructure<T> = Nd4j.create(*shape.asArray())!!.wrap()
         struct.indicesIterator().forEach { struct[it] = elementAlgebra.initializer(it) }
         return struct
     }
@@ -45,23 +48,23 @@ public sealed interface Nd4jArrayAlgebra<T, out C : Algebra<T>> : AlgebraND<T, C
         return newStruct
     }
 
-    @OptIn(PerformancePitfall::class)
+    @OptIn(PerformancePitfall::class, UnsafeKMathAPI::class)
     override fun StructureND<T>.mapIndexed(
         transform: C.(index: IntArray, T) -> T,
     ): Nd4jArrayStructure<T> {
-        val new = Nd4j.create(*shape).wrap()
+        val new = Nd4j.create(*shape.asArray()).wrap()
         new.indicesIterator().forEach { idx -> new[idx] = elementAlgebra.transform(idx, this[idx]) }
         return new
     }
 
-    @OptIn(PerformancePitfall::class)
+    @OptIn(PerformancePitfall::class, UnsafeKMathAPI::class)
     override fun zip(
         left: StructureND<T>,
         right: StructureND<T>,
         transform: C.(T, T) -> T,
     ): Nd4jArrayStructure<T> {
         require(left.shape.contentEquals(right.shape)) { "Can't zip tow structures of shape ${left.shape} and ${right.shape}" }
-        val new = Nd4j.create(*left.shape).wrap()
+        val new = Nd4j.create(*left.shape.asArray()).wrap()
         new.indicesIterator().forEach { idx -> new[idx] = elementAlgebra.transform(left[idx], right[idx]) }
         return new
     }
@@ -192,11 +195,11 @@ public open class DoubleNd4jArrayFieldOps : Nd4jArrayExtendedFieldOps<Double, Do
 
     override fun INDArray.wrap(): Nd4jArrayStructure<Double> = asDoubleStructure()
 
-    @OptIn(PerformancePitfall::class)
+    @OptIn(PerformancePitfall::class, UnsafeKMathAPI::class)
     override val StructureND<Double>.ndArray: INDArray
         get() = when (this) {
             is Nd4jArrayStructure<Double> -> ndArray
-            else -> Nd4j.zeros(*shape).also {
+            else -> Nd4j.zeros(*shape.asArray()).also {
                 elements().forEach { (idx, value) -> it.putScalar(idx, value) }
             }
         }
@@ -222,10 +225,10 @@ public open class DoubleNd4jArrayFieldOps : Nd4jArrayExtendedFieldOps<Double, Do
 
 public val DoubleField.nd4j: DoubleNd4jArrayFieldOps get() = DoubleNd4jArrayFieldOps
 
-public class DoubleNd4jArrayField(override val shape: Shape) : DoubleNd4jArrayFieldOps(), FieldND<Double, DoubleField>
+public class DoubleNd4jArrayField(override val shape: ShapeND) : DoubleNd4jArrayFieldOps(), FieldND<Double, DoubleField>
 
 public fun DoubleField.nd4j(shapeFirst: Int, vararg shapeRest: Int): DoubleNd4jArrayField =
-    DoubleNd4jArrayField(intArrayOf(shapeFirst, * shapeRest))
+    DoubleNd4jArrayField(ShapeND(shapeFirst, * shapeRest))
 
 
 /**
@@ -236,11 +239,11 @@ public open class FloatNd4jArrayFieldOps : Nd4jArrayExtendedFieldOps<Float, Floa
 
     override fun INDArray.wrap(): Nd4jArrayStructure<Float> = asFloatStructure()
 
-    @OptIn(PerformancePitfall::class)
+    @OptIn(PerformancePitfall::class, UnsafeKMathAPI::class)
     override val StructureND<Float>.ndArray: INDArray
         get() = when (this) {
             is Nd4jArrayStructure<Float> -> ndArray
-            else -> Nd4j.zeros(*shape).also {
+            else -> Nd4j.zeros(*shape.asArray()).also {
                 elements().forEach { (idx, value) -> it.putScalar(idx, value) }
             }
         }
@@ -269,12 +272,12 @@ public open class FloatNd4jArrayFieldOps : Nd4jArrayExtendedFieldOps<Float, Floa
     public companion object : FloatNd4jArrayFieldOps()
 }
 
-public class FloatNd4jArrayField(override val shape: Shape) : FloatNd4jArrayFieldOps(), RingND<Float, FloatField>
+public class FloatNd4jArrayField(override val shape: ShapeND) : FloatNd4jArrayFieldOps(), RingND<Float, FloatField>
 
 public val FloatField.nd4j: FloatNd4jArrayFieldOps get() = FloatNd4jArrayFieldOps
 
 public fun FloatField.nd4j(shapeFirst: Int, vararg shapeRest: Int): FloatNd4jArrayField =
-    FloatNd4jArrayField(intArrayOf(shapeFirst, * shapeRest))
+    FloatNd4jArrayField(ShapeND(shapeFirst, * shapeRest))
 
 /**
  * Represents [RingND] over [Nd4jArrayIntStructure].
@@ -284,11 +287,11 @@ public open class IntNd4jArrayRingOps : Nd4jArrayRingOps<Int, IntRing> {
 
     override fun INDArray.wrap(): Nd4jArrayStructure<Int> = asIntStructure()
 
-    @OptIn(PerformancePitfall::class)
+    @OptIn(PerformancePitfall::class, UnsafeKMathAPI::class)
     override val StructureND<Int>.ndArray: INDArray
         get() = when (this) {
             is Nd4jArrayStructure<Int> -> ndArray
-            else -> Nd4j.zeros(*shape).also {
+            else -> Nd4j.zeros(*shape.asArray()).also {
                 elements().forEach { (idx, value) -> it.putScalar(idx, value) }
             }
         }
@@ -310,7 +313,7 @@ public open class IntNd4jArrayRingOps : Nd4jArrayRingOps<Int, IntRing> {
 
 public val IntRing.nd4j: IntNd4jArrayRingOps get() = IntNd4jArrayRingOps
 
-public class IntNd4jArrayRing(override val shape: Shape) : IntNd4jArrayRingOps(), RingND<Int, IntRing>
+public class IntNd4jArrayRing(override val shape: ShapeND) : IntNd4jArrayRingOps(), RingND<Int, IntRing>
 
 public fun IntRing.nd4j(shapeFirst: Int, vararg shapeRest: Int): IntNd4jArrayRing =
-    IntNd4jArrayRing(intArrayOf(shapeFirst, * shapeRest))
+    IntNd4jArrayRing(ShapeND(shapeFirst, * shapeRest))

kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jArrayStructure.kt

@@ -7,8 +7,7 @@ package space.kscience.kmath.nd4j
 
 import org.nd4j.linalg.api.ndarray.INDArray
 import space.kscience.kmath.misc.PerformancePitfall
-import space.kscience.kmath.nd.MutableStructureND
-import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.nd.*
 
 /**
  * Represents a [StructureND] wrapping an [INDArray] object.
@@ -22,7 +21,7 @@ public sealed class Nd4jArrayStructure<T> : MutableStructureND<T> {
      */
     public abstract val ndArray: INDArray
 
-    override val shape: IntArray get() = ndArray.shape().toIntArray()
+    override val shape: ShapeND get() = ShapeND(ndArray.shape().toIntArray())
 
     internal abstract fun elementsIterator(): Iterator<Pair<IntArray, T>>
     internal fun indicesIterator(): Iterator<IntArray> = ndArray.indicesIterator()
@@ -31,20 +30,31 @@ public sealed class Nd4jArrayStructure<T> : MutableStructureND<T> {
     override fun elements(): Sequence<Pair<IntArray, T>> = Sequence(::elementsIterator)
 }
 
-private data class Nd4jArrayIntStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Int>() {
+public data class Nd4jArrayIntStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Int>(), StructureNDOfInt {
     override fun elementsIterator(): Iterator<Pair<IntArray, Int>> = ndArray.intIterator()
+
+    @OptIn(PerformancePitfall::class)
     override fun get(index: IntArray): Int = ndArray.getInt(*index)
+
+    override fun getInt(index: IntArray): Int = ndArray.getInt(*index)
+
+    @OptIn(PerformancePitfall::class)
     override fun set(index: IntArray, value: Int): Unit = run { ndArray.putScalar(index, value) }
 }
 
 /**
  * Wraps this [INDArray] to [Nd4jArrayStructure].
  */
-public fun INDArray.asIntStructure(): Nd4jArrayStructure<Int> = Nd4jArrayIntStructure(this)
+public fun INDArray.asIntStructure(): Nd4jArrayIntStructure = Nd4jArrayIntStructure(this)
 
-private data class Nd4jArrayDoubleStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Double>() {
+public data class Nd4jArrayDoubleStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Double>(), StructureNDOfDouble {
     override fun elementsIterator(): Iterator<Pair<IntArray, Double>> = ndArray.realIterator()
+    @OptIn(PerformancePitfall::class)
     override fun get(index: IntArray): Double = ndArray.getDouble(*index)
+
+    override fun getDouble(index: IntArray): Double = ndArray.getDouble(*index)
+
+    @OptIn(PerformancePitfall::class)
     override fun set(index: IntArray, value: Double): Unit = run { ndArray.putScalar(index, value) }
 }
 
@@ -53,9 +63,12 @@ private data class Nd4jArrayDoubleStructure(override val ndArray: INDArray) : Nd
  */
 public fun INDArray.asDoubleStructure(): Nd4jArrayStructure<Double> = Nd4jArrayDoubleStructure(this)
 
-private data class Nd4jArrayFloatStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Float>() {
+public data class Nd4jArrayFloatStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Float>() {
     override fun elementsIterator(): Iterator<Pair<IntArray, Float>> = ndArray.floatIterator()
+    @PerformancePitfall
     override fun get(index: IntArray): Float = ndArray.getFloat(*index)
+
+    @PerformancePitfall
     override fun set(index: IntArray, value: Float): Unit = run { ndArray.putScalar(index, value) }
 }
 

kmath-nd4j/src/main/kotlin/space/kscience/kmath/nd4j/Nd4jTensorAlgebra.kt

@@ -13,9 +13,8 @@ import org.nd4j.linalg.factory.Nd4j
 import org.nd4j.linalg.factory.ops.NDBase
 import org.nd4j.linalg.ops.transforms.Transforms
 import space.kscience.kmath.misc.PerformancePitfall
-import space.kscience.kmath.nd.DefaultStrides
-import space.kscience.kmath.nd.Shape
-import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.misc.UnsafeKMathAPI
+import space.kscience.kmath.nd.*
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.Field
 import space.kscience.kmath.tensors.api.AnalyticTensorAlgebra
@@ -38,7 +37,7 @@ public sealed interface Nd4jTensorAlgebra<T : Number, A : Field<T>> : AnalyticTe
      */
     public val StructureND<T>.ndArray: INDArray
 
-    override fun structureND(shape: Shape, initializer: A.(IntArray) -> T): Nd4jArrayStructure<T>
+    override fun structureND(shape: ShapeND, initializer: A.(IntArray) -> T): Nd4jArrayStructure<T>
 
     @OptIn(PerformancePitfall::class)
     override fun StructureND<T>.map(transform: A.(T) -> T): Nd4jArrayStructure<T> =
@@ -96,7 +95,7 @@ public sealed interface Nd4jTensorAlgebra<T : Number, A : Field<T>> : AnalyticTe
 
     override fun StructureND<T>.unaryMinus(): Nd4jArrayStructure<T> = ndArray.neg().wrap()
     override fun Tensor<T>.getTensor(i: Int): Nd4jArrayStructure<T> = ndArray.slice(i.toLong()).wrap()
-    override fun Tensor<T>.transposed(i: Int, j: Int): Nd4jArrayStructure<T> = ndArray.swapAxes(i, j).wrap()
+    override fun StructureND<T>.transposed(i: Int, j: Int): Nd4jArrayStructure<T> = ndArray.swapAxes(i, j).wrap()
     override fun StructureND<T>.dot(other: StructureND<T>): Nd4jArrayStructure<T> = ndArray.mmul(other.ndArray).wrap()
 
     override fun StructureND<T>.min(dim: Int, keepDim: Boolean): Nd4jArrayStructure<T> =
@@ -108,7 +107,9 @@ public sealed interface Nd4jTensorAlgebra<T : Number, A : Field<T>> : AnalyticTe
     override fun StructureND<T>.max(dim: Int, keepDim: Boolean): Nd4jArrayStructure<T> =
         ndArray.max(keepDim, dim).wrap()
 
-    override fun Tensor<T>.view(shape: IntArray): Nd4jArrayStructure<T> = ndArray.reshape(shape).wrap()
+    @OptIn(UnsafeKMathAPI::class)
+    override fun Tensor<T>.view(shape: ShapeND): Nd4jArrayStructure<T> = ndArray.reshape(shape.asArray()).wrap()
+
     override fun Tensor<T>.viewAs(other: StructureND<T>): Nd4jArrayStructure<T> = view(other.shape)
 
     override fun StructureND<T>.argMin(dim: Int, keepDim: Boolean): Tensor<Int> =
@@ -117,35 +118,35 @@ public sealed interface Nd4jTensorAlgebra<T : Number, A : Field<T>> : AnalyticTe
     override fun StructureND<T>.argMax(dim: Int, keepDim: Boolean): Tensor<Int> =
         ndBase.get().argmax(ndArray, keepDim, dim).asIntStructure()
 
-    override fun StructureND<T>.mean(dim: Int, keepDim: Boolean): Nd4jArrayStructure<T> =
-        ndArray.mean(keepDim, dim).wrap()
+    override fun mean(structureND: StructureND<T>, dim: Int, keepDim: Boolean): Tensor<T> =
+        structureND.ndArray.mean(keepDim, dim).wrap()
 
-    override fun StructureND<T>.exp(): Nd4jArrayStructure<T> = Transforms.exp(ndArray).wrap()
-    override fun StructureND<T>.ln(): Nd4jArrayStructure<T> = Transforms.log(ndArray).wrap()
-    override fun StructureND<T>.sqrt(): Nd4jArrayStructure<T> = Transforms.sqrt(ndArray).wrap()
-    override fun StructureND<T>.cos(): Nd4jArrayStructure<T> = Transforms.cos(ndArray).wrap()
-    override fun StructureND<T>.acos(): Nd4jArrayStructure<T> = Transforms.acos(ndArray).wrap()
-    override fun StructureND<T>.cosh(): Nd4jArrayStructure<T> = Transforms.cosh(ndArray).wrap()
+    override fun exp(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.exp(arg.ndArray).wrap()
+    override fun ln(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.log(arg.ndArray).wrap()
+    override fun sqrt(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.sqrt(arg.ndArray).wrap()
+    override fun cos(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.cos(arg.ndArray).wrap()
+    override fun acos(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.acos(arg.ndArray).wrap()
+    override fun cosh(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.cosh(arg.ndArray).wrap()
 
-    override fun StructureND<T>.acosh(): Nd4jArrayStructure<T> =
-        Nd4j.getExecutioner().exec(ACosh(ndArray, ndArray.ulike())).wrap()
+    override fun acosh(arg: StructureND<T>): Nd4jArrayStructure<T> =
+        Nd4j.getExecutioner().exec(ACosh(arg.ndArray, arg.ndArray.ulike())).wrap()
 
-    override fun StructureND<T>.sin(): Nd4jArrayStructure<T> = Transforms.sin(ndArray).wrap()
-    override fun StructureND<T>.asin(): Nd4jArrayStructure<T> = Transforms.asin(ndArray).wrap()
-    override fun StructureND<T>.sinh(): Tensor<T> = Transforms.sinh(ndArray).wrap()
+    override fun sin(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.sin(arg.ndArray).wrap()
+    override fun asin(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.asin(arg.ndArray).wrap()
+    override fun sinh(arg: StructureND<T>): Tensor<T> = Transforms.sinh(arg.ndArray).wrap()
 
-    override fun StructureND<T>.asinh(): Nd4jArrayStructure<T> =
-        Nd4j.getExecutioner().exec(ASinh(ndArray, ndArray.ulike())).wrap()
+    override fun asinh(arg: StructureND<T>): Nd4jArrayStructure<T> =
+        Nd4j.getExecutioner().exec(ASinh(arg.ndArray, arg.ndArray.ulike())).wrap()
 
-    override fun StructureND<T>.tan(): Nd4jArrayStructure<T> = Transforms.tan(ndArray).wrap()
-    override fun StructureND<T>.atan(): Nd4jArrayStructure<T> = Transforms.atan(ndArray).wrap()
-    override fun StructureND<T>.tanh(): Nd4jArrayStructure<T> = Transforms.tanh(ndArray).wrap()
-    override fun StructureND<T>.atanh(): Nd4jArrayStructure<T> = Transforms.atanh(ndArray).wrap()
+    override fun tan(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.tan(arg.ndArray).wrap()
+    override fun atan(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.atan(arg.ndArray).wrap()
+    override fun tanh(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.tanh(arg.ndArray).wrap()
+    override fun atanh(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.atanh(arg.ndArray).wrap()
     override fun power(arg: StructureND<T>, pow: Number): StructureND<T> = Transforms.pow(arg.ndArray, pow).wrap()
-    override fun StructureND<T>.ceil(): Nd4jArrayStructure<T> = Transforms.ceil(ndArray).wrap()
-    override fun StructureND<T>.floor(): Nd4jArrayStructure<T> = Transforms.floor(ndArray).wrap()
-    override fun StructureND<T>.std(dim: Int, keepDim: Boolean): Nd4jArrayStructure<T> =
-        ndArray.std(true, keepDim, dim).wrap()
+    override fun ceil(arg: StructureND<T>): Nd4jArrayStructure<T> = Transforms.ceil(arg.ndArray).wrap()
+    override fun floor(structureND: StructureND<T>): Nd4jArrayStructure<T> = Transforms.floor(structureND.ndArray).wrap()
+    override fun std(structureND: StructureND<T>, dim: Int, keepDim: Boolean): Tensor<T> =
+        structureND.ndArray.std(true, keepDim, dim).wrap()
 
     override fun T.div(arg: StructureND<T>): Nd4jArrayStructure<T> = arg.ndArray.rdiv(this).wrap()
     override fun StructureND<T>.div(arg: T): Nd4jArrayStructure<T> = ndArray.div(arg).wrap()
@@ -159,8 +160,8 @@ public sealed interface Nd4jTensorAlgebra<T : Number, A : Field<T>> : AnalyticTe
         ndArray.divi(arg.ndArray)
     }
 
-    override fun StructureND<T>.variance(dim: Int, keepDim: Boolean): Nd4jArrayStructure<T> =
-        Nd4j.getExecutioner().exec(Variance(ndArray, true, true, dim)).wrap()
+    override fun variance(structureND: StructureND<T>, dim: Int, keepDim: Boolean): Tensor<T> =
+        Nd4j.getExecutioner().exec(Variance(structureND.ndArray, true, true, dim)).wrap()
 
     private companion object {
         private val ndBase: ThreadLocal<NDBase> = ThreadLocal.withInitial(::NDBase)
@@ -176,9 +177,10 @@ public object DoubleNd4jTensorAlgebra : Nd4jTensorAlgebra<Double, DoubleField> {
 
     override fun INDArray.wrap(): Nd4jArrayStructure<Double> = asDoubleStructure()
 
-    override fun structureND(shape: Shape, initializer: DoubleField.(IntArray) -> Double): Nd4jArrayStructure<Double> {
-        val array: INDArray = Nd4j.zeros(*shape)
-        val indices = DefaultStrides(shape)
+    @OptIn(UnsafeKMathAPI::class)
+    override fun structureND(shape: ShapeND, initializer: DoubleField.(IntArray) -> Double): Nd4jArrayStructure<Double> {
+        val array: INDArray = Nd4j.zeros(*shape.asArray())
+        val indices = ColumnStrides(shape)
         indices.asSequence().forEach { index ->
             array.putScalar(index, elementAlgebra.initializer(index))
         }
@@ -186,21 +188,21 @@ public object DoubleNd4jTensorAlgebra : Nd4jTensorAlgebra<Double, DoubleField> {
     }
 
 
-    @OptIn(PerformancePitfall::class)
+    @OptIn(PerformancePitfall::class, UnsafeKMathAPI::class)
     override val StructureND<Double>.ndArray: INDArray
         get() = when (this) {
             is Nd4jArrayStructure<Double> -> ndArray
-            else -> Nd4j.zeros(*shape).also {
+            else -> Nd4j.zeros(*shape.asArray()).also {
                 elements().forEach { (idx, value) -> it.putScalar(idx, value) }
             }
         }
 
     override fun StructureND<Double>.valueOrNull(): Double? =
-        if (shape contentEquals intArrayOf(1)) ndArray.getDouble(0) else null
+        if (shape contentEquals ShapeND(1)) ndArray.getDouble(0) else null
 
     // TODO rewrite
     override fun diagonalEmbedding(
-        diagonalEntries: Tensor<Double>,
+        diagonalEntries: StructureND<Double>,
         offset: Int,
         dim1: Int,
         dim2: Int,
@@ -209,7 +211,7 @@ public object DoubleNd4jTensorAlgebra : Nd4jTensorAlgebra<Double, DoubleField> {
     override fun StructureND<Double>.sum(): Double = ndArray.sumNumber().toDouble()
     override fun StructureND<Double>.min(): Double = ndArray.minNumber().toDouble()
     override fun StructureND<Double>.max(): Double = ndArray.maxNumber().toDouble()
-    override fun StructureND<Double>.mean(): Double = ndArray.meanNumber().toDouble()
-    override fun StructureND<Double>.std(): Double = ndArray.stdNumber().toDouble()
-    override fun StructureND<Double>.variance(): Double = ndArray.varNumber().toDouble()
+    override fun mean(structureND: StructureND<Double>): Double = structureND.ndArray.meanNumber().toDouble()
+    override fun std(structureND: StructureND<Double>): Double = structureND.ndArray.stdNumber().toDouble()
+    override fun variance(structureND: StructureND<Double>): Double = structureND.ndArray.varNumber().toDouble()
 }

kmath-nd4j/src/test/kotlin/space/kscience/kmath/nd4j/Nd4jArrayStructureTest.kt

@@ -7,6 +7,7 @@ package space.kscience.kmath.nd4j
 
 import org.nd4j.linalg.factory.Nd4j
 import space.kscience.kmath.misc.PerformancePitfall
+import space.kscience.kmath.nd.asList
 import space.kscience.kmath.nd.get
 import kotlin.test.Test
 import kotlin.test.assertEquals
@@ -27,7 +28,7 @@ internal class Nd4jArrayStructureTest {
     fun testShape() {
         val nd = Nd4j.rand(10, 2, 3, 6) ?: fail()
         val struct = nd.asDoubleStructure()
-        assertEquals(intArrayOf(10, 2, 3, 6).toList(), struct.shape.toList())
+        assertEquals(intArrayOf(10, 2, 3, 6).toList(), struct.shape.asList())
     }
 
     @Test

kmath-optimization/src/commonMain/tmp/minuit/HessianGradientCalculator.kt

@@ -33,7 +33,7 @@ internal class HessianGradientCalculator(fcn: MnFcn, par: MnUserTransformation,
         val g2: RealVector = gradient.getGradientDerivative()
         val gstep: RealVector = gradient.getStep()
         val fcnmin: Double = par.fval()
-        //   std::cout<<"fval: "<<fcnmin<<std::endl;
+        //   standardDiviation::cout<<"fval: "<<fcnmin<<standardDiviation::endl;
         val dfmin: Double = 4.0 * precision().eps2() * (abs(fcnmin) + theFcn.errorDef())
         val n: Int = x.getDimension()
         val dgrd: RealVector = ArrayRealVector(n)

kmath-optimization/src/commonMain/tmp/minuit/MnPosDef.kt

@@ -36,7 +36,7 @@ internal object MnPosDef {
         if (err.size() === 1 && err[0, 0] > prec.eps()) {
             return e
         }
-        //   std::cout<<"MnPosDef init matrix= "<<err<<std::endl;
+        //   standardDiviation::cout<<"MnPosDef init matrix= "<<err<<standardDiviation::endl;
         val epspdf: Double = max(1e-6, prec.eps2())
         var dgmin: Double = err[0, 0]
         for (i in 0 until err.nrow()) {
@@ -66,11 +66,11 @@ internal object MnPosDef {
             }
         }
 
-        //   std::cout<<"MnPosDef p: "<<p<<std::endl;
+        //   standardDiviation::cout<<"MnPosDef p: "<<p<<standardDiviation::endl;
         val eval: RealVector = p.eigenvalues()
         val pmin: Double = eval.getEntry(0)
         var pmax: Double = eval.getEntry(eval.getDimension() - 1)
-        //   std::cout<<"pmin= "<<pmin<<" pmax= "<<pmax<<std::endl;
+        //   standardDiviation::cout<<"pmin= "<<pmin<<" pmax= "<<pmax<<standardDiviation::endl;
         pmax = max(abs(pmax), 1.0)
         if (pmin > epspdf * pmax) {
             return e
@@ -81,9 +81,9 @@ internal object MnPosDef {
             err[i, i] = err[i, i] * (1.0 + padd)
             MINUITPlugin.logStatic(java.lang.Double.toString(eval.getEntry(i)))
         }
-        //   std::cout<<"MnPosDef final matrix: "<<err<<std::endl;
+        //   standardDiviation::cout<<"MnPosDef final matrix: "<<err<<standardDiviation::endl;
         MINUITPlugin.logStatic("matrix forced pos-def by adding $padd to diagonal")
-        //   std::cout<<"eigenvalues: "<<eval<<std::endl;
+        //   standardDiviation::cout<<"eigenvalues: "<<eval<<standardDiviation::endl;
         return MinimumError(err, MnMadePosDef())
     }
 }

kmath-polynomial/README.md

@@ -1,50 +0,0 @@
-# Module kmath-polynomial
-
-Polynomials, rational functions, and utilities
-
-## Features
-
- - [polynomial abstraction](src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt) : Abstraction for polynomial spaces.
- - [rational function abstraction](src/commonMain/kotlin/space/kscience/kmath/functions/RationalFunction.kt) : Abstraction for rational functions spaces.
- - ["list" polynomials](src/commonMain/kotlin/space/kscience/kmath/functions/ListRationalFunction.kt) : List implementation of univariate polynomials.
- - ["list" rational functions](src/commonMain/kotlin/space/kscience/kmath/functions/ListPolynomial.kt) : List implementation of univariate rational functions.
- - ["list" polynomials and rational functions constructors](src/commonMain/kotlin/space/kscience/kmath/functions/listConstructors.kt) : Constructors for list polynomials and rational functions.
- - ["list" polynomials and rational functions utilities](src/commonMain/kotlin/space/kscience/kmath/functions/listUtil.kt) : Utilities for list polynomials and rational functions.
- - ["numbered" polynomials](src/commonMain/kotlin/space/kscience/kmath/functions/NumberedRationalFunction.kt) : Numbered implementation of multivariate polynomials.
- - ["numbered" rational functions](src/commonMain/kotlin/space/kscience/kmath/functions/NumberedPolynomial.kt) : Numbered implementation of multivariate rational functions.
- - ["numbered" polynomials and rational functions constructors](src/commonMain/kotlin/space/kscience/kmath/functions/numberedConstructors.kt) : Constructors for numbered polynomials and rational functions.
- - ["numbered" polynomials and rational functions utilities](src/commonMain/kotlin/space/kscience/kmath/functions/numberedUtil.kt) : Utilities for numbered polynomials and rational functions.
- - ["labeled" polynomials](src/commonMain/kotlin/space/kscience/kmath/functions/LabeledRationalFunction.kt) : Labeled implementation of multivariate polynomials.
- - ["labeled" rational functions](src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt) : Labeled implementation of multivariate rational functions.
- - ["labeled" polynomials and rational functions constructors](src/commonMain/kotlin/space/kscience/kmath/functions/labeledConstructors.kt) : Constructors for labeled polynomials and rational functions.
- - ["labeled" polynomials and rational functions utilities](src/commonMain/kotlin/space/kscience/kmath/functions/labeledUtil.kt) : Utilities for labeled polynomials and rational functions.
-
-
-## Usage
-
-## Artifact:
-
-The Maven coordinates of this project are `space.kscience:kmath-polynomial:0.3.1-dev-1`.
-
-**Gradle Groovy:**
-```groovy
-repositories {
-    maven { url 'https://repo.kotlin.link' }
-    mavenCentral()
-}
-
-dependencies {
-    implementation 'space.kscience:kmath-polynomial:0.3.1-dev-1'
-}
-```
-**Gradle Kotlin DSL:**
-```kotlin
-repositories {
-    maven("https://repo.kotlin.link")
-    mavenCentral()
-}
-
-dependencies {
-    implementation("space.kscience:kmath-polynomial:0.3.1-dev-1")
-}
-```

kmath-polynomial/build.gradle.kts

@@ -1,69 +0,0 @@
-plugins {
-    id("space.kscience.gradle.mpp")
-}
-
-kscience{
-    native()
-}
-
-description = "Polynomials, rational functions, and utilities"
-
-kotlin.sourceSets {
-    commonMain {
-        dependencies {
-            api(projects.kmathCore)
-        }
-    }
-}
-
-dependencies {
-    dokkaPlugin("org.jetbrains.dokka:mathjax-plugin:${npmlibs.versions.dokka.get()}")
-}
-
-readme {
-    maturity = space.kscience.gradle.Maturity.PROTOTYPE
-    propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md"))
-
-    feature("polynomial abstraction", "src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt") {
-        "Abstraction for polynomial spaces."
-    }
-    feature("rational function abstraction", "src/commonMain/kotlin/space/kscience/kmath/functions/RationalFunction.kt") {
-        "Abstraction for rational functions spaces."
-    }
-    feature("\"list\" polynomials", "src/commonMain/kotlin/space/kscience/kmath/functions/ListRationalFunction.kt") {
-        "List implementation of univariate polynomials."
-    }
-    feature("\"list\" rational functions", "src/commonMain/kotlin/space/kscience/kmath/functions/ListPolynomial.kt") {
-        "List implementation of univariate rational functions."
-    }
-    feature("\"list\" polynomials and rational functions constructors", "src/commonMain/kotlin/space/kscience/kmath/functions/listConstructors.kt") {
-        "Constructors for list polynomials and rational functions."
-    }
-    feature("\"list\" polynomials and rational functions utilities", "src/commonMain/kotlin/space/kscience/kmath/functions/listUtil.kt") {
-        "Utilities for list polynomials and rational functions."
-    }
-    feature("\"numbered\" polynomials", "src/commonMain/kotlin/space/kscience/kmath/functions/NumberedRationalFunction.kt") {
-        "Numbered implementation of multivariate polynomials."
-    }
-    feature("\"numbered\" rational functions", "src/commonMain/kotlin/space/kscience/kmath/functions/NumberedPolynomial.kt") {
-        "Numbered implementation of multivariate rational functions."
-    }
-    feature("\"numbered\" polynomials and rational functions constructors", "src/commonMain/kotlin/space/kscience/kmath/functions/numberedConstructors.kt") {
-        "Constructors for numbered polynomials and rational functions."
-    }
-    feature("\"numbered\" polynomials and rational functions utilities", "src/commonMain/kotlin/space/kscience/kmath/functions/numberedUtil.kt") {
-        "Utilities for numbered polynomials and rational functions."
-    }
-    feature("\"labeled\" polynomials", "src/commonMain/kotlin/space/kscience/kmath/functions/LabeledRationalFunction.kt") {
-        "Labeled implementation of multivariate polynomials."
-    }
-    feature("\"labeled\" rational functions", "src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt") {
-        "Labeled implementation of multivariate rational functions."
-    }
-    feature("\"labeled\" polynomials and rational functions constructors", "src/commonMain/kotlin/space/kscience/kmath/functions/labeledConstructors.kt") {
-        "Constructors for labeled polynomials and rational functions."
-    }
-    feature("\"labeled\" polynomials and rational functions utilities", "src/commonMain/kotlin/space/kscience/kmath/functions/labeledUtil.kt") {
-        "Utilities for labeled polynomials and rational functions."
-    }
-}

kmath-polynomial/src/commonMain/kotlin/space/kscience/kmath/functions/LabeledPolynomial.kt

@@ -1,529 +0,0 @@
-/*
- * Copyright 2018-2022 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("NOTHING_TO_INLINE", "KotlinRedundantDiagnosticSuppress")
-
-package space.kscience.kmath.functions
-
-import space.kscience.kmath.expressions.Symbol
-import space.kscience.kmath.operations.Ring
-import kotlin.jvm.JvmName
-import kotlin.math.max
-
-
-/**
- * Represents multivariate polynomial that stores its coefficients in a [Map] and terms' signatures in a [Map] that
- * associates variables (of type [Symbol]) with their degree.
- *
- * @param C the type of constants.
- */
-public data class LabeledPolynomial<out C>
-@PublishedApi
-internal constructor(
-    /**
-     * Map that contains coefficients of the polynomial.
-     *
-     * Every monomial \(a x_1^{d_1} ... x_n^{d_n}\) is stored as a pair "key-value" in the map, where the value is the
-     * coefficient \(a\) and the key is a map that associates variables in the monomial with their degree in the monomial.
-     * For example, coefficients of a polynomial \(5 a^2 c^3 - 6 b\) can be represented as
-     * ```
-     * mapOf(
-     *      mapOf(
-     *          a to 2,
-     *          c to 3
-     *      ) to 5,
-     *      mapOf(
-     *          b to 1
-     *      ) to (-6)
-     * )
-     * ```
-     * and also as
-     * ```
-     * mapOf(
-     *      mapOf(
-     *          a to 2,
-     *          c to 3
-     *      ) to 5,
-     *      mapOf(
-     *          b to 1
-     *      ) to (-6),
-     *      mapOf(
-     *          b to 1,
-     *          c to 1
-     *      ) to 0
-     * )
-     * ```
-     * where \(a\), \(b\) and \(c\) are corresponding [Symbol] objects.
-     *
-     * It is not prohibited to put extra zero monomials into the map (as for \(0 b c\) in the example). But the
-     * bigger the coefficients map the worse performance of arithmetical operations performed on it. Thus, it is
-     * recommended not to put (or even to remove) extra (or useless) monomials in the coefficients map.
-     * @usesMathJax
-     */
-    public val coefficients: Map<Map<Symbol, UInt>, C>
-) {
-    override fun toString(): String = "LabeledPolynomial$coefficients"
-}
-
-/**
- * Arithmetic context for multivariate polynomials with coefficients stored as a [Map] and terms' signatures stored as a
- * [Map] constructed with the provided [ring] of constants.
- *
- * @param C the type of constants. Polynomials have them a coefficients in their terms.
- * @param A type of provided underlying ring of constants. It's [Ring] of [C].
- * @param ring underlying ring of constants of type [A].
- */
-public class LabeledPolynomialSpace<C, out A : Ring<C>>(
-    public override val ring: A,
-) : MultivariatePolynomialSpace<C, Symbol, LabeledPolynomial<C>>, PolynomialSpaceOverRing<C, LabeledPolynomial<C>, A> {
-    /**
-     * Returns sum of the variable represented as a monic monomial and the integer represented as a constant polynomial.
-     */
-    public override operator fun Symbol.plus(other: Int): LabeledPolynomial<C> =
-        if (other == 0) LabeledPolynomialAsIs(
-            mapOf(this@plus to 1U) to constantOne,
-        )
-        else LabeledPolynomialAsIs(
-            mapOf(this@plus to 1U) to constantOne,
-            emptyMap<Symbol, UInt>() to other.asConstant(),
-        )
-    /**
-     * Returns difference between the variable represented as a monic monomial and the integer represented as a constant polynomial.
-     */
-    public override operator fun Symbol.minus(other: Int): LabeledPolynomial<C> =
-        if (other == 0) LabeledPolynomialAsIs(
-            mapOf(this@minus to 1U) to constantOne,
-        )
-        else LabeledPolynomialAsIs(
-            mapOf(this@minus to 1U) to constantOne,
-            emptyMap<Symbol, UInt>() to (-other).asConstant(),
-        )
-    /**
-     * Returns product of the variable represented as a monic monomial and the integer represented as a constant polynomial.
-     */
-    public override operator fun Symbol.times(other: Int): LabeledPolynomial<C> =
-        if (other == 0) zero
-        else LabeledPolynomialAsIs(
-            mapOf(this to 1U) to other.asConstant(),
-        )
-
-    /**
-     * Returns sum of the integer represented as a constant polynomial and the variable represented as a monic monomial.
-     */
-    public override operator fun Int.plus(other: Symbol): LabeledPolynomial<C> =
-        if (this == 0) LabeledPolynomialAsIs(
-            mapOf(other to 1U) to constantOne,
-        )
-        else LabeledPolynomialAsIs(
-            mapOf(other to 1U) to constantOne,
-            emptyMap<Symbol, UInt>() to this@plus.asConstant(),
-        )
-    /**
-     * Returns difference between the integer represented as a constant polynomial and the variable represented as a monic monomial.
-     */
-    public override operator fun Int.minus(other: Symbol): LabeledPolynomial<C> =
-        if (this == 0) LabeledPolynomialAsIs(
-            mapOf(other to 1U) to -constantOne,
-        )
-        else LabeledPolynomialAsIs(
-            mapOf(other to 1U) to -constantOne,
-            emptyMap<Symbol, UInt>() to constantOne * this@minus,
-        )
-    /**
-     * Returns product of the integer represented as a constant polynomial and the variable represented as a monic monomial.
-     */
-    public override operator fun Int.times(other: Symbol): LabeledPolynomial<C> =
-        if (this == 0) zero
-        else LabeledPolynomialAsIs(
-            mapOf(other to 1U) to this@times.asConstant(),
-        )
-
-    /**
-     * Returns sum of the polynomial and the integer represented as a polynomial.
-     *
-     * The operation is equivalent to adding [other] copies of unit polynomial to [this].
-     */
-    public override operator fun LabeledPolynomial<C>.plus(other: Int): LabeledPolynomial<C> =
-        when {
-            other == 0 -> this
-            coefficients.isEmpty() -> other.asPolynomial()
-            else -> LabeledPolynomialAsIs(
-                coefficients.withPutOrChanged(emptyMap(), other.asConstant()) { it -> it + other }
-            )
-        }
-    /**
-     * Returns difference between the polynomial and the integer represented as a polynomial.
-     *
-     * The operation is equivalent to subtraction [other] copies of unit polynomial from [this].
-     */
-    public override operator fun LabeledPolynomial<C>.minus(other: Int): LabeledPolynomial<C> =
-        when {
-            other == 0 -> this
-            coefficients.isEmpty() -> other.asPolynomial()
-            else -> LabeledPolynomialAsIs(
-                coefficients.withPutOrChanged(emptyMap(), (-other).asConstant()) { it -> it - other }
-            )
-        }
-    /**
-     * Returns product of the polynomial and the integer represented as a polynomial.
-     *
-     * The operation is equivalent to sum of [other] copies of [this].
-     */
-    public override operator fun LabeledPolynomial<C>.times(other: Int): LabeledPolynomial<C> =
-        when(other) {
-            0 -> zero
-            1 -> this
-            else -> LabeledPolynomialAsIs(
-                coefficients.mapValues { (_, value) -> value * other }
-            )
-        }
-
-    /**
-     * Returns sum of the integer represented as a polynomial and the polynomial.
-     *
-     * The operation is equivalent to adding [this] copies of unit polynomial to [other].
-     */
-    public override operator fun Int.plus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        when {
-            this == 0 -> other
-            other.coefficients.isEmpty() -> this@plus.asPolynomial()
-            else -> LabeledPolynomialAsIs(
-                other.coefficients.withPutOrChanged(emptyMap(), this@plus.asConstant()) { it -> this@plus + it }
-            )
-        }
-    /**
-     * Returns difference between the integer represented as a polynomial and the polynomial.
-     *
-     * The operation is equivalent to subtraction [this] copies of unit polynomial from [other].
-     */
-    public override operator fun Int.minus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        when {
-            this == 0 -> -other
-            other.coefficients.isEmpty() -> this@minus.asPolynomial()
-            else -> LabeledPolynomialAsIs(
-                buildMap(other.coefficients.size + 1) {
-                    put(emptyMap(), asConstant())
-                    other.coefficients.copyMapToBy(this, { _, c -> -c }, { currentC, newC -> currentC - newC })
-                }
-            )
-        }
-    /**
-     * Returns product of the integer represented as a polynomial and the polynomial.
-     *
-     * The operation is equivalent to sum of [this] copies of [other].
-     */
-    public override operator fun Int.times(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        when(this) {
-            0 -> zero
-            1 -> other
-            else -> LabeledPolynomialAsIs(
-                other.coefficients.mapValues { (_, value) -> this@times * value }
-            )
-        }
-
-    /**
-     * Returns sum of the variable represented as a monic monomial and the constant represented as a constant polynomial.
-     */
-    public override operator fun Symbol.plus(other: C): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            mapOf(this@plus to 1U) to constantOne,
-            emptyMap<Symbol, UInt>() to other,
-        )
-    /**
-     * Returns difference between the variable represented as a monic monomial and the constant represented as a constant polynomial.
-     */
-    public override operator fun Symbol.minus(other: C): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            mapOf(this@minus to 1U) to constantOne,
-            emptyMap<Symbol, UInt>() to -other,
-        )
-    /**
-     * Returns product of the variable represented as a monic monomial and the constant represented as a constant polynomial.
-     */
-    public override operator fun Symbol.times(other: C): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            mapOf(this@times to 1U) to other,
-        )
-
-    /**
-     * Returns sum of the constant represented as a constant polynomial and the variable represented as a monic monomial.
-     */
-    public override operator fun C.plus(other: Symbol): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            mapOf(other to 1U) to constantOne,
-            emptyMap<Symbol, UInt>() to this@plus,
-        )
-    /**
-     * Returns difference between the constant represented as a constant polynomial and the variable represented as a monic monomial.
-     */
-    public override operator fun C.minus(other: Symbol): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            mapOf(other to 1U) to -constantOne,
-            emptyMap<Symbol, UInt>() to this@minus,
-        )
-    /**
-     * Returns product of the constant represented as a constant polynomial and the variable represented as a monic monomial.
-     */
-    public override operator fun C.times(other: Symbol): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            mapOf(other to 1U) to this@times,
-        )
-
-    /**
-     * Returns sum of the constant represented as a polynomial and the polynomial.
-     */
-    override operator fun C.plus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        if (other.coefficients.isEmpty()) this@plus.asLabeledPolynomial()
-        else LabeledPolynomialAsIs(
-            other.coefficients.withPutOrChanged(emptyMap(), this@plus) { it -> this@plus + it }
-        )
-    /**
-     * Returns difference between the constant represented as a polynomial and the polynomial.
-     */
-    override operator fun C.minus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        if (other.coefficients.isEmpty()) this@minus.asPolynomial()
-        else LabeledPolynomialAsIs(
-            buildMap(other.coefficients.size + 1) {
-                put(emptyMap(), this@minus)
-                other.coefficients.copyMapToBy(this, { _, c -> -c }, { currentC, newC -> currentC - newC })
-            }
-        )
-    /**
-     * Returns product of the constant represented as a polynomial and the polynomial.
-     */
-    override operator fun C.times(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            other.coefficients.mapValues { this@times * it.value }
-        )
-
-    /**
-     * Returns sum of the constant represented as a polynomial and the polynomial.
-     */
-    override operator fun LabeledPolynomial<C>.plus(other: C): LabeledPolynomial<C> =
-        if (coefficients.isEmpty()) other.asLabeledPolynomial()
-        else LabeledPolynomialAsIs(
-            coefficients.withPutOrChanged(emptyMap(), other) { it -> it + other }
-        )
-    /**
-     * Returns difference between the constant represented as a polynomial and the polynomial.
-     */
-    override operator fun LabeledPolynomial<C>.minus(other: C): LabeledPolynomial<C> =
-        if (coefficients.isEmpty()) other.asLabeledPolynomial()
-        else LabeledPolynomialAsIs(
-            coefficients.withPutOrChanged(emptyMap(), -other) { it -> it - other }
-        )
-    /**
-     * Returns product of the constant represented as a polynomial and the polynomial.
-     */
-    override operator fun LabeledPolynomial<C>.times(other: C): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            coefficients.mapValues { it.value * other }
-        )
-
-    /**
-     * Converts the constant [value] to polynomial.
-     */
-    public override fun number(value: C): LabeledPolynomial<C> = value.asLabeledPolynomial()
-
-    /**
-     * Represents the variable as a monic monomial.
-     */
-    public override operator fun Symbol.unaryPlus(): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            mapOf(this to 1U) to constantOne,
-        )
-    /**
-     * Returns negation of representation of the variable as a monic monomial.
-     */
-    public override operator fun Symbol.unaryMinus(): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            mapOf(this to 1U) to -constantOne,
-        )
-    /**
-     * Returns sum of the variables represented as monic monomials.
-     */
-    public override operator fun Symbol.plus(other: Symbol): LabeledPolynomial<C> =
-        if (this == other) LabeledPolynomialAsIs(
-            mapOf(this to 1U) to constantOne * 2
-        )
-        else LabeledPolynomialAsIs(
-            mapOf(this to 1U) to constantOne,
-            mapOf(other to 1U) to constantOne,
-        )
-    /**
-     * Returns difference between the variables represented as monic monomials.
-     */
-    public override operator fun Symbol.minus(other: Symbol): LabeledPolynomial<C> =
-        if (this == other) zero
-        else LabeledPolynomialAsIs(
-            mapOf(this to 1U) to constantOne,
-            mapOf(other to 1U) to -constantOne,
-        )
-    /**
-     * Returns product of the variables represented as monic monomials.
-     */
-    public override operator fun Symbol.times(other: Symbol): LabeledPolynomial<C> =
-        if (this == other) LabeledPolynomialAsIs(
-            mapOf(this to 2U) to constantOne
-        )
-        else LabeledPolynomialAsIs(
-            mapOf(this to 1U, other to 1U) to constantOne,
-        )
-
-    /**
-     * Returns sum of the variable represented as a monic monomial and the polynomial.
-     */
-    public override operator fun Symbol.plus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        if (other.coefficients.isEmpty()) this@plus.asPolynomial()
-        else LabeledPolynomialAsIs(
-            other.coefficients.withPutOrChanged(mapOf(this@plus to 1U), constantOne) { it -> constantOne + it }
-        )
-    /**
-     * Returns difference between the variable represented as a monic monomial and the polynomial.
-     */
-    public override operator fun Symbol.minus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        if (other.coefficients.isEmpty()) this@minus.asPolynomial()
-        else LabeledPolynomialAsIs(
-            buildMap(other.coefficients.size + 1) {
-                put(mapOf(this@minus to 1U), constantOne)
-                other.coefficients.copyMapToBy(this, { _, c -> -c }) { currentC, newC -> currentC - newC }
-            }
-        )
-    /**
-     * Returns product of the variable represented as a monic monomial and the polynomial.
-     */
-    public override operator fun Symbol.times(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            other.coefficients
-                .mapKeys { (degs, _) -> degs.withPutOrChanged(this, 1u) { it -> it + 1u } }
-        )
-
-    /**
-     * Returns sum of the polynomial and the variable represented as a monic monomial.
-     */
-    public override operator fun LabeledPolynomial<C>.plus(other: Symbol): LabeledPolynomial<C> =
-        if (coefficients.isEmpty()) other.asPolynomial()
-        else LabeledPolynomialAsIs(
-            coefficients.withPutOrChanged(mapOf(other to 1U), constantOne) { it -> it + constantOne }
-        )
-    /**
-     * Returns difference between the polynomial and the variable represented as a monic monomial.
-     */
-    public override operator fun LabeledPolynomial<C>.minus(other: Symbol): LabeledPolynomial<C> =
-        if (coefficients.isEmpty()) other.asPolynomial()
-        else LabeledPolynomialAsIs(
-            coefficients.withPutOrChanged(mapOf(other to 1U), -constantOne) { it -> it - constantOne }
-        )
-    /**
-     * Returns product of the polynomial and the variable represented as a monic monomial.
-     */
-    public override operator fun LabeledPolynomial<C>.times(other: Symbol): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            coefficients
-                .mapKeys { (degs, _) -> degs.withPutOrChanged(other, 1u) { it -> it + 1u } }
-        )
-
-    /**
-     * Returns negation of the polynomial.
-     */
-    override fun LabeledPolynomial<C>.unaryMinus(): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            coefficients.mapValues { -it.value }
-        )
-    /**
-     * Returns sum of the polynomials.
-     */
-    override operator fun LabeledPolynomial<C>.plus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            mergeBy(coefficients, other.coefficients) { c1, c2 -> c1 + c2 }
-        )
-    /**
-     * Returns difference of the polynomials.
-     */
-    override operator fun LabeledPolynomial<C>.minus(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            buildMap(coefficients.size + other.coefficients.size) {
-                coefficients.copyTo(this)
-                other.coefficients.copyMapToBy(this, { _, c -> -c }, { currentC, newC -> currentC - newC })
-            }
-        )
-    /**
-     * Returns product of the polynomials.
-     */
-    override operator fun LabeledPolynomial<C>.times(other: LabeledPolynomial<C>): LabeledPolynomial<C> =
-        LabeledPolynomialAsIs(
-            buildMap(coefficients.size * other.coefficients.size) {
-                for ((degs1, c1) in coefficients) for ((degs2, c2) in other.coefficients) {
-                    val degs = mergeBy(degs1, degs2) { deg1, deg2 -> deg1 + deg2 }
-                    val c = c1 * c2
-                    this.putOrChange(degs, c) { it -> it + c }
-                }
-            }
-        )
-
-    /**
-     * Instance of zero polynomial (zero of the polynomial ring).
-     */
-    override val zero: LabeledPolynomial<C> = LabeledPolynomialAsIs()
-    /**
-     * Instance of unit polynomial (unit of the polynomial ring).
-     */
-    override val one: LabeledPolynomial<C> = constantOne.asLabeledPolynomial()
-
-    /**
-     * Degree of the polynomial, [see also](https://en.wikipedia.org/wiki/Degree_of_a_polynomial). If the polynomial is
-     * zero, degree is -1.
-     */
-    override val LabeledPolynomial<C>.degree: Int
-        get() = coefficients.entries.maxOfOrNull { (degs, _) -> degs.values.sum().toInt() } ?: -1
-    /**
-     * Map that associates variables (that appear in the polynomial in positive exponents) with their most exponents
-     * in which they are appeared in the polynomial.
-     *
-     * As consequence all values in the map are positive integers. Also, if the polynomial is constant, the map is empty.
-     * And keys of the map is the same as in [variables].
-     */
-    public override val LabeledPolynomial<C>.degrees: Map<Symbol, UInt>
-        get() =
-            buildMap {
-                coefficients.keys.forEach { degs ->
-                    degs.copyToBy(this, ::max)
-                }
-            }
-    /**
-     * Counts degree of the polynomial by the specified [variable].
-     */
-    public override fun LabeledPolynomial<C>.degreeBy(variable: Symbol): UInt =
-        coefficients.entries.maxOfOrNull { (degs, _) -> degs.getOrElse(variable) { 0u } } ?: 0u
-    /**
-     * Counts degree of the polynomial by the specified [variables].
-     */
-    public override fun LabeledPolynomial<C>.degreeBy(variables: Collection<Symbol>): UInt =
-        coefficients.entries.maxOfOrNull { (degs, _) -> degs.filterKeys { it in variables }.values.sum() } ?: 0u
-    /**
-     * Set of all variables that appear in the polynomial in positive exponents.
-     */
-    public override val LabeledPolynomial<C>.variables: Set<Symbol>
-        get() =
-            buildSet {
-                coefficients.entries.forEach { (degs, _) -> addAll(degs.keys) }
-            }
-    /**
-     * Count of all variables that appear in the polynomial in positive exponents.
-     */
-    public override val LabeledPolynomial<C>.countOfVariables: Int get() = variables.size
-
-    // TODO: When context receivers will be ready move all of this substitutions and invocations to utilities with
-    //  [ListPolynomialSpace] as a context receiver
-    /**
-     * Substitutes provided arguments [arguments] into [this] polynomial.
-     */
-    public inline fun LabeledPolynomial<C>.substitute(arguments: Map<Symbol, C>): LabeledPolynomial<C> = substitute(ring, arguments)
-    /**
-     * Substitutes provided arguments [arguments] into [this] polynomial.
-     */
-    @JvmName("substitutePolynomial")
-    public inline fun LabeledPolynomial<C>.substitute(arguments: Map<Symbol, LabeledPolynomial<C>>) : LabeledPolynomial<C> = substitute(ring, arguments)
-}