From 061398b009ae4819adbe95c570bc95fccfb46f82 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Fri, 22 Jan 2021 23:22:24 +0300
Subject: [PATCH] [WIP] Refactor NDStructures

---
 CHANGELOG.md                                  |   1 +
 README.md                                     |   2 +-
 build.gradle.kts                              |   2 +-
 .../kscience/kmath/benchmarks/DotBenchmark.kt |   4 +-
 .../benchmarks/LinearAlgebraBenchmark.kt      |   1 -
 .../kmath/benchmarks/NDFieldBenchmark.kt      |   9 +-
 .../kmath/benchmarks/ViktorBenchmark.kt       |   4 +-
 .../kscience/kmath/operations/ComplexDemo.kt  |   3 +-
 .../kscience/kmath/structures/ComplexND.kt    |  13 +-
 .../kscience/kmath/structures/NDField.kt      |  17 +-
 .../structures/StructureReadBenchmark.kt      |   4 +-
 .../structures/StructureWriteBenchmark.kt     |   2 +
 .../kscience/kmath/commons/linear/CMMatrix.kt |   3 +-
 .../kscience/kmath/commons/linear/CMSolver.kt |   2 +-
 kmath-core/README.md                          |   2 +-
 .../kmath/expressions/SymbolIndexer.kt        |   2 +-
 .../kscience/kmath/linear/BufferMatrix.kt     |  13 +-
 .../kscience/kmath/linear/LinearAlgebra.kt    |   1 -
 .../kscience/kmath/linear/LupDecomposition.kt |   1 +
 .../kscience/kmath/linear/MatrixBuilder.kt    |   9 +-
 .../kscience/kmath/linear/MatrixContext.kt    |   1 -
 .../kscience/kmath/linear/MatrixFeatures.kt   |   2 -
 .../kscience/kmath/linear/MatrixWrapper.kt    |   8 +-
 .../kmath/linear/RealMatrixContext.kt         |   1 -
 .../kscience/kmath/linear/VirtualMatrix.kt    |   2 -
 .../kscience/kmath/nd/BufferNDAlgebra.kt      | 128 +++++++++
 .../kscience/kmath/nd/ComplexNDField.kt       | 113 ++++++++
 .../kotlin/kscience/kmath/nd/NDAlgebra.kt     | 262 ++++++++++++++++++
 .../kmath/{structures => nd}/NDStructure.kt   |  96 ++++---
 .../kotlin/kscience/kmath/nd/RealNDField.kt   | 111 ++++++++
 .../kotlin/kscience/kmath/nd/ShortNDRing.kt   |  36 +++
 .../kmath/{structures => nd}/Structure1D.kt   |   7 +-
 .../kmath/{structures => nd}/Structure2D.kt   |  12 +-
 .../kscience/kmath/operations/BigInt.kt       |   1 +
 .../kmath/structures/BoxingNDField.kt         |  81 ------
 .../kscience/kmath/structures/BoxingNDRing.kt |  71 -----
 .../kmath/structures/BufferAccessor2D.kt      |   5 +
 .../kmath/structures/BufferedNDAlgebra.kt     |  43 ---
 .../kmath/structures/BufferedNDElement.kt     |  86 ------
 .../kscience/kmath/structures/Buffers.kt      |   2 +-
 .../kmath/structures/ComplexNDField.kt        | 158 -----------
 .../kmath/structures/ExtendedNDField.kt       |  44 ---
 .../kmath/structures/FlaggedBuffer.kt         |   2 +-
 .../kscience/kmath/structures/MemoryBuffer.kt |   2 +-
 .../kscience/kmath/structures/NDAlgebra.kt    | 259 -----------------
 .../kscience/kmath/structures/NDElement.kt    | 134 ---------
 .../kscience/kmath/structures/RealNDField.kt  | 140 ----------
 .../kscience/kmath/structures/ShortNDRing.kt  |  93 -------
 .../kscience/kmath/linear/MatrixTest.kt       |   5 +-
 .../kscience/kmath/linear/RealLUSolverTest.kt |   1 -
 .../kscience/kmath/structures/NDFieldTest.kt  |   1 +
 .../kmath/structures/NumberNDFieldTest.kt     |   3 +-
 .../kmath/structures/LazyNDStructure.kt       |   4 +-
 .../kscience/kmath/dimensions/Wrappers.kt     |   3 +-
 .../kotlin/kscience/kmath/ejml/EjmlMatrix.kt  |   3 +-
 .../kscience/kmath/ejml/EjmlMatrixContext.kt  |   5 +-
 .../kscience/kmath/ejml/EjmlMatrixTest.kt     |   2 +-
 .../kotlin/kscience/kmath/real/RealMatrix.kt  |   6 +-
 .../kotlin/kscience/kmath/real/realND.kt      |  31 +++
 .../kaceince/kmath/real/RealMatrixTest.kt     |   4 +-
 .../kmath/interpolation/XYPointSet.kt         |   2 +-
 .../kscience/kmath/histogram/RealHistogram.kt |   8 +-
 .../kscience.kmath.nd4j/Nd4jArrayAlgebra.kt   | 132 ++++-----
 .../kscience.kmath.nd4j/Nd4jArrayStructure.kt |   4 +-
 .../kmath/nd4j/Nd4jArrayStructureTest.kt      |   2 +-
 .../kmath/viktor/ViktorNDStructure.kt         |   8 +-
 66 files changed, 889 insertions(+), 1330 deletions(-)
 create mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/nd/BufferNDAlgebra.kt
 create mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/nd/ComplexNDField.kt
 create mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDAlgebra.kt
 rename kmath-core/src/commonMain/kotlin/kscience/kmath/{structures => nd}/NDStructure.kt (84%)
 create mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/nd/RealNDField.kt
 create mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/nd/ShortNDRing.kt
 rename kmath-core/src/commonMain/kotlin/kscience/kmath/{structures => nd}/Structure1D.kt (89%)
 rename kmath-core/src/commonMain/kotlin/kscience/kmath/{structures => nd}/Structure2D.kt (92%)
 delete mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BoxingNDField.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BoxingNDRing.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferedNDAlgebra.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferedNDElement.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ComplexNDField.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ExtendedNDField.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDAlgebra.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDElement.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/structures/RealNDField.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ShortNDRing.kt
 create mode 100644 kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/realND.kt

diff --git a/CHANGELOG.md b/CHANGELOG.md
index aa70e6116..7536690b2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -35,6 +35,7 @@
 - Matrix LUP `inverse` renamed to `inverseWithLUP`
 - `NumericAlgebra` moved outside of regular algebra chain (`Ring` no longer implements it).
 - Features moved to NDStructure and became transparent.
+- Refactored `NDStructure` algebra to be more simple, preferring under-the-hood conversion to explicit NDStructure types
 
 ### Deprecated
 
diff --git a/README.md b/README.md
index 0899f77cc..258d79be5 100644
--- a/README.md
+++ b/README.md
@@ -114,7 +114,7 @@ submit a feature request if you want something to be implemented first.
 >
 > **Features:**
 > - [algebras](kmath-core/src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt) : Algebraic structures: contexts and elements
-> - [nd](kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt) : Many-dimensional structures
+> - [nd](kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDStructure.kt) : Many-dimensional structures
 > - [buffers](kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : One-dimensional structure
 > - [expressions](kmath-core/src/commonMain/kotlin/kscience/kmath/expressions) : Functional Expressions
 > - [domains](kmath-core/src/commonMain/kotlin/kscience/kmath/domains) : Domains
diff --git a/build.gradle.kts b/build.gradle.kts
index d171bd608..0572217af 100644
--- a/build.gradle.kts
+++ b/build.gradle.kts
@@ -4,7 +4,7 @@ plugins {
     id("ru.mipt.npm.project")
 }
 
-internal val kmathVersion: String by extra("0.2.0-dev-5")
+internal val kmathVersion: String by extra("0.2.0-dev-6")
 internal val bintrayRepo: String by extra("kscience")
 internal val githubProject: String by extra("kmath")
 
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/DotBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/DotBenchmark.kt
index 5c59afaee..4a0c96f78 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/DotBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/DotBenchmark.kt
@@ -3,14 +3,12 @@ package kscience.kmath.benchmarks
 import kotlinx.benchmark.Benchmark
 import kscience.kmath.commons.linear.CMMatrixContext
 import kscience.kmath.ejml.EjmlMatrixContext
-
 import kscience.kmath.linear.BufferMatrixContext
 import kscience.kmath.linear.RealMatrixContext
-import kscience.kmath.linear.real
+import kscience.kmath.nd.Matrix
 import kscience.kmath.operations.RealField
 import kscience.kmath.operations.invoke
 import kscience.kmath.structures.Buffer
-import kscience.kmath.structures.Matrix
 import org.openjdk.jmh.annotations.Scope
 import org.openjdk.jmh.annotations.State
 import kotlin.random.Random
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
index 5ff43ef80..9e5f76112 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
@@ -8,7 +8,6 @@ import kscience.kmath.commons.linear.inverse
 import kscience.kmath.ejml.EjmlMatrixContext
 import kscience.kmath.ejml.inverse
 import kscience.kmath.operations.invoke
-import kscience.kmath.structures.Matrix
 import org.openjdk.jmh.annotations.Scope
 import org.openjdk.jmh.annotations.State
 import kotlin.random.Random
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt
index 1be8e7236..3f103dd37 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt
@@ -1,5 +1,8 @@
 package kscience.kmath.benchmarks
 
+import kscience.kmath.nd.AbstractNDBuffer
+import kscience.kmath.nd.NDField
+import kscience.kmath.nd.RealNDField
 import kscience.kmath.operations.RealField
 import kscience.kmath.operations.invoke
 import kscience.kmath.structures.*
@@ -12,7 +15,7 @@ internal class NDFieldBenchmark {
     @Benchmark
     fun autoFieldAdd() {
         bufferedField {
-            var res: NDBuffer<Double> = one
+            var res: AbstractNDBuffer<Double> = one
             repeat(n) { res += one }
         }
     }
@@ -26,7 +29,7 @@ internal class NDFieldBenchmark {
     @Benchmark
     fun specializedFieldAdd() {
         specializedField {
-            var res: NDBuffer<Double> = one
+            var res: AbstractNDBuffer<Double> = one
             repeat(n) { res += 1.0 }
         }
     }
@@ -35,7 +38,7 @@ internal class NDFieldBenchmark {
     @Benchmark
     fun boxingFieldAdd() {
         genericField {
-            var res: NDBuffer<Double> = one
+            var res: AbstractNDBuffer<Double> = one
             repeat(n) { res += one }
         }
     }
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt
index 8663e353c..49b30cf7e 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt
@@ -1,10 +1,10 @@
 package kscience.kmath.benchmarks
 
+import kscience.kmath.nd.NDField
+import kscience.kmath.nd.RealNDField
 import kscience.kmath.operations.RealField
 import kscience.kmath.operations.invoke
 import kscience.kmath.structures.BufferedNDField
-import kscience.kmath.structures.NDField
-import kscience.kmath.structures.RealNDField
 import kscience.kmath.viktor.ViktorNDField
 import org.jetbrains.bio.viktor.F64Array
 import org.openjdk.jmh.annotations.Benchmark
diff --git a/examples/src/main/kotlin/kscience/kmath/operations/ComplexDemo.kt b/examples/src/main/kotlin/kscience/kmath/operations/ComplexDemo.kt
index e84fd8df3..f90221582 100644
--- a/examples/src/main/kotlin/kscience/kmath/operations/ComplexDemo.kt
+++ b/examples/src/main/kotlin/kscience/kmath/operations/ComplexDemo.kt
@@ -1,8 +1,7 @@
 package kscience.kmath.operations
 
+import kscience.kmath.nd.NDField
 import kscience.kmath.structures.NDElement
-import kscience.kmath.structures.NDField
-import kscience.kmath.structures.complex
 
 fun main() {
     // 2d element
diff --git a/examples/src/main/kotlin/kscience/kmath/structures/ComplexND.kt b/examples/src/main/kotlin/kscience/kmath/structures/ComplexND.kt
index b69590473..7dd725f0e 100644
--- a/examples/src/main/kotlin/kscience/kmath/structures/ComplexND.kt
+++ b/examples/src/main/kotlin/kscience/kmath/structures/ComplexND.kt
@@ -1,6 +1,7 @@
 package kscience.kmath.structures
 
 import kscience.kmath.linear.transpose
+import kscience.kmath.nd.*
 import kscience.kmath.operations.Complex
 import kscience.kmath.operations.ComplexField
 import kscience.kmath.operations.invoke
@@ -10,12 +11,12 @@ fun main() {
     val dim = 1000
     val n = 1000
 
-    val realField = NDField.real(dim, dim)
-    val complexField: ComplexNDField = NDField.complex(dim, dim)
+    val realField = NDAlgebra.real(dim, dim)
+    val complexField: ComplexNDField = NDAlgebra.complex(dim, dim)
 
     val realTime = measureTimeMillis {
         realField {
-            var res: NDBuffer<Double> = one
+            var res: NDStructure<Double> = one
             repeat(n) {
                 res += 1.0
             }
@@ -26,8 +27,10 @@ fun main() {
 
     val complexTime = measureTimeMillis {
         complexField {
-            var res: NDBuffer<Complex> = one
-            repeat(n) { res += 1.0 }
+            var res: NDStructure<Complex> = one
+            repeat(n) {
+                res += 1.0
+            }
         }
     }
 
diff --git a/examples/src/main/kotlin/kscience/kmath/structures/NDField.kt b/examples/src/main/kotlin/kscience/kmath/structures/NDField.kt
index b5130c92b..e26b8cbce 100644
--- a/examples/src/main/kotlin/kscience/kmath/structures/NDField.kt
+++ b/examples/src/main/kotlin/kscience/kmath/structures/NDField.kt
@@ -1,6 +1,7 @@
 package kscience.kmath.structures
 
 import kotlinx.coroutines.GlobalScope
+import kscience.kmath.nd.*
 import kscience.kmath.nd4j.Nd4jArrayField
 import kscience.kmath.operations.RealField
 import kscience.kmath.operations.invoke
@@ -22,29 +23,29 @@ fun main() {
     val n = 1000
 
     // automatically build context most suited for given type.
-    val autoField = NDField.auto(RealField, dim, dim)
+    val autoField = NDAlgebra.field(RealField, Buffer.Companion::auto, dim, dim)
     // specialized nd-field for Double. It works as generic Double field as well
-    val specializedField = NDField.real(dim, dim)
+    val specializedField = NDAlgebra.real(dim, dim)
     //A generic boxing field. It should be used for objects, not primitives.
-    val genericField = NDField.boxing(RealField, dim, dim)
+    val boxingField =  NDAlgebra.field(RealField, Buffer.Companion::boxing, dim, dim)
     // Nd4j specialized field.
     val nd4jField = Nd4jArrayField.real(dim, dim)
 
     measureAndPrint("Automatic field addition") {
         autoField {
-            var res: NDBuffer<Double> = one
+            var res: NDStructure<Double> = one
             repeat(n) { res += 1.0 }
         }
     }
 
     measureAndPrint("Element addition") {
-        var res = genericField.one
+        var res: NDStructure<Double> = boxingField.one
         repeat(n) { res += 1.0 }
     }
 
     measureAndPrint("Specialized addition") {
         specializedField {
-            var res: NDBuffer<Double> = one
+            var res: NDStructure<Double> = one
             repeat(n) { res += 1.0 }
         }
     }
@@ -70,8 +71,8 @@ fun main() {
 
     measureAndPrint("Generic addition") {
         //genericField.run(action)
-        genericField {
-            var res: NDBuffer<Double> = one
+        boxingField {
+            var res: NDStructure<Double> = one
             repeat(n) {
                 res += 1.0 // couldn't avoid using `one` due to resolution ambiguity }
             }
diff --git a/examples/src/main/kotlin/kscience/kmath/structures/StructureReadBenchmark.kt b/examples/src/main/kotlin/kscience/kmath/structures/StructureReadBenchmark.kt
index 51fd4f956..ba7bcb5f2 100644
--- a/examples/src/main/kotlin/kscience/kmath/structures/StructureReadBenchmark.kt
+++ b/examples/src/main/kotlin/kscience/kmath/structures/StructureReadBenchmark.kt
@@ -1,5 +1,7 @@
 package kscience.kmath.structures
 
+import kscience.kmath.nd.DefaultStrides
+import kscience.kmath.nd.NDBuffer
 import kotlin.system.measureTimeMillis
 
 fun main() {
@@ -7,7 +9,7 @@ fun main() {
     val array = DoubleArray(n * n) { 1.0 }
     val buffer = RealBuffer(array)
     val strides = DefaultStrides(intArrayOf(n, n))
-    val structure = BufferNDStructure(strides, buffer)
+    val structure = NDBuffer(strides, buffer)
 
     measureTimeMillis {
         var res = 0.0
diff --git a/examples/src/main/kotlin/kscience/kmath/structures/StructureWriteBenchmark.kt b/examples/src/main/kotlin/kscience/kmath/structures/StructureWriteBenchmark.kt
index db55b454f..307b068df 100644
--- a/examples/src/main/kotlin/kscience/kmath/structures/StructureWriteBenchmark.kt
+++ b/examples/src/main/kotlin/kscience/kmath/structures/StructureWriteBenchmark.kt
@@ -1,5 +1,7 @@
 package kscience.kmath.structures
 
+import kscience.kmath.nd.NDStructure
+import kscience.kmath.nd.mapToBuffer
 import kotlin.system.measureTimeMillis
 
 fun main() {
diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
index 850446afa..b9fc8b72a 100644
--- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
+++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMMatrix.kt
@@ -3,9 +3,8 @@ package kscience.kmath.commons.linear
 import kscience.kmath.linear.DiagonalFeature
 import kscience.kmath.linear.MatrixContext
 import kscience.kmath.linear.Point
-import kscience.kmath.linear.origin
 import kscience.kmath.misc.UnstableKMathAPI
-import kscience.kmath.structures.Matrix
+import kscience.kmath.nd.Matrix
 import org.apache.commons.math3.linear.*
 import kotlin.reflect.KClass
 import kotlin.reflect.cast
diff --git a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMSolver.kt b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMSolver.kt
index 210014e1a..aa0bf4e1a 100644
--- a/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMSolver.kt
+++ b/kmath-commons/src/main/kotlin/kscience/kmath/commons/linear/CMSolver.kt
@@ -1,7 +1,7 @@
 package kscience.kmath.commons.linear
 
 import kscience.kmath.linear.Point
-import kscience.kmath.structures.Matrix
+import kscience.kmath.nd.Matrix
 import org.apache.commons.math3.linear.*
 
 public enum class CMDecomposition {
diff --git a/kmath-core/README.md b/kmath-core/README.md
index 7882e5252..9ed54b9eb 100644
--- a/kmath-core/README.md
+++ b/kmath-core/README.md
@@ -3,7 +3,7 @@
 The core features of KMath:
 
  - [algebras](src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt) : Algebraic structures: contexts and elements
- - [nd](src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt) : Many-dimensional structures
+ - [nd](src/commonMain/kotlin/kscience/kmath/nd/NDStructure.kt) : Many-dimensional structures
  - [buffers](src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : One-dimensional structure
  - [expressions](src/commonMain/kotlin/kscience/kmath/expressions) : Functional Expressions
  - [domains](src/commonMain/kotlin/kscience/kmath/domains) : Domains
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/SymbolIndexer.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/SymbolIndexer.kt
index 6c61c7c7d..fd0661080 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/SymbolIndexer.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/SymbolIndexer.kt
@@ -1,8 +1,8 @@
 package kscience.kmath.expressions
 
 import kscience.kmath.linear.Point
+import kscience.kmath.nd.Structure2D
 import kscience.kmath.structures.BufferFactory
-import kscience.kmath.structures.Structure2D
 
 /**
  * An environment to easy transform indexed variables to symbols and back.
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
index a74d948fc..bcd1c9976 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/BufferMatrix.kt
@@ -1,7 +1,18 @@
 package kscience.kmath.linear
 
+import kscience.kmath.nd.NDStructure
+import kscience.kmath.nd.Structure2D
 import kscience.kmath.operations.Ring
-import kscience.kmath.structures.*
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.BufferFactory
+import kscience.kmath.structures.asSequence
+
+/**
+ * Alias for [Structure2D] with more familiar name.
+ *
+ * @param T the type of items.
+ */
+public typealias Matrix<T> = Structure2D<T>
 
 /**
  * Basic implementation of Matrix space based on [NDStructure]
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LinearAlgebra.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LinearAlgebra.kt
index 034decc2f..d9201d0f1 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LinearAlgebra.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LinearAlgebra.kt
@@ -1,7 +1,6 @@
 package kscience.kmath.linear
 
 import kscience.kmath.structures.Buffer
-import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.VirtualBuffer
 
 public typealias Point<T> = Buffer<T>
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LupDecomposition.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LupDecomposition.kt
index 5cf7c8f70..ddf7c3fd2 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LupDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/LupDecomposition.kt
@@ -1,6 +1,7 @@
 package kscience.kmath.linear
 
 import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.nd.getFeature
 import kscience.kmath.operations.*
 import kscience.kmath.structures.*
 
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixBuilder.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixBuilder.kt
index c0c209248..cb01939cc 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixBuilder.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixBuilder.kt
@@ -1,6 +1,9 @@
 package kscience.kmath.linear
 
-import kscience.kmath.structures.*
+import kscience.kmath.nd.Structure2D
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.BufferFactory
+import kscience.kmath.structures.asBuffer
 
 public class MatrixBuilder(public val rows: Int, public val columns: Int) {
     public operator fun <T : Any> invoke(vararg elements: T): Matrix<T> {
@@ -22,7 +25,7 @@ public fun <T : Any> Structure2D.Companion.row(vararg values: T): Matrix<T> {
 public inline fun <reified T : Any> Structure2D.Companion.row(
     size: Int,
     factory: BufferFactory<T> = Buffer.Companion::auto,
-    noinline builder: (Int) -> T
+    noinline builder: (Int) -> T,
 ): Matrix<T> {
     val buffer = factory(size, builder)
     return BufferMatrix(1, size, buffer)
@@ -36,7 +39,7 @@ public fun <T : Any> Structure2D.Companion.column(vararg values: T): Matrix<T> {
 public inline fun <reified T : Any> Structure2D.Companion.column(
     size: Int,
     factory: BufferFactory<T> = Buffer.Companion::auto,
-    noinline builder: (Int) -> T
+    noinline builder: (Int) -> T,
 ): Matrix<T> {
     val buffer = factory(size, builder)
     return BufferMatrix(size, 1, buffer)
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
index 59a41f840..8766f2e5f 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixContext.kt
@@ -6,7 +6,6 @@ import kscience.kmath.operations.invoke
 import kscience.kmath.operations.sum
 import kscience.kmath.structures.Buffer
 import kscience.kmath.structures.BufferFactory
-import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.asSequence
 
 /**
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixFeatures.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixFeatures.kt
index e61feec6c..250e0cc04 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixFeatures.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixFeatures.kt
@@ -1,7 +1,5 @@
 package kscience.kmath.linear
 
-import kscience.kmath.structures.Matrix
-
 /**
  * A marker interface representing some properties of matrices or additional transformations of them. Features are used
  * to optimize matrix operations performance in some cases or retrieve the APIs.
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixWrapper.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixWrapper.kt
index 362db1fe7..726a91e97 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixWrapper.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/MatrixWrapper.kt
@@ -1,11 +1,10 @@
 package kscience.kmath.linear
 
 import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.nd.Structure2D
+import kscience.kmath.nd.getFeature
 import kscience.kmath.operations.Ring
-import kscience.kmath.structures.Matrix
-import kscience.kmath.structures.Structure2D
 import kscience.kmath.structures.asBuffer
-import kscience.kmath.structures.getFeature
 import kotlin.math.sqrt
 import kotlin.reflect.KClass
 import kotlin.reflect.safeCast
@@ -39,7 +38,8 @@ public class MatrixWrapper<T : Any> internal  constructor(
  * Origin does not necessary store all features.
  */
 @UnstableKMathAPI
-public val <T : Any> Matrix<T>.origin: Matrix<T> get() = (this as? MatrixWrapper)?.origin ?: this
+public val <T : Any> Matrix<T>.origin: Matrix<T>
+    get() = (this as? MatrixWrapper)?.origin ?: this
 
 /**
  * Add a single feature to a [Matrix]
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt
index 8e197672f..3d544a9af 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/RealMatrixContext.kt
@@ -1,6 +1,5 @@
 package kscience.kmath.linear
 
-import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.RealBuffer
 
 @Suppress("OVERRIDE_BY_INLINE")
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/VirtualMatrix.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/VirtualMatrix.kt
index 0269a64d1..df4e4de09 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/VirtualMatrix.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/linear/VirtualMatrix.kt
@@ -1,7 +1,5 @@
 package kscience.kmath.linear
 
-import kscience.kmath.structures.Matrix
-
 public class VirtualMatrix<T : Any>(
     override val rowNum: Int,
     override val colNum: Int,
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/BufferNDAlgebra.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/BufferNDAlgebra.kt
new file mode 100644
index 000000000..a92e08c60
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/BufferNDAlgebra.kt
@@ -0,0 +1,128 @@
+package kscience.kmath.nd
+
+import kscience.kmath.nd.*
+import kscience.kmath.operations.*
+import kscience.kmath.structures.BufferFactory
+import kotlin.contracts.InvocationKind
+import kotlin.contracts.contract
+
+public interface BufferNDAlgebra<T, C> : NDAlgebra<T, C> {
+    public val strides: Strides
+    public val bufferFactory: BufferFactory<T>
+
+    override fun produce(initializer: C.(IntArray) -> T): NDBuffer<T> = NDBuffer(
+        strides,
+        bufferFactory(strides.linearSize) { offset ->
+            elementContext.initializer(strides.index(offset))
+        }
+    )
+
+    public val NDStructure<T>.ndBuffer: NDBuffer<T>
+        get() = when {
+            !shape.contentEquals(this@BufferNDAlgebra.shape) -> throw ShapeMismatchException(
+                this@BufferNDAlgebra.shape,
+                shape
+            )
+            this is NDBuffer && this.strides == this@BufferNDAlgebra.strides -> this
+            else -> produce { this@ndBuffer[it] }
+        }
+
+    override fun map(arg: NDStructure<T>, transform: C.(T) -> T): NDBuffer<T> {
+        val argAsBuffer = arg.ndBuffer
+        val buffer = bufferFactory(strides.linearSize) { offset ->
+            elementContext.transform(argAsBuffer.buffer[offset])
+        }
+        return NDBuffer(strides, buffer)
+    }
+
+    override fun mapIndexed(arg: NDStructure<T>, transform: C.(index: IntArray, T) -> T): NDStructure<T> {
+        val argAsBuffer = arg.ndBuffer
+        val buffer = bufferFactory(strides.linearSize) { offset ->
+            elementContext.transform(
+                strides.index(offset),
+                argAsBuffer[offset]
+            )
+        }
+        return NDBuffer(strides, buffer)
+    }
+
+    override fun combine(a: NDStructure<T>, b: NDStructure<T>, transform: C.(T, T) -> T): NDStructure<T> {
+        val aBuffer = a.ndBuffer
+        val bBuffer = b.ndBuffer
+        val buffer = bufferFactory(strides.linearSize) { offset ->
+            elementContext.transform(aBuffer.buffer[offset], bBuffer[offset])
+        }
+        return NDBuffer(strides, buffer)
+    }
+}
+
+public open class BufferedNDSpace<T, R : Space<T>>(
+    final override val shape: IntArray,
+    final override val elementContext: R,
+    final override val bufferFactory: BufferFactory<T>,
+) : NDSpace<T, R>, BufferNDAlgebra<T, R> {
+    override val strides: Strides = DefaultStrides(shape)
+    override val zero: NDBuffer<T> by lazy { produce { zero } }
+}
+
+public open class BufferedNDRing<T, R : Ring<T>>(
+    shape: IntArray,
+    elementContext: R,
+    bufferFactory: BufferFactory<T>,
+) : BufferedNDSpace<T, R>(shape, elementContext, bufferFactory), NDRing<T, R> {
+    override val one: NDBuffer<T> by lazy { produce { one } }
+}
+
+public open class BufferedNDField<T, R : Field<T>>(
+    shape: IntArray,
+    elementContext: R,
+    bufferFactory: BufferFactory<T>,
+) : BufferedNDRing<T, R>(shape, elementContext, bufferFactory), NDField<T, R>
+
+// space factories
+public fun <T, A : Space<T>> NDAlgebra.Companion.space(
+    space: A,
+    bufferFactory: BufferFactory<T>,
+    vararg shape: Int,
+): BufferedNDSpace<T, A> = BufferedNDSpace(shape, space, bufferFactory)
+
+public inline fun <T, A : Space<T>, R> A.nd(
+    noinline bufferFactory: BufferFactory<T>,
+    vararg shape: Int,
+    action: BufferedNDSpace<T, A>.() -> R,
+): R {
+    contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
+    return NDAlgebra.space(this, bufferFactory, *shape).run(action)
+}
+
+//ring factories
+public fun <T, A : Ring<T>> NDAlgebra.Companion.ring(
+    ring: A,
+    bufferFactory: BufferFactory<T>,
+    vararg shape: Int,
+): BufferedNDRing<T, A> = BufferedNDRing(shape, ring, bufferFactory)
+
+public inline fun <T, A : Ring<T>, R> A.nd(
+    noinline bufferFactory: BufferFactory<T>,
+    vararg shape: Int,
+    action: BufferedNDRing<T, A>.() -> R,
+): R {
+    contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
+    return NDAlgebra.ring(this, bufferFactory, *shape).run(action)
+}
+
+//field factories
+public fun <T, A : Field<T>> NDAlgebra.Companion.field(
+    field: A,
+    bufferFactory: BufferFactory<T>,
+    vararg shape: Int,
+): BufferedNDField<T, A> = BufferedNDField(shape, field, bufferFactory)
+
+public inline fun <T, A : Field<T>, R> A.nd(
+    noinline bufferFactory: BufferFactory<T>,
+    vararg shape: Int,
+    action: BufferedNDField<T, A>.() -> R,
+): R {
+    contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
+    return NDAlgebra.field(this, bufferFactory, *shape).run(action)
+}
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/ComplexNDField.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/ComplexNDField.kt
new file mode 100644
index 000000000..dbcce5440
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/ComplexNDField.kt
@@ -0,0 +1,113 @@
+package kscience.kmath.nd
+
+import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.operations.*
+import kscience.kmath.structures.Buffer
+import kotlin.contracts.InvocationKind
+import kotlin.contracts.contract
+
+
+/**
+ * An optimized nd-field for complex numbers
+ */
+@OptIn(UnstableKMathAPI::class)
+public open class ComplexNDField(
+    shape: IntArray,
+) : BufferedNDField<Complex, ComplexField>(shape, ComplexField, Buffer.Companion::complex),
+    RingWithNumbers<NDStructure<Complex>>,
+    ExtendedField<NDStructure<Complex>> {
+
+    override val zero: NDBuffer<Complex> by lazy { produce { zero } }
+    override val one: NDBuffer<Complex> by lazy { produce { one } }
+
+    override fun number(value: Number): NDBuffer<Complex> {
+        val d = value.toComplex() // minimize conversions
+        return produce { d }
+    }
+//
+//    @Suppress("OVERRIDE_BY_INLINE")
+//    override inline fun map(
+//        arg: AbstractNDBuffer<Double>,
+//        transform: RealField.(Double) -> Double,
+//    ): RealNDElement {
+//        check(arg)
+//        val array = RealBuffer(arg.strides.linearSize) { offset -> RealField.transform(arg.buffer[offset]) }
+//        return BufferedNDFieldElement(this, array)
+//    }
+//
+//    @Suppress("OVERRIDE_BY_INLINE")
+//    override inline fun produce(initializer: RealField.(IntArray) -> Double): RealNDElement {
+//        val array = RealBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) }
+//        return BufferedNDFieldElement(this, array)
+//    }
+//
+//    @Suppress("OVERRIDE_BY_INLINE")
+//    override inline fun mapIndexed(
+//        arg: AbstractNDBuffer<Double>,
+//        transform: RealField.(index: IntArray, Double) -> Double,
+//    ): RealNDElement {
+//        check(arg)
+//        return BufferedNDFieldElement(
+//            this,
+//            RealBuffer(arg.strides.linearSize) { offset ->
+//                elementContext.transform(
+//                    arg.strides.index(offset),
+//                    arg.buffer[offset]
+//                )
+//            })
+//    }
+//
+//    @Suppress("OVERRIDE_BY_INLINE")
+//    override inline fun combine(
+//        a: AbstractNDBuffer<Double>,
+//        b: AbstractNDBuffer<Double>,
+//        transform: RealField.(Double, Double) -> Double,
+//    ): RealNDElement {
+//        check(a, b)
+//        val buffer = RealBuffer(strides.linearSize) { offset ->
+//            elementContext.transform(a.buffer[offset], b.buffer[offset])
+//        }
+//        return BufferedNDFieldElement(this, buffer)
+//    }
+
+    override fun power(arg: NDStructure<Complex>, pow: Number): NDBuffer<Complex> = map(arg) { power(it, pow) }
+
+    override fun exp(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { exp(it) }
+
+    override fun ln(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { ln(it) }
+
+    override fun sin(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { sin(it) }
+    override fun cos(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { cos(it) }
+    override fun tan(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { tan(it) }
+    override fun asin(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { asin(it) }
+    override fun acos(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { acos(it) }
+    override fun atan(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { atan(it) }
+
+    override fun sinh(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { sinh(it) }
+    override fun cosh(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { cosh(it) }
+    override fun tanh(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { tanh(it) }
+    override fun asinh(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { asinh(it) }
+    override fun acosh(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { acosh(it) }
+    override fun atanh(arg: NDStructure<Complex>): NDBuffer<Complex> = map(arg) { atanh(it) }
+}
+
+
+/**
+ * Fast element production using function inlining
+ */
+public inline fun BufferedNDField<Complex, ComplexField>.produceInline(initializer: ComplexField.(Int) -> Complex): NDBuffer<Complex> {
+    contract { callsInPlace(initializer, InvocationKind.EXACTLY_ONCE) }
+    val buffer = Buffer.complex(strides.linearSize) { offset -> ComplexField.initializer(offset) }
+    return NDBuffer(strides, buffer)
+}
+
+
+public fun NDAlgebra.Companion.complex(vararg shape: Int): ComplexNDField = ComplexNDField(shape)
+
+/**
+ * Produce a context for n-dimensional operations inside this real field
+ */
+public inline fun <R> ComplexField.nd(vararg shape: Int, action: ComplexNDField.() -> R): R {
+    contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
+    return ComplexNDField(shape).action()
+}
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDAlgebra.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDAlgebra.kt
new file mode 100644
index 000000000..facc7eb38
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDAlgebra.kt
@@ -0,0 +1,262 @@
+package kscience.kmath.nd
+
+import kscience.kmath.operations.Field
+import kscience.kmath.operations.Ring
+import kscience.kmath.operations.Space
+import kscience.kmath.structures.*
+
+/**
+ * An exception is thrown when the expected ans actual shape of NDArray differs.
+ *
+ * @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()}.")
+
+/**
+ * The base interface for all ND-algebra implementations.
+ *
+ * @param T the type of ND-structure element.
+ * @param C the type of the element context.
+ * @param N the type of the structure.
+ */
+public interface NDAlgebra<T, C> {
+    /**
+     * The shape of ND-structures this algebra operates on.
+     */
+    public val shape: IntArray
+
+    /**
+     * The algebra over elements of ND structure.
+     */
+    public val elementContext: C
+
+    /**
+     * Produces a new [N] structure using given initializer function.
+     */
+    public fun produce(initializer: C.(IntArray) -> T): NDStructure<T>
+
+    /**
+     * Maps elements from one structure to another one by applying [transform] to them.
+     */
+    public fun map(arg: NDStructure<T>, transform: C.(T) -> T): NDStructure<T>
+
+    /**
+     * Maps elements from one structure to another one by applying [transform] to them alongside with their indices.
+     */
+    public fun mapIndexed(arg: NDStructure<T>, transform: C.(index: IntArray, T) -> T): NDStructure<T>
+
+    /**
+     * Combines two structures into one.
+     */
+    public fun combine(a: NDStructure<T>, b: NDStructure<T>, transform: C.(T, T) -> T): NDStructure<T>
+
+    /**
+     * Element-wise invocation of function working on [T] on a [NDStructure].
+     */
+    public operator fun Function1<T, T>.invoke(structure: NDStructure<T>): NDStructure<T> =
+        map(structure) { value -> this@invoke(value) }
+
+    public companion object
+}
+
+/**
+ * Checks if given elements are consistent with this context.
+ *
+ * @param structures the structures to check.
+ * @return the array of valid structures.
+ */
+internal fun <T, C> NDAlgebra<T, C>.checkShape(vararg structures: NDStructure<T>): Array<out NDStructure<T>> = structures
+    .map(NDStructure<T>::shape)
+    .singleOrNull { !shape.contentEquals(it) }
+    ?.let<IntArray, Array<out NDStructure<T>>> { throw ShapeMismatchException(shape, it) }
+    ?: structures
+
+/**
+ * Checks if given element is consistent with this context.
+ *
+ * @param element the structure to check.
+ * @return the valid structure.
+ */
+internal fun <T, C> NDAlgebra<T, C>.checkShape(element: NDStructure<T>): NDStructure<T> {
+    if (!element.shape.contentEquals(shape)) throw ShapeMismatchException(shape, element.shape)
+    return element
+}
+
+/**
+ * Space of [NDStructure].
+ *
+ * @param T the type of the element contained in ND structure.
+ * @param N the type of ND structure.
+ * @param S the type of space of structure elements.
+ */
+public interface NDSpace<T, S : Space<T>> : Space<NDStructure<T>>, NDAlgebra<T, S> {
+    /**
+     * Element-wise addition.
+     *
+     * @param a the addend.
+     * @param b the augend.
+     * @return the sum.
+     */
+    public override fun add(a: NDStructure<T>, b: NDStructure<T>): NDStructure<T> =
+        combine(a, b) { aValue, bValue -> add(aValue, bValue) }
+
+    /**
+     * Element-wise multiplication by scalar.
+     *
+     * @param a the multiplicand.
+     * @param k the multiplier.
+     * @return the product.
+     */
+    public override fun multiply(a: NDStructure<T>, k: Number): NDStructure<T> = map(a) { multiply(it, k) }
+
+    // TODO move to extensions after KEEP-176
+
+    /**
+     * Adds an ND structure to an element of it.
+     *
+     * @receiver the addend.
+     * @param arg the augend.
+     * @return the sum.
+     */
+    public operator fun NDStructure<T>.plus(arg: T): NDStructure<T> = map(this) { value -> add(arg, value) }
+
+    /**
+     * Subtracts an element from ND structure of it.
+     *
+     * @receiver the dividend.
+     * @param arg the divisor.
+     * @return the quotient.
+     */
+    public operator fun NDStructure<T>.minus(arg: T): NDStructure<T> = map(this) { value -> add(arg, -value) }
+
+    /**
+     * Adds an element to ND structure of it.
+     *
+     * @receiver the addend.
+     * @param arg the augend.
+     * @return the sum.
+     */
+    public operator fun T.plus(arg: NDStructure<T>): NDStructure<T> = map(arg) { value -> add(this@plus, value) }
+
+    /**
+     * Subtracts an ND structure from an element of it.
+     *
+     * @receiver the dividend.
+     * @param arg the divisor.
+     * @return the quotient.
+     */
+    public operator fun T.minus(arg: NDStructure<T>): NDStructure<T> = map(arg) { value -> add(-this@minus, value) }
+
+    public companion object
+}
+
+/**
+ * Ring of [NDStructure].
+ *
+ * @param T the type of the element contained in ND structure.
+ * @param N the type of ND structure.
+ * @param R the type of ring of structure elements.
+ */
+public interface NDRing<T, R : Ring<T>> : Ring<NDStructure<T>>, NDSpace<T, R> {
+    /**
+     * Element-wise multiplication.
+     *
+     * @param a the multiplicand.
+     * @param b the multiplier.
+     * @return the product.
+     */
+    public override fun multiply(a: NDStructure<T>, b: NDStructure<T>): NDStructure<T> =
+        combine(a, b) { aValue, bValue -> multiply(aValue, bValue) }
+
+    //TODO move to extensions after KEEP-176
+
+    /**
+     * Multiplies an ND structure by an element of it.
+     *
+     * @receiver the multiplicand.
+     * @param arg the multiplier.
+     * @return the product.
+     */
+    public operator fun NDStructure<T>.times(arg: T): NDStructure<T> = map(this) { value -> multiply(arg, value) }
+
+    /**
+     * Multiplies an element by a ND structure of it.
+     *
+     * @receiver the multiplicand.
+     * @param arg the multiplier.
+     * @return the product.
+     */
+    public operator fun T.times(arg: NDStructure<T>): NDStructure<T> = map(arg) { value -> multiply(this@times, value) }
+
+    public companion object
+}
+
+/**
+ * Field of [NDStructure].
+ *
+ * @param T the type of the element contained in ND structure.
+ * @param N the type of ND structure.
+ * @param F the type field of structure elements.
+ */
+public interface NDField<T, F : Field<T>> : Field<NDStructure<T>>, NDRing<T, F> {
+    /**
+     * Element-wise division.
+     *
+     * @param a the dividend.
+     * @param b the divisor.
+     * @return the quotient.
+     */
+    public override fun divide(a: NDStructure<T>, b: NDStructure<T>): NDStructure<T> =
+        combine(a, b) { aValue, bValue -> divide(aValue, bValue) }
+
+    //TODO move to extensions after KEEP-176
+    /**
+     * Divides an ND structure by an element of it.
+     *
+     * @receiver the dividend.
+     * @param arg the divisor.
+     * @return the quotient.
+     */
+    public operator fun NDStructure<T>.div(arg: T): NDStructure<T> = map(this) { value -> divide(arg, value) }
+
+    /**
+     * Divides an element by an ND structure of it.
+     *
+     * @receiver the dividend.
+     * @param arg the divisor.
+     * @return the quotient.
+     */
+    public operator fun T.div(arg: NDStructure<T>): NDStructure<T> = map(arg) { divide(it, this@div) }
+
+//    @ThreadLocal
+//    public companion object {
+//        private val realNDFieldCache: MutableMap<IntArray, RealNDField> = hashMapOf()
+//
+//        /**
+//         * Create a nd-field for [Double] values or pull it from cache if it was created previously.
+//         */
+//        public fun real(vararg shape: Int): RealNDField = realNDFieldCache.getOrPut(shape) { RealNDField(shape) }
+//
+//        /**
+//         * Create an ND field with boxing generic buffer.
+//         */
+//        public fun <T : Any, F : Field<T>> boxing(
+//            field: F,
+//            vararg shape: Int,
+//            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
+//        ): BufferedNDField<T, F> = BufferedNDField(shape, field, bufferFactory)
+//
+//        /**
+//         * Create a most suitable implementation for nd-field using reified class.
+//         */
+//        @Suppress("UNCHECKED_CAST")
+//        public inline fun <reified T : Any, F : Field<T>> auto(field: F, vararg shape: Int): NDField<T, F> =
+//            when {
+//                T::class == Double::class -> real(*shape) as NDField<T, F>
+//                T::class == Complex::class -> complex(*shape) as BufferedNDField<T, F>
+//                else -> BoxingNDField(shape, field, Buffer.Companion::auto)
+//            }
+//    }
+}
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDStructure.kt
similarity index 84%
rename from kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt
rename to kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDStructure.kt
index e7d89ca7e..09c6fa198 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDStructure.kt
@@ -1,6 +1,10 @@
-package kscience.kmath.structures
+package kscience.kmath.nd
 
 import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.BufferFactory
+import kscience.kmath.structures.MutableBuffer
+import kscience.kmath.structures.asSequence
 import kotlin.jvm.JvmName
 import kotlin.native.concurrent.ThreadLocal
 import kotlin.reflect.KClass
@@ -74,8 +78,8 @@ public interface NDStructure<T> {
             strides: Strides,
             bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
             initializer: (IntArray) -> T,
-        ): BufferNDStructure<T> =
-            BufferNDStructure(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
+        ): NDBuffer<T> =
+            NDBuffer(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
 
         /**
          * Inline create NDStructure with non-boxing buffer implementation if it is possible
@@ -83,40 +87,40 @@ public interface NDStructure<T> {
         public inline fun <reified T : Any> auto(
             strides: Strides,
             crossinline initializer: (IntArray) -> T,
-        ): BufferNDStructure<T> =
-            BufferNDStructure(strides, Buffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) })
+        ): NDBuffer<T> =
+            NDBuffer(strides, Buffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) })
 
         public inline fun <T : Any> auto(
             type: KClass<T>,
             strides: Strides,
             crossinline initializer: (IntArray) -> T,
-        ): BufferNDStructure<T> =
-            BufferNDStructure(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) })
+        ): NDBuffer<T> =
+            NDBuffer(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) })
 
         public fun <T> build(
             shape: IntArray,
             bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
             initializer: (IntArray) -> T,
-        ): BufferNDStructure<T> = build(DefaultStrides(shape), bufferFactory, initializer)
+        ): NDBuffer<T> = build(DefaultStrides(shape), bufferFactory, initializer)
 
         public inline fun <reified T : Any> auto(
             shape: IntArray,
             crossinline initializer: (IntArray) -> T,
-        ): BufferNDStructure<T> =
+        ): NDBuffer<T> =
             auto(DefaultStrides(shape), initializer)
 
         @JvmName("autoVarArg")
         public inline fun <reified T : Any> auto(
             vararg shape: Int,
             crossinline initializer: (IntArray) -> T,
-        ): BufferNDStructure<T> =
+        ): NDBuffer<T> =
             auto(DefaultStrides(shape), initializer)
 
         public inline fun <T : Any> auto(
             type: KClass<T>,
             vararg shape: Int,
             crossinline initializer: (IntArray) -> T,
-        ): BufferNDStructure<T> =
+        ): NDBuffer<T> =
             auto(type, DefaultStrides(shape), initializer)
     }
 }
@@ -156,7 +160,7 @@ public inline fun <T> MutableNDStructure<T>.mapInPlace(action: (IntArray, T) ->
  */
 public interface Strides {
     /**
-     * Shape of NDstructure
+     * Shape of NDStructure
      */
     public val shape: IntArray
 
@@ -185,7 +189,9 @@ public interface Strides {
     /**
      * Iterate over ND indices in a natural order
      */
-    public fun indices(): Sequence<IntArray> = (0 until linearSize).asSequence().map { index(it) }
+    public fun indices(): Sequence<IntArray> = (0 until linearSize).asSequence().map {
+        index(it)
+    }
 }
 
 /**
@@ -211,7 +217,7 @@ public class DefaultStrides private constructor(override val shape: IntArray) :
     }
 
     override fun offset(index: IntArray): Int = index.mapIndexed { i, value ->
-        if (value < 0 || value >= this.shape[i])
+        if (value < 0 || value >= shape[i])
             throw IndexOutOfBoundsException("Index $value out of shape bounds: (0,${this.shape[i]})")
 
         value * strides[i]
@@ -256,23 +262,29 @@ public class DefaultStrides private constructor(override val shape: IntArray) :
  * Represents [NDStructure] over [Buffer].
  *
  * @param T the type of items.
+ * @param strides The strides to access elements of [Buffer] by linear indices.
+ * @param buffer The underlying buffer.
  */
-public abstract class NDBuffer<T> : NDStructure<T> {
-    /**
-     * The underlying buffer.
-     */
-    public abstract val buffer: Buffer<T>
+public open class NDBuffer<T>(
+    public val strides: Strides,
+    buffer: Buffer<T>,
+) : NDStructure<T> {
 
-    /**
-     * The strides to access elements of [Buffer] by linear indices.
-     */
-    public abstract val strides: Strides
+    init {
+        if (strides.linearSize != buffer.size) {
+            error("Expected buffer side of ${strides.linearSize}, but found ${buffer.size}")
+        }
+    }
+
+    public open val buffer: Buffer<T> = buffer
 
     override operator fun get(index: IntArray): T = buffer[strides.offset(index)]
 
     override val shape: IntArray get() = strides.shape
 
-    override fun elements(): Sequence<Pair<IntArray, T>> = strides.indices().map { it to this[it] }
+    override fun elements(): Sequence<Pair<IntArray, T>> = strides.indices().map {
+        it to this[it]
+    }
 
     override fun equals(other: Any?): Boolean {
         return NDStructure.contentEquals(this, other as? NDStructure<*> ?: return false)
@@ -297,46 +309,30 @@ public abstract class NDBuffer<T> : NDStructure<T> {
         }
         return "NDBuffer(shape=${shape.contentToString()}, buffer=$bufferRepr)"
     }
-
-
 }
 
 /**
- * Boxing generic [NDStructure]
- */
-public class BufferNDStructure<T>(
-    override val strides: Strides,
-    override val buffer: Buffer<T>,
-) : NDBuffer<T>() {
-    init {
-        if (strides.linearSize != buffer.size) {
-            error("Expected buffer side of ${strides.linearSize}, but found ${buffer.size}")
-        }
-    }
-}
-
-/**
- * Transform structure to a new structure using provided [BufferFactory] and optimizing if argument is [BufferNDStructure]
+ * Transform structure to a new structure using provided [BufferFactory] and optimizing if argument is [NDBuffer]
  */
 public inline fun <T, reified R : Any> NDStructure<T>.mapToBuffer(
     factory: BufferFactory<R> = Buffer.Companion::auto,
     crossinline transform: (T) -> R,
-): BufferNDStructure<R> {
-    return if (this is BufferNDStructure<T>)
-        BufferNDStructure(this.strides, factory.invoke(strides.linearSize) { transform(buffer[it]) })
+): NDBuffer<R> {
+    return if (this is NDBuffer<T>)
+        NDBuffer(this.strides, factory.invoke(strides.linearSize) { transform(buffer[it]) })
     else {
         val strides = DefaultStrides(shape)
-        BufferNDStructure(strides, factory.invoke(strides.linearSize) { transform(get(strides.index(it))) })
+        NDBuffer(strides, factory.invoke(strides.linearSize) { transform(get(strides.index(it))) })
     }
 }
 
 /**
  * Mutable ND buffer based on linear [MutableBuffer].
  */
-public class MutableBufferNDStructure<T>(
-    override val strides: Strides,
-    override val buffer: MutableBuffer<T>,
-) : NDBuffer<T>(), MutableNDStructure<T> {
+public class MutableNDBuffer<T>(
+    strides: Strides,
+    buffer: MutableBuffer<T>,
+) : NDBuffer<T>(strides,buffer), MutableNDStructure<T> {
 
     init {
         require(strides.linearSize == buffer.size) {
@@ -344,6 +340,8 @@ public class MutableBufferNDStructure<T>(
         }
     }
 
+    override val buffer: MutableBuffer<T> = super.buffer as MutableBuffer<T>
+
     override operator fun set(index: IntArray, value: T): Unit = buffer.set(strides.offset(index), value)
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/RealNDField.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/RealNDField.kt
new file mode 100644
index 000000000..feb3c509a
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/RealNDField.kt
@@ -0,0 +1,111 @@
+package kscience.kmath.nd
+
+import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.operations.ExtendedField
+import kscience.kmath.operations.RealField
+import kscience.kmath.operations.RingWithNumbers
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.RealBuffer
+import kotlin.contracts.InvocationKind
+import kotlin.contracts.contract
+
+@OptIn(UnstableKMathAPI::class)
+public open class RealNDField(
+    shape: IntArray,
+) : BufferedNDField<Double, RealField>(shape, RealField, Buffer.Companion::real),
+    RingWithNumbers<NDStructure<Double>>,
+    ExtendedField<NDStructure<Double>> {
+
+    override val zero: NDBuffer<Double> by lazy { produce { zero } }
+    override val one: NDBuffer<Double> by lazy { produce { one } }
+
+    override fun number(value: Number): NDBuffer<Double> {
+        val d = value.toDouble() // minimize conversions
+        return produce { d }
+    }
+//
+//    @Suppress("OVERRIDE_BY_INLINE")
+//    override inline fun map(
+//        arg: AbstractNDBuffer<Double>,
+//        transform: RealField.(Double) -> Double,
+//    ): RealNDElement {
+//        check(arg)
+//        val array = RealBuffer(arg.strides.linearSize) { offset -> RealField.transform(arg.buffer[offset]) }
+//        return BufferedNDFieldElement(this, array)
+//    }
+//
+//    @Suppress("OVERRIDE_BY_INLINE")
+//    override inline fun produce(initializer: RealField.(IntArray) -> Double): RealNDElement {
+//        val array = RealBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) }
+//        return BufferedNDFieldElement(this, array)
+//    }
+//
+//    @Suppress("OVERRIDE_BY_INLINE")
+//    override inline fun mapIndexed(
+//        arg: AbstractNDBuffer<Double>,
+//        transform: RealField.(index: IntArray, Double) -> Double,
+//    ): RealNDElement {
+//        check(arg)
+//        return BufferedNDFieldElement(
+//            this,
+//            RealBuffer(arg.strides.linearSize) { offset ->
+//                elementContext.transform(
+//                    arg.strides.index(offset),
+//                    arg.buffer[offset]
+//                )
+//            })
+//    }
+//
+//    @Suppress("OVERRIDE_BY_INLINE")
+//    override inline fun combine(
+//        a: AbstractNDBuffer<Double>,
+//        b: AbstractNDBuffer<Double>,
+//        transform: RealField.(Double, Double) -> Double,
+//    ): RealNDElement {
+//        check(a, b)
+//        val buffer = RealBuffer(strides.linearSize) { offset ->
+//            elementContext.transform(a.buffer[offset], b.buffer[offset])
+//        }
+//        return BufferedNDFieldElement(this, buffer)
+//    }
+
+    override fun power(arg: NDStructure<Double>, pow: Number): NDBuffer<Double> = map(arg) { power(it, pow) }
+
+    override fun exp(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { exp(it) }
+
+    override fun ln(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { ln(it) }
+
+    override fun sin(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { sin(it) }
+    override fun cos(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { cos(it) }
+    override fun tan(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { tan(it) }
+    override fun asin(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { asin(it) }
+    override fun acos(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { acos(it) }
+    override fun atan(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { atan(it) }
+
+    override fun sinh(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { sinh(it) }
+    override fun cosh(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { cosh(it) }
+    override fun tanh(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { tanh(it) }
+    override fun asinh(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { asinh(it) }
+    override fun acosh(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { acosh(it) }
+    override fun atanh(arg: NDStructure<Double>): NDBuffer<Double> = map(arg) { atanh(it) }
+}
+
+
+/**
+ * Fast element production using function inlining
+ */
+public inline fun BufferedNDField<Double, RealField>.produceInline(crossinline initializer: RealField.(Int) -> Double): NDBuffer<Double> {
+    contract { callsInPlace(initializer, InvocationKind.EXACTLY_ONCE) }
+    val array = DoubleArray(strides.linearSize) { offset -> RealField.initializer(offset) }
+    return NDBuffer(strides, RealBuffer(array))
+}
+
+public fun NDAlgebra.Companion.real(vararg shape: Int): RealNDField = RealNDField(shape)
+
+/**
+ * Produce a context for n-dimensional operations inside this real field
+ */
+public inline fun <R> RealField.nd(vararg shape: Int, action: RealNDField.() -> R): R {
+    contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
+    return RealNDField(shape).run(action)
+}
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/ShortNDRing.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/ShortNDRing.kt
new file mode 100644
index 000000000..8663c4249
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/ShortNDRing.kt
@@ -0,0 +1,36 @@
+package kscience.kmath.nd
+
+import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.operations.RingWithNumbers
+import kscience.kmath.operations.ShortRing
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.ShortBuffer
+import kotlin.contracts.InvocationKind
+import kotlin.contracts.contract
+
+@OptIn(UnstableKMathAPI::class)
+public open class ShortNDRing(
+    shape: IntArray,
+) : BufferedNDRing<Short, ShortRing>(shape, ShortRing, Buffer.Companion::auto),
+    RingWithNumbers<NDStructure<Short>> {
+
+    override val zero: NDBuffer<Short> by lazy { produce { zero } }
+    override val one: NDBuffer<Short> by lazy { produce { one } }
+
+    override fun number(value: Number): NDBuffer<Short> {
+        val d = value.toShort() // minimize conversions
+        return produce { d }
+    }
+}
+
+/**
+ * Fast element production using function inlining.
+ */
+public inline fun BufferedNDRing<Short, ShortRing>.produceInline(crossinline initializer: ShortRing.(Int) -> Short): NDBuffer<Short> {
+    return NDBuffer(strides, ShortBuffer(ShortArray(strides.linearSize) { offset -> ShortRing.initializer(offset) }))
+}
+
+public inline fun <R> ShortRing.nd(vararg shape: Int, action: ShortNDRing.() -> R): R {
+    contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
+    return ShortNDRing(shape).run(action)
+}
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure1D.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/Structure1D.kt
similarity index 89%
rename from kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure1D.kt
rename to kmath-core/src/commonMain/kotlin/kscience/kmath/nd/Structure1D.kt
index 95422ac60..d69c23fb4 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure1D.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/Structure1D.kt
@@ -1,4 +1,7 @@
-package kscience.kmath.structures
+package kscience.kmath.nd
+
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.asSequence
 
 /**
  * A structure that is guaranteed to be one-dimensional
@@ -34,7 +37,7 @@ private inline class Buffer1DWrapper<T>(val buffer: Buffer<T>) : Structure1D<T>
     override val size: Int get() = buffer.size
 
     override fun elements(): Sequence<Pair<IntArray, T>> =
-        asSequence().mapIndexed { index, value -> intArrayOf(index) to value }
+        buffer.asSequence().mapIndexed { index, value -> intArrayOf(index) to value }
 
     override operator fun get(index: Int): T = buffer[index]
 }
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure2D.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/Structure2D.kt
similarity index 92%
rename from kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure2D.kt
rename to kmath-core/src/commonMain/kotlin/kscience/kmath/nd/Structure2D.kt
index d20e9e53b..4c730877c 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Structure2D.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/Structure2D.kt
@@ -1,4 +1,7 @@
-package kscience.kmath.structures
+package kscience.kmath.nd
+
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.VirtualBuffer
 
 /**
  * A structure that is guaranteed to be two-dimensional.
@@ -73,10 +76,3 @@ public fun <T> NDStructure<T>.as2D(): Structure2D<T> = if (shape.size == 2)
     Structure2DWrapper(this)
 else
     error("Can't create 2d-structure from ${shape.size}d-structure")
-
-/**
- * Alias for [Structure2D] with more familiar name.
- *
- * @param T the type of items.
- */
-public typealias Matrix<T> = Structure2D<T>
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/BigInt.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/BigInt.kt
index 0be72e80c..9599401bc 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/BigInt.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/operations/BigInt.kt
@@ -1,6 +1,7 @@
 package kscience.kmath.operations
 
 import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.nd.NDAlgebra
 import kscience.kmath.operations.BigInt.Companion.BASE
 import kscience.kmath.operations.BigInt.Companion.BASE_SIZE
 import kscience.kmath.structures.*
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BoxingNDField.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BoxingNDField.kt
deleted file mode 100644
index dc65b12c4..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BoxingNDField.kt
+++ /dev/null
@@ -1,81 +0,0 @@
-package kscience.kmath.structures
-
-import kscience.kmath.operations.Field
-import kscience.kmath.operations.FieldElement
-
-public class BoxingNDField<T, F : Field<T>>(
-    public override val shape: IntArray,
-    public override val elementContext: F,
-    public val bufferFactory: BufferFactory<T>
-) : BufferedNDField<T, F> {
-    public override val zero: BufferedNDFieldElement<T, F> by lazy { produce { zero } }
-    public override val one: BufferedNDFieldElement<T, F> by lazy { produce { one } }
-    public override val strides: Strides = DefaultStrides(shape)
-
-    public fun buildBuffer(size: Int, initializer: (Int) -> T): Buffer<T> =
-        bufferFactory(size, initializer)
-
-    public override fun check(vararg elements: NDBuffer<T>): Array<out NDBuffer<T>> {
-        require(elements.all { it.strides == strides }) { "Element strides are not the same as context strides" }
-        return elements
-    }
-
-    public override fun produce(initializer: F.(IntArray) -> T): BufferedNDFieldElement<T, F> =
-        BufferedNDFieldElement(
-            this,
-            buildBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) })
-
-    public override fun map(arg: NDBuffer<T>, transform: F.(T) -> T): BufferedNDFieldElement<T, F> {
-        check(arg)
-
-        return BufferedNDFieldElement(
-            this,
-            buildBuffer(arg.strides.linearSize) { offset -> elementContext.transform(arg.buffer[offset]) })
-
-//        val buffer = arg.buffer.transform { _, value -> elementContext.transform(value) }
-//        return BufferedNDFieldElement(this, buffer)
-
-    }
-
-    public override fun mapIndexed(
-        arg: NDBuffer<T>,
-        transform: F.(index: IntArray, T) -> T
-    ): BufferedNDFieldElement<T, F> {
-        check(arg)
-        return BufferedNDFieldElement(
-            this,
-            buildBuffer(arg.strides.linearSize) { offset ->
-                elementContext.transform(
-                    arg.strides.index(offset),
-                    arg.buffer[offset]
-                )
-            })
-
-//        val buffer =
-//            arg.buffer.transform { offset, value -> elementContext.transform(arg.strides.index(offset), value) }
-//        return BufferedNDFieldElement(this, buffer)
-    }
-
-    public override fun combine(
-        a: NDBuffer<T>,
-        b: NDBuffer<T>,
-        transform: F.(T, T) -> T
-    ): BufferedNDFieldElement<T, F> {
-        check(a, b)
-        return BufferedNDFieldElement(
-            this,
-            buildBuffer(strides.linearSize) { offset -> elementContext.transform(a.buffer[offset], b.buffer[offset]) })
-    }
-
-    public override fun NDBuffer<T>.toElement(): FieldElement<NDBuffer<T>, *, out BufferedNDField<T, F>> =
-        BufferedNDFieldElement(this@BoxingNDField, buffer)
-}
-
-public inline fun <T : Any, F : Field<T>, R> F.nd(
-    noinline bufferFactory: BufferFactory<T>,
-    vararg shape: Int,
-    action: NDField<T, F, *>.() -> R
-): R {
-    val ndfield = NDField.boxing(this, *shape, bufferFactory = bufferFactory)
-    return ndfield.action()
-}
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BoxingNDRing.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BoxingNDRing.kt
deleted file mode 100644
index b6794984c..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BoxingNDRing.kt
+++ /dev/null
@@ -1,71 +0,0 @@
-package kscience.kmath.structures
-
-import kscience.kmath.operations.Ring
-import kscience.kmath.operations.RingElement
-
-public class BoxingNDRing<T, R : Ring<T>>(
-    override val shape: IntArray,
-    override val elementContext: R,
-    public val bufferFactory: BufferFactory<T>
-) : BufferedNDRing<T, R> {
-    override val strides: Strides = DefaultStrides(shape)
-    override val zero: BufferedNDRingElement<T, R> by lazy { produce { zero } }
-    override val one: BufferedNDRingElement<T, R> by lazy { produce { one } }
-
-    public fun buildBuffer(size: Int, initializer: (Int) -> T): Buffer<T> = bufferFactory(size, initializer)
-
-    override fun check(vararg elements: NDBuffer<T>): Array<out NDBuffer<T>> {
-        if (!elements.all { it.strides == this.strides }) error("Element strides are not the same as context strides")
-        return elements
-    }
-
-    override fun produce(initializer: R.(IntArray) -> T): BufferedNDRingElement<T, R> =
-        BufferedNDRingElement(
-            this,
-            buildBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) })
-
-    override fun map(arg: NDBuffer<T>, transform: R.(T) -> T): BufferedNDRingElement<T, R> {
-        check(arg)
-        return BufferedNDRingElement(
-            this,
-            buildBuffer(arg.strides.linearSize) { offset -> elementContext.transform(arg.buffer[offset]) })
-
-//        val buffer = arg.buffer.transform { _, value -> elementContext.transform(value) }
-//        return BufferedNDFieldElement(this, buffer)
-
-    }
-
-    override fun mapIndexed(
-        arg: NDBuffer<T>,
-        transform: R.(index: IntArray, T) -> T
-    ): BufferedNDRingElement<T, R> {
-        check(arg)
-        return BufferedNDRingElement(
-            this,
-            buildBuffer(arg.strides.linearSize) { offset ->
-                elementContext.transform(
-                    arg.strides.index(offset),
-                    arg.buffer[offset]
-                )
-            })
-
-//        val buffer =
-//            arg.buffer.transform { offset, value -> elementContext.transform(arg.strides.index(offset), value) }
-//        return BufferedNDFieldElement(this, buffer)
-    }
-
-    override fun combine(
-        a: NDBuffer<T>,
-        b: NDBuffer<T>,
-        transform: R.(T, T) -> T
-    ): BufferedNDRingElement<T, R> {
-        check(a, b)
-
-        return BufferedNDRingElement(
-            this,
-            buildBuffer(strides.linearSize) { offset -> elementContext.transform(a.buffer[offset], b.buffer[offset]) })
-    }
-
-    override fun NDBuffer<T>.toElement(): RingElement<NDBuffer<T>, *, out BufferedNDRing<T, R>> =
-        BufferedNDRingElement(this@BoxingNDRing, buffer)
-}
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt
index 5d7ba611f..3c53d6c51 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferAccessor2D.kt
@@ -1,5 +1,10 @@
 package kscience.kmath.structures
 
+import kscience.kmath.nd.DefaultStrides
+import kscience.kmath.nd.NDStructure
+import kscience.kmath.nd.Structure2D
+import kscience.kmath.nd.as2D
+
 /**
  * A context that allows to operate on a [MutableBuffer] as on 2d array
  */
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferedNDAlgebra.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferedNDAlgebra.kt
deleted file mode 100644
index 3dcd0322c..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferedNDAlgebra.kt
+++ /dev/null
@@ -1,43 +0,0 @@
-package kscience.kmath.structures
-
-import kscience.kmath.operations.*
-
-public interface BufferedNDAlgebra<T, C> : NDAlgebra<T, C, NDBuffer<T>> {
-    public val strides: Strides
-
-    public override fun check(vararg elements: NDBuffer<T>): Array<out NDBuffer<T>> {
-        require(elements.all { it.strides == strides }) { "Strides mismatch" }
-        return elements
-    }
-
-    /**
-     * Convert any [NDStructure] to buffered structure using strides from this context.
-     * If the structure is already [NDBuffer], conversion is free. If not, it could be expensive because iteration over
-     * indices.
-     *
-     * If the argument is [NDBuffer] with different strides structure, the new element will be produced.
-     */
-    public fun NDStructure<T>.toBuffer(): NDBuffer<T> =
-        if (this is NDBuffer<T> && this.strides == this@BufferedNDAlgebra.strides)
-            this
-        else
-            produce { index -> this@toBuffer[index] }
-
-    /**
-     * Convert a buffer to element of this algebra
-     */
-    public fun NDBuffer<T>.toElement(): MathElement<out BufferedNDAlgebra<T, C>>
-}
-
-
-public interface BufferedNDSpace<T, S : Space<T>> : NDSpace<T, S, NDBuffer<T>>, BufferedNDAlgebra<T, S> {
-    public override fun NDBuffer<T>.toElement(): SpaceElement<NDBuffer<T>, *, out BufferedNDSpace<T, S>>
-}
-
-public interface BufferedNDRing<T, R : Ring<T>> : NDRing<T, R, NDBuffer<T>>, BufferedNDSpace<T, R> {
-    override fun NDBuffer<T>.toElement(): RingElement<NDBuffer<T>, *, out BufferedNDRing<T, R>>
-}
-
-public interface BufferedNDField<T, F : Field<T>> : NDField<T, F, NDBuffer<T>>, BufferedNDRing<T, F> {
-    override fun NDBuffer<T>.toElement(): FieldElement<NDBuffer<T>, *, out BufferedNDField<T, F>>
-}
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferedNDElement.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferedNDElement.kt
deleted file mode 100644
index d53702566..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/BufferedNDElement.kt
+++ /dev/null
@@ -1,86 +0,0 @@
-package kscience.kmath.structures
-
-import kscience.kmath.operations.*
-
-/**
- * Base class for an element with context, containing strides
- */
-public abstract class BufferedNDElement<T, C> : NDBuffer<T>(), NDElement<T, C, NDBuffer<T>> {
-    abstract override val context: BufferedNDAlgebra<T, C>
-
-    override val strides: Strides get() = context.strides
-
-    override val shape: IntArray get() = context.shape
-}
-
-public class BufferedNDSpaceElement<T, S : Space<T>>(
-    override val context: BufferedNDSpace<T, S>,
-    override val buffer: Buffer<T>
-) : BufferedNDElement<T, S>(), SpaceElement<NDBuffer<T>, BufferedNDSpaceElement<T, S>, BufferedNDSpace<T, S>> {
-
-    override fun unwrap(): NDBuffer<T> = this
-
-    override fun NDBuffer<T>.wrap(): BufferedNDSpaceElement<T, S> {
-        context.check(this)
-        return BufferedNDSpaceElement(context, buffer)
-    }
-}
-
-public class BufferedNDRingElement<T, R : Ring<T>>(
-    override val context: BufferedNDRing<T, R>,
-    override val buffer: Buffer<T>
-) : BufferedNDElement<T, R>(), RingElement<NDBuffer<T>, BufferedNDRingElement<T, R>, BufferedNDRing<T, R>> {
-    override fun unwrap(): NDBuffer<T> = this
-
-    override fun NDBuffer<T>.wrap(): BufferedNDRingElement<T, R> {
-        context.check(this)
-        return BufferedNDRingElement(context, buffer)
-    }
-}
-
-public class BufferedNDFieldElement<T, F : Field<T>>(
-    override val context: BufferedNDField<T, F>,
-    override val buffer: Buffer<T>
-) : BufferedNDElement<T, F>(), FieldElement<NDBuffer<T>, BufferedNDFieldElement<T, F>, BufferedNDField<T, F>> {
-    override fun unwrap(): NDBuffer<T> = this
-
-    override fun NDBuffer<T>.wrap(): BufferedNDFieldElement<T, F> {
-        context.check(this)
-        return BufferedNDFieldElement(context, buffer)
-    }
-}
-
-
-/**
- * Element by element application of any operation on elements to the whole array. Just like in numpy.
- */
-public operator fun <T : Any, F : Field<T>> Function1<T, T>.invoke(ndElement: BufferedNDElement<T, F>): MathElement<out BufferedNDAlgebra<T, F>> =
-    ndElement.context.run { map(ndElement) { invoke(it) }.toElement() }
-
-/* plus and minus */
-
-/**
- * Summation operation for [BufferedNDElement] and single element
- */
-public operator fun <T : Any, F : Space<T>> BufferedNDElement<T, F>.plus(arg: T): NDElement<T, F, NDBuffer<T>> =
-    context.map(this) { it + arg }.wrap()
-
-/**
- * Subtraction operation between [BufferedNDElement] and single element
- */
-public operator fun <T : Any, F : Space<T>> BufferedNDElement<T, F>.minus(arg: T): NDElement<T, F, NDBuffer<T>> =
-    context.map(this) { it - arg }.wrap()
-
-/* prod and div */
-
-/**
- * Product operation for [BufferedNDElement] and single element
- */
-public operator fun <T : Any, F : Ring<T>> BufferedNDElement<T, F>.times(arg: T): NDElement<T, F, NDBuffer<T>> =
-    context.map(this) { it * arg }.wrap()
-
-/**
- * Division operation between [BufferedNDElement] and single element
- */
-public operator fun <T : Any, F : Field<T>> BufferedNDElement<T, F>.div(arg: T): NDElement<T, F, NDBuffer<T>> =
-    context.map(this) { it / arg }.wrap()
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt
index bfec6f871..4e7c0d1b1 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt
@@ -43,7 +43,7 @@ public interface Buffer<T> {
      * Checks content equality with another buffer.
      */
     public fun contentEquals(other: Buffer<*>): Boolean =
-        asSequence().mapIndexed { index, value -> value == other[index] }.all { it }
+        kscience.kmath.nd.mapIndexed { index, value -> value == other[index] }.all { it }
 
     public companion object {
         /**
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ComplexNDField.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ComplexNDField.kt
deleted file mode 100644
index 6de69cabe..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ComplexNDField.kt
+++ /dev/null
@@ -1,158 +0,0 @@
-package kscience.kmath.structures
-
-import kscience.kmath.misc.UnstableKMathAPI
-import kscience.kmath.operations.*
-import kotlin.contracts.InvocationKind
-import kotlin.contracts.contract
-
-public typealias ComplexNDElement = BufferedNDFieldElement<Complex, ComplexField>
-
-/**
- * An optimized nd-field for complex numbers
- */
-@OptIn(UnstableKMathAPI::class)
-public class ComplexNDField(override val shape: IntArray) :
-    BufferedNDField<Complex, ComplexField>,
-    ExtendedNDField<Complex, ComplexField, NDBuffer<Complex>>,
-    RingWithNumbers<NDBuffer<Complex>>{
-
-    override val strides: Strides = DefaultStrides(shape)
-    override val elementContext: ComplexField get() = ComplexField
-    override val zero: ComplexNDElement by lazy { produce { zero } }
-    override val one: ComplexNDElement by lazy { produce { one } }
-
-    override fun number(value: Number): NDBuffer<Complex> {
-        val c = value.toComplex()
-        return produce { c }
-    }
-
-    public inline fun buildBuffer(size: Int, crossinline initializer: (Int) -> Complex): Buffer<Complex> =
-        Buffer.complex(size) { initializer(it) }
-
-    /**
-     * Inline transform an NDStructure to another structure
-     */
-    override fun map(
-        arg: NDBuffer<Complex>,
-        transform: ComplexField.(Complex) -> Complex,
-    ): ComplexNDElement {
-        check(arg)
-        val array = buildBuffer(arg.strides.linearSize) { offset -> ComplexField.transform(arg.buffer[offset]) }
-        return BufferedNDFieldElement(this, array)
-    }
-
-    override fun produce(initializer: ComplexField.(IntArray) -> Complex): ComplexNDElement {
-        val array = buildBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) }
-        return BufferedNDFieldElement(this, array)
-    }
-
-    override fun mapIndexed(
-        arg: NDBuffer<Complex>,
-        transform: ComplexField.(index: IntArray, Complex) -> Complex,
-    ): ComplexNDElement {
-        check(arg)
-
-        return BufferedNDFieldElement(
-            this,
-            buildBuffer(arg.strides.linearSize) { offset ->
-                elementContext.transform(
-                    arg.strides.index(offset),
-                    arg.buffer[offset]
-                )
-            })
-    }
-
-    override fun combine(
-        a: NDBuffer<Complex>,
-        b: NDBuffer<Complex>,
-        transform: ComplexField.(Complex, Complex) -> Complex,
-    ): ComplexNDElement {
-        check(a, b)
-
-        return BufferedNDFieldElement(
-            this,
-            buildBuffer(strides.linearSize) { offset -> elementContext.transform(a.buffer[offset], b.buffer[offset]) })
-    }
-
-    override fun NDBuffer<Complex>.toElement(): FieldElement<NDBuffer<Complex>, *, out BufferedNDField<Complex, ComplexField>> =
-        BufferedNDFieldElement(this@ComplexNDField, buffer)
-
-    override fun power(arg: NDBuffer<Complex>, pow: Number): ComplexNDElement =
-        map(arg) { power(it, pow) }
-
-    override fun exp(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { exp(it) }
-    override fun ln(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { ln(it) }
-
-    override fun sin(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { sin(it) }
-    override fun cos(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { cos(it) }
-    override fun tan(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { tan(it) }
-    override fun asin(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { asin(it) }
-    override fun acos(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { acos(it) }
-    override fun atan(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { atan(it) }
-
-    override fun sinh(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { sinh(it) }
-    override fun cosh(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { cosh(it) }
-    override fun tanh(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { tanh(it) }
-    override fun asinh(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { asinh(it) }
-    override fun acosh(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { acosh(it) }
-    override fun atanh(arg: NDBuffer<Complex>): ComplexNDElement = map(arg) { atanh(it) }
-}
-
-
-/**
- * Fast element production using function inlining
- */
-public inline fun BufferedNDField<Complex, ComplexField>.produceInline(initializer: ComplexField.(Int) -> Complex): ComplexNDElement {
-    val buffer = Buffer.complex(strides.linearSize) { offset -> ComplexField.initializer(offset) }
-    return BufferedNDFieldElement(this, buffer)
-}
-
-/**
- * Map one [ComplexNDElement] using function with indices.
- */
-public inline fun ComplexNDElement.mapIndexed(transform: ComplexField.(index: IntArray, Complex) -> Complex): ComplexNDElement =
-    context.produceInline { offset -> transform(strides.index(offset), buffer[offset]) }
-
-/**
- * Map one [ComplexNDElement] using function without indices.
- */
-public inline fun ComplexNDElement.map(transform: ComplexField.(Complex) -> Complex): ComplexNDElement {
-    val buffer = Buffer.complex(strides.linearSize) { offset -> ComplexField.transform(buffer[offset]) }
-    return BufferedNDFieldElement(context, buffer)
-}
-
-/**
- * Element by element application of any operation on elements to the whole array. Just like in numpy
- */
-public operator fun Function1<Complex, Complex>.invoke(ndElement: ComplexNDElement): ComplexNDElement =
-    ndElement.map { this@invoke(it) }
-
-/* plus and minus */
-
-/**
- * Summation operation for [BufferedNDElement] and single element
- */
-public operator fun ComplexNDElement.plus(arg: Complex): ComplexNDElement = map { it + arg }
-
-/**
- * Subtraction operation between [BufferedNDElement] and single element
- */
-public operator fun ComplexNDElement.minus(arg: Complex): ComplexNDElement = map { it - arg }
-
-public operator fun ComplexNDElement.plus(arg: Double): ComplexNDElement = map { it + arg }
-public operator fun ComplexNDElement.minus(arg: Double): ComplexNDElement = map { it - arg }
-
-public fun NDField.Companion.complex(vararg shape: Int): ComplexNDField = ComplexNDField(shape)
-
-public fun NDElement.Companion.complex(
-    vararg shape: Int,
-    initializer: ComplexField.(IntArray) -> Complex,
-): ComplexNDElement = NDField.complex(*shape).produce(initializer)
-
-/**
- * Produce a context for n-dimensional operations inside this real field
- */
-public inline fun <R> ComplexField.nd(vararg shape: Int, action: ComplexNDField.() -> R): R {
-    contract { callsInPlace(action, InvocationKind.EXACTLY_ONCE) }
-    return NDField.complex(*shape).action()
-}
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ExtendedNDField.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ExtendedNDField.kt
deleted file mode 100644
index a9fa2763b..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ExtendedNDField.kt
+++ /dev/null
@@ -1,44 +0,0 @@
-package kscience.kmath.structures
-
-import kscience.kmath.operations.ExtendedField
-
-/**
- * [ExtendedField] over [NDStructure].
- *
- * @param T the type of the element contained in ND structure.
- * @param N the type of ND structure.
- * @param F the extended field of structure elements.
- */
-public interface ExtendedNDField<T : Any, F : ExtendedField<T>, N : NDStructure<T>> : NDField<T, F, N>, ExtendedField<N>
-
-///**
-// * NDField that supports [ExtendedField] operations on its elements
-// */
-//class ExtendedNDFieldWrapper<T : Any, F : ExtendedField<T>, N : NDStructure<T>>(private val ndField: NDField<T, F, N>) :
-//    ExtendedNDField<T, F, N>, NDField<T, F, N> by ndField {
-//
-//    override val shape: IntArray get() = ndField.shape
-//    override val elementContext: F get() = ndField.elementContext
-//
-//    override fun produce(initializer: F.(IntArray) -> T) = ndField.produce(initializer)
-//
-//    override fun power(arg: N, pow: Double): N {
-//        return produce { with(elementContext) { power(arg[it], pow) } }
-//    }
-//
-//    override fun exp(arg: N): N {
-//        return produce { with(elementContext) { exp(arg[it]) } }
-//    }
-//
-//    override fun ln(arg: N): N {
-//        return produce { with(elementContext) { ln(arg[it]) } }
-//    }
-//
-//    override fun sin(arg: N): N {
-//        return produce { with(elementContext) { sin(arg[it]) } }
-//    }
-//
-//    override fun cos(arg: N): N {
-//        return produce { with(elementContext) { cos(arg[it]) } }
-//    }
-//}
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/FlaggedBuffer.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/FlaggedBuffer.kt
index 4965e37cf..7828cea10 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/FlaggedBuffer.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/FlaggedBuffer.kt
@@ -60,7 +60,7 @@ public class FlaggedRealBuffer(public val values: DoubleArray, public val flags:
 
     override operator fun get(index: Int): Double? = if (isValid(index)) values[index] else null
 
-    override operator fun iterator(): Iterator<Double?> = values.indices.asSequence().map {
+    override operator fun iterator(): Iterator<Double?> = kscience.kmath.nd.map {
         if (isValid(it)) values[it] else null
     }.iterator()
 }
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/MemoryBuffer.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/MemoryBuffer.kt
index 66c9212cf..22331ae81 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/MemoryBuffer.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/MemoryBuffer.kt
@@ -15,7 +15,7 @@ public open class MemoryBuffer<T : Any>(protected val memory: Memory, protected
     private val reader: MemoryReader = memory.reader()
 
     override operator fun get(index: Int): T = reader.read(spec, spec.objectSize * index)
-    override operator fun iterator(): Iterator<T> = (0 until size).asSequence().map { get(it) }.iterator()
+    override operator fun iterator(): Iterator<T> = kscience.kmath.nd.map { get(it) }.iterator()
 
     public companion object {
         public fun <T : Any> create(spec: MemorySpec<T>, size: Int): MemoryBuffer<T> =
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDAlgebra.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDAlgebra.kt
deleted file mode 100644
index d7b019c65..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDAlgebra.kt
+++ /dev/null
@@ -1,259 +0,0 @@
-package kscience.kmath.structures
-
-import kscience.kmath.operations.Complex
-import kscience.kmath.operations.Field
-import kscience.kmath.operations.Ring
-import kscience.kmath.operations.Space
-import kotlin.native.concurrent.ThreadLocal
-
-/**
- * An exception is thrown when the expected ans actual shape of NDArray differs.
- *
- * @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()}.")
-
-/**
- * The base interface for all ND-algebra implementations.
- *
- * @param T the type of ND-structure element.
- * @param C the type of the element context.
- * @param N the type of the structure.
- */
-public interface NDAlgebra<T, C, N : NDStructure<T>> {
-    /**
-     * The shape of ND-structures this algebra operates on.
-     */
-    public val shape: IntArray
-
-    /**
-     * The algebra over elements of ND structure.
-     */
-    public val elementContext: C
-
-    /**
-     * Produces a new [N] structure using given initializer function.
-     */
-    public fun produce(initializer: C.(IntArray) -> T): N
-
-    /**
-     * Maps elements from one structure to another one by applying [transform] to them.
-     */
-    public fun map(arg: N, transform: C.(T) -> T): N
-
-    /**
-     * Maps elements from one structure to another one by applying [transform] to them alongside with their indices.
-     */
-    public fun mapIndexed(arg: N, transform: C.(index: IntArray, T) -> T): N
-
-    /**
-     * Combines two structures into one.
-     */
-    public fun combine(a: N, b: N, transform: C.(T, T) -> T): N
-
-    /**
-     * Checks if given element is consistent with this context.
-     *
-     * @param element the structure to check.
-     * @return the valid structure.
-     */
-    public fun check(element: N): N {
-        if (!element.shape.contentEquals(shape)) throw ShapeMismatchException(shape, element.shape)
-        return element
-    }
-
-    /**
-     * Checks if given elements are consistent with this context.
-     *
-     * @param elements the structures to check.
-     * @return the array of valid structures.
-     */
-    public fun check(vararg elements: N): Array<out N> = elements
-        .map(NDStructure<T>::shape)
-        .singleOrNull { !shape.contentEquals(it) }
-        ?.let<IntArray, Array<out N>> { throw ShapeMismatchException(shape, it) }
-        ?: elements
-
-    /**
-     * Element-wise invocation of function working on [T] on a [NDStructure].
-     */
-    public operator fun Function1<T, T>.invoke(structure: N): N = map(structure) { value -> this@invoke(value) }
-
-    public companion object
-}
-
-/**
- * Space of [NDStructure].
- *
- * @param T the type of the element contained in ND structure.
- * @param N the type of ND structure.
- * @param S the type of space of structure elements.
- */
-public interface NDSpace<T, S : Space<T>, N : NDStructure<T>> : Space<N>, NDAlgebra<T, S, N> {
-    /**
-     * Element-wise addition.
-     *
-     * @param a the addend.
-     * @param b the augend.
-     * @return the sum.
-     */
-    public override fun add(a: N, b: N): N = combine(a, b) { aValue, bValue -> add(aValue, bValue) }
-
-    /**
-     * Element-wise multiplication by scalar.
-     *
-     * @param a the multiplicand.
-     * @param k the multiplier.
-     * @return the product.
-     */
-    public override fun multiply(a: N, k: Number): N = map(a) { multiply(it, k) }
-
-    // TODO move to extensions after KEEP-176
-
-    /**
-     * Adds an ND structure to an element of it.
-     *
-     * @receiver the addend.
-     * @param arg the augend.
-     * @return the sum.
-     */
-    public operator fun N.plus(arg: T): N = map(this) { value -> add(arg, value) }
-
-    /**
-     * Subtracts an element from ND structure of it.
-     *
-     * @receiver the dividend.
-     * @param arg the divisor.
-     * @return the quotient.
-     */
-    public operator fun N.minus(arg: T): N = map(this) { value -> add(arg, -value) }
-
-    /**
-     * Adds an element to ND structure of it.
-     *
-     * @receiver the addend.
-     * @param arg the augend.
-     * @return the sum.
-     */
-    public operator fun T.plus(arg: N): N = map(arg) { value -> add(this@plus, value) }
-
-    /**
-     * Subtracts an ND structure from an element of it.
-     *
-     * @receiver the dividend.
-     * @param arg the divisor.
-     * @return the quotient.
-     */
-    public operator fun T.minus(arg: N): N = map(arg) { value -> add(-this@minus, value) }
-
-    public companion object
-}
-
-/**
- * Ring of [NDStructure].
- *
- * @param T the type of the element contained in ND structure.
- * @param N the type of ND structure.
- * @param R the type of ring of structure elements.
- */
-public interface NDRing<T, R : Ring<T>, N : NDStructure<T>> : Ring<N>, NDSpace<T, R, N> {
-    /**
-     * Element-wise multiplication.
-     *
-     * @param a the multiplicand.
-     * @param b the multiplier.
-     * @return the product.
-     */
-    public override fun multiply(a: N, b: N): N = combine(a, b) { aValue, bValue -> multiply(aValue, bValue) }
-
-    //TODO move to extensions after KEEP-176
-
-    /**
-     * Multiplies an ND structure by an element of it.
-     *
-     * @receiver the multiplicand.
-     * @param arg the multiplier.
-     * @return the product.
-     */
-    public operator fun N.times(arg: T): N = map(this) { value -> multiply(arg, value) }
-
-    /**
-     * Multiplies an element by a ND structure of it.
-     *
-     * @receiver the multiplicand.
-     * @param arg the multiplier.
-     * @return the product.
-     */
-    public operator fun T.times(arg: N): N = map(arg) { value -> multiply(this@times, value) }
-
-    public companion object
-}
-
-/**
- * Field of [NDStructure].
- *
- * @param T the type of the element contained in ND structure.
- * @param N the type of ND structure.
- * @param F the type field of structure elements.
- */
-public interface NDField<T, F : Field<T>, N : NDStructure<T>> : Field<N>, NDRing<T, F, N> {
-    /**
-     * Element-wise division.
-     *
-     * @param a the dividend.
-     * @param b the divisor.
-     * @return the quotient.
-     */
-    public override fun divide(a: N, b: N): N = combine(a, b) { aValue, bValue -> divide(aValue, bValue) }
-
-    //TODO move to extensions after KEEP-176
-    /**
-     * Divides an ND structure by an element of it.
-     *
-     * @receiver the dividend.
-     * @param arg the divisor.
-     * @return the quotient.
-     */
-    public operator fun N.div(arg: T): N = map(this) { value -> divide(arg, value) }
-
-    /**
-     * Divides an element by an ND structure of it.
-     *
-     * @receiver the dividend.
-     * @param arg the divisor.
-     * @return the quotient.
-     */
-    public operator fun T.div(arg: N): N = map(arg) { divide(it, this@div) }
-
-    @ThreadLocal
-    public companion object {
-        private val realNDFieldCache: MutableMap<IntArray, RealNDField> = hashMapOf()
-
-        /**
-         * Create a nd-field for [Double] values or pull it from cache if it was created previously.
-         */
-        public fun real(vararg shape: Int): RealNDField = realNDFieldCache.getOrPut(shape) { RealNDField(shape) }
-
-        /**
-         * Create an ND field with boxing generic buffer.
-         */
-        public fun <T : Any, F : Field<T>> boxing(
-            field: F,
-            vararg shape: Int,
-            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
-        ): BoxingNDField<T, F> = BoxingNDField(shape, field, bufferFactory)
-
-        /**
-         * Create a most suitable implementation for nd-field using reified class.
-         */
-        @Suppress("UNCHECKED_CAST")
-        public inline fun <reified T : Any, F : Field<T>> auto(field: F, vararg shape: Int): BufferedNDField<T, F> =
-            when {
-                T::class == Double::class -> real(*shape) as BufferedNDField<T, F>
-                T::class == Complex::class -> complex(*shape) as BufferedNDField<T, F>
-                else -> BoxingNDField(shape, field, Buffer.Companion::auto)
-            }
-    }
-}
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDElement.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDElement.kt
deleted file mode 100644
index f2f565064..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDElement.kt
+++ /dev/null
@@ -1,134 +0,0 @@
-package kscience.kmath.structures
-
-import kscience.kmath.operations.Field
-import kscience.kmath.operations.RealField
-import kscience.kmath.operations.Ring
-import kscience.kmath.operations.Space
-
-/**
- * The root for all [NDStructure] based algebra elements. Does not implement algebra element root because of problems with recursive self-types
- * @param T  the type of the element of the structure
- * @param C the type of the context for the element
- * @param N the type of the underlying [NDStructure]
- */
-public interface NDElement<T, C, N : NDStructure<T>> : NDStructure<T> {
-    public val context: NDAlgebra<T, C, N>
-
-    public fun unwrap(): N
-
-    public fun N.wrap(): NDElement<T, C, N>
-
-    public companion object {
-        /**
-         * Create a optimized NDArray of doubles
-         */
-        public fun real(shape: IntArray, initializer: RealField.(IntArray) -> Double = { 0.0 }): RealNDElement =
-            NDField.real(*shape).produce(initializer)
-
-        public inline fun real1D(dim: Int, crossinline initializer: (Int) -> Double = { _ -> 0.0 }): RealNDElement =
-            real(intArrayOf(dim)) { initializer(it[0]) }
-
-        public inline fun real2D(
-            dim1: Int,
-            dim2: Int,
-            crossinline initializer: (Int, Int) -> Double = { _, _ -> 0.0 }
-        ): RealNDElement = real(intArrayOf(dim1, dim2)) { initializer(it[0], it[1]) }
-
-        public inline fun real3D(
-            dim1: Int,
-            dim2: Int,
-            dim3: Int,
-            crossinline initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }
-        ): RealNDElement = real(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
-
-
-        /**
-         * Simple boxing NDArray
-         */
-        public fun <T : Any, F : Field<T>> boxing(
-            shape: IntArray,
-            field: F,
-            initializer: F.(IntArray) -> T
-        ): BufferedNDElement<T, F> {
-            val ndField = BoxingNDField(shape, field, Buffer.Companion::boxing)
-            return ndField.produce(initializer)
-        }
-
-        public inline fun <reified T : Any, F : Field<T>> auto(
-            shape: IntArray,
-            field: F,
-            noinline initializer: F.(IntArray) -> T
-        ): BufferedNDFieldElement<T, F> {
-            val ndField = NDField.auto(field, *shape)
-            return BufferedNDFieldElement(ndField, ndField.produce(initializer).buffer)
-        }
-    }
-}
-
-public fun <T, C, N : NDStructure<T>> NDElement<T, C, N>.mapIndexed(transform: C.(index: IntArray, T) -> T): NDElement<T, C, N> =
-    context.mapIndexed(unwrap(), transform).wrap()
-
-public fun <T, C, N : NDStructure<T>> NDElement<T, C, N>.map(transform: C.(T) -> T): NDElement<T, C, N> =
-    context.map(unwrap(), transform).wrap()
-
-/**
- * Element by element application of any operation on elements to the whole [NDElement]
- */
-public operator fun <T, C, N : NDStructure<T>> Function1<T, T>.invoke(ndElement: NDElement<T, C, N>): NDElement<T, C, N> =
-    ndElement.map { value -> this@invoke(value) }
-
-/* plus and minus */
-
-/**
- * Summation operation for [NDElement] and single element
- */
-public operator fun <T, S : Space<T>, N : NDStructure<T>> NDElement<T, S, N>.plus(arg: T): NDElement<T, S, N> =
-    map { value -> arg + value }
-
-/**
- * Subtraction operation between [NDElement] and single element
- */
-public operator fun <T, S : Space<T>, N : NDStructure<T>> NDElement<T, S, N>.minus(arg: T): NDElement<T, S, N> =
-    map { value -> arg - value }
-
-/* prod and div */
-
-/**
- * Product operation for [NDElement] and single element
- */
-public operator fun <T, R : Ring<T>, N : NDStructure<T>> NDElement<T, R, N>.times(arg: T): NDElement<T, R, N> =
-    map { value -> arg * value }
-
-/**
- * Division operation between [NDElement] and single element
- */
-public operator fun <T, F : Field<T>, N : NDStructure<T>> NDElement<T, F, N>.div(arg: T): NDElement<T, F, N> =
-    map { value -> arg / value }
-
-//    /**
-//     * Reverse sum operation
-//     */
-//    operator fun T.plus(arg: NDStructure<T>): NDElement<T, F> = produce { index ->
-//        field.run { this@plus + arg[index] }
-//    }
-//
-//    /**
-//     * Reverse minus operation
-//     */
-//    operator fun T.minus(arg: NDStructure<T>): NDElement<T, F> = produce { index ->
-//        field.run { this@minus - arg[index] }
-//    }
-//
-//    /**
-//     * Reverse product operation
-//     */
-//    operator fun T.times(arg: NDStructure<T>): NDElement<T, F> = produce { index ->
-//        field.run { this@times * arg[index] }
-//    }
-//
-//    /**
-//     * Reverse division operation
-//     */
-//    operator fun T.div(arg: NDStructure<T>): NDElement<T, F> = produce { index ->
-//        field.run { this@div / arg[index] }
-//    }
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/RealNDField.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/RealNDField.kt
deleted file mode 100644
index 60e6de440..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/RealNDField.kt
+++ /dev/null
@@ -1,140 +0,0 @@
-package kscience.kmath.structures
-
-import kscience.kmath.misc.UnstableKMathAPI
-import kscience.kmath.operations.FieldElement
-import kscience.kmath.operations.RealField
-import kscience.kmath.operations.RingWithNumbers
-
-public typealias RealNDElement = BufferedNDFieldElement<Double, RealField>
-
-@OptIn(UnstableKMathAPI::class)
-public class RealNDField(override val shape: IntArray) :
-    BufferedNDField<Double, RealField>,
-    ExtendedNDField<Double, RealField, NDBuffer<Double>>,
-    RingWithNumbers<NDBuffer<Double>> {
-
-    override val strides: Strides = DefaultStrides(shape)
-
-    override val elementContext: RealField get() = RealField
-    override val zero: RealNDElement by lazy { produce { zero } }
-    override val one: RealNDElement by lazy { produce { one } }
-
-    override fun number(value: Number): NDBuffer<Double> {
-        val d = value.toDouble()
-        return produce { d }
-    }
-
-    @Suppress("OVERRIDE_BY_INLINE")
-    override inline fun map(
-        arg: NDBuffer<Double>,
-        transform: RealField.(Double) -> Double,
-    ): RealNDElement {
-        check(arg)
-        val array = RealBuffer(arg.strides.linearSize) { offset -> RealField.transform(arg.buffer[offset]) }
-        return BufferedNDFieldElement(this, array)
-    }
-
-    @Suppress("OVERRIDE_BY_INLINE")
-    override inline fun produce(initializer: RealField.(IntArray) -> Double): RealNDElement {
-        val array = RealBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) }
-        return BufferedNDFieldElement(this, array)
-    }
-
-    @Suppress("OVERRIDE_BY_INLINE")
-    override inline fun mapIndexed(
-        arg: NDBuffer<Double>,
-        transform: RealField.(index: IntArray, Double) -> Double,
-    ): RealNDElement {
-        check(arg)
-        return BufferedNDFieldElement(
-            this,
-            RealBuffer(arg.strides.linearSize) { offset ->
-                elementContext.transform(
-                    arg.strides.index(offset),
-                    arg.buffer[offset]
-                )
-            })
-    }
-
-    @Suppress("OVERRIDE_BY_INLINE")
-    override inline fun combine(
-        a: NDBuffer<Double>,
-        b: NDBuffer<Double>,
-        transform: RealField.(Double, Double) -> Double,
-    ): RealNDElement {
-        check(a, b)
-        val buffer = RealBuffer(strides.linearSize) { offset ->
-            elementContext.transform(a.buffer[offset], b.buffer[offset])
-        }
-        return BufferedNDFieldElement(this, buffer)
-    }
-
-    override fun NDBuffer<Double>.toElement(): FieldElement<NDBuffer<Double>, *, out BufferedNDField<Double, RealField>> =
-        BufferedNDFieldElement(this@RealNDField, buffer)
-
-    override fun power(arg: NDBuffer<Double>, pow: Number): RealNDElement = map(arg) { power(it, pow) }
-
-    override fun exp(arg: NDBuffer<Double>): RealNDElement = map(arg) { exp(it) }
-
-    override fun ln(arg: NDBuffer<Double>): RealNDElement = map(arg) { ln(it) }
-
-    override fun sin(arg: NDBuffer<Double>): RealNDElement = map(arg) { sin(it) }
-    override fun cos(arg: NDBuffer<Double>): RealNDElement = map(arg) { cos(it) }
-    override fun tan(arg: NDBuffer<Double>): RealNDElement = map(arg) { tan(it) }
-    override fun asin(arg: NDBuffer<Double>): RealNDElement = map(arg) { asin(it) }
-    override fun acos(arg: NDBuffer<Double>): RealNDElement = map(arg) { acos(it) }
-    override fun atan(arg: NDBuffer<Double>): RealNDElement = map(arg) { atan(it) }
-
-    override fun sinh(arg: NDBuffer<Double>): RealNDElement = map(arg) { sinh(it) }
-    override fun cosh(arg: NDBuffer<Double>): RealNDElement = map(arg) { cosh(it) }
-    override fun tanh(arg: NDBuffer<Double>): RealNDElement = map(arg) { tanh(it) }
-    override fun asinh(arg: NDBuffer<Double>): RealNDElement = map(arg) { asinh(it) }
-    override fun acosh(arg: NDBuffer<Double>): RealNDElement = map(arg) { acosh(it) }
-    override fun atanh(arg: NDBuffer<Double>): RealNDElement = map(arg) { atanh(it) }
-}
-
-
-/**
- * Fast element production using function inlining
- */
-public inline fun BufferedNDField<Double, RealField>.produceInline(crossinline initializer: RealField.(Int) -> Double): RealNDElement {
-    val array = DoubleArray(strides.linearSize) { offset -> RealField.initializer(offset) }
-    return BufferedNDFieldElement(this, RealBuffer(array))
-}
-
-/**
- * Map one [RealNDElement] using function with indices.
- */
-public inline fun RealNDElement.mapIndexed(crossinline transform: RealField.(index: IntArray, Double) -> Double): RealNDElement =
-    context.produceInline { offset -> transform(strides.index(offset), buffer[offset]) }
-
-/**
- * Map one [RealNDElement] using function without indices.
- */
-public inline fun RealNDElement.map(crossinline transform: RealField.(Double) -> Double): RealNDElement {
-    val array = DoubleArray(strides.linearSize) { offset -> RealField.transform(buffer[offset]) }
-    return BufferedNDFieldElement(context, RealBuffer(array))
-}
-
-/**
- * Element by element application of any operation on elements to the whole array. Just like in numpy.
- */
-public operator fun Function1<Double, Double>.invoke(ndElement: RealNDElement): RealNDElement =
-    ndElement.map { this@invoke(it) }
-
-/* plus and minus */
-
-/**
- * Summation operation for [BufferedNDElement] and single element
- */
-public operator fun RealNDElement.plus(arg: Double): RealNDElement = map { it + arg }
-
-/**
- * Subtraction operation between [BufferedNDElement] and single element
- */
-public operator fun RealNDElement.minus(arg: Double): RealNDElement = map { it - arg }
-
-/**
- * Produce a context for n-dimensional operations inside this real field
- */
-public inline fun <R> RealField.nd(vararg shape: Int, action: RealNDField.() -> R): R = NDField.real(*shape).run(action)
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ShortNDRing.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ShortNDRing.kt
deleted file mode 100644
index 3b506a26a..000000000
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/ShortNDRing.kt
+++ /dev/null
@@ -1,93 +0,0 @@
-package kscience.kmath.structures
-
-import kscience.kmath.operations.RingElement
-import kscience.kmath.operations.ShortRing
-
-public typealias ShortNDElement = BufferedNDRingElement<Short, ShortRing>
-
-public class ShortNDRing(override val shape: IntArray) :
-    BufferedNDRing<Short, ShortRing> {
-
-    override val strides: Strides = DefaultStrides(shape)
-    override val elementContext: ShortRing get() = ShortRing
-    override val zero: ShortNDElement by lazy { produce { zero } }
-    override val one: ShortNDElement by lazy { produce { one } }
-
-    public inline fun buildBuffer(size: Int, crossinline initializer: (Int) -> Short): Buffer<Short> =
-        ShortBuffer(ShortArray(size) { initializer(it) })
-
-    /**
-     * Inline transform an NDStructure to
-     */
-    override fun map(
-        arg: NDBuffer<Short>,
-        transform: ShortRing.(Short) -> Short
-    ): ShortNDElement {
-        check(arg)
-        val array = buildBuffer(arg.strides.linearSize) { offset -> ShortRing.transform(arg.buffer[offset]) }
-        return BufferedNDRingElement(this, array)
-    }
-
-    override fun produce(initializer: ShortRing.(IntArray) -> Short): ShortNDElement {
-        val array = buildBuffer(strides.linearSize) { offset -> elementContext.initializer(strides.index(offset)) }
-        return BufferedNDRingElement(this, array)
-    }
-
-    override fun mapIndexed(
-        arg: NDBuffer<Short>,
-        transform: ShortRing.(index: IntArray, Short) -> Short
-    ): ShortNDElement {
-        check(arg)
-
-        return BufferedNDRingElement(
-            this,
-            buildBuffer(arg.strides.linearSize) { offset ->
-                elementContext.transform(
-                    arg.strides.index(offset),
-                    arg.buffer[offset]
-                )
-            })
-    }
-
-    override fun combine(
-        a: NDBuffer<Short>,
-        b: NDBuffer<Short>,
-        transform: ShortRing.(Short, Short) -> Short
-    ): ShortNDElement {
-        check(a, b)
-        return BufferedNDRingElement(
-            this,
-            buildBuffer(strides.linearSize) { offset -> elementContext.transform(a.buffer[offset], b.buffer[offset]) })
-    }
-
-    override fun NDBuffer<Short>.toElement(): RingElement<NDBuffer<Short>, *, out BufferedNDRing<Short, ShortRing>> =
-        BufferedNDRingElement(this@ShortNDRing, buffer)
-}
-
-
-/**
- * Fast element production using function inlining.
- */
-public inline fun BufferedNDRing<Short, ShortRing>.produceInline(crossinline initializer: ShortRing.(Int) -> Short): ShortNDElement =
-    BufferedNDRingElement(this, ShortBuffer(ShortArray(strides.linearSize) { offset -> ShortRing.initializer(offset) }))
-
-/**
- * Element by element application of any operation on elements to the whole array.
- */
-public operator fun Function1<Short, Short>.invoke(ndElement: ShortNDElement): ShortNDElement =
-    ndElement.context.produceInline { i -> invoke(ndElement.buffer[i]) }
-
-
-/* plus and minus */
-
-/**
- * Summation operation for [ShortNDElement] and single element.
- */
-public operator fun ShortNDElement.plus(arg: Short): ShortNDElement =
-    context.produceInline { i -> (buffer[i] + arg).toShort() }
-
-/**
- * Subtraction operation between [ShortNDElement] and single element.
- */
-public operator fun ShortNDElement.minus(arg: Short): ShortNDElement =
-    context.produceInline { i -> (buffer[i] - arg).toShort() }
diff --git a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt
index d7755dcb5..a0ac92172 100644
--- a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt
+++ b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/MatrixTest.kt
@@ -1,9 +1,8 @@
 package kscience.kmath.linear
 
+import kscience.kmath.nd.NDStructure
+import kscience.kmath.nd.as2D
 import kscience.kmath.operations.invoke
-import kscience.kmath.structures.Matrix
-import kscience.kmath.structures.NDStructure
-import kscience.kmath.structures.as2D
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
diff --git a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt
index 28dfe46ec..522e711b6 100644
--- a/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt
+++ b/kmath-core/src/commonTest/kotlin/kscience/kmath/linear/RealLUSolverTest.kt
@@ -1,6 +1,5 @@
 package kscience.kmath.linear
 
-import kscience.kmath.structures.Matrix
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
diff --git a/kmath-core/src/commonTest/kotlin/kscience/kmath/structures/NDFieldTest.kt b/kmath-core/src/commonTest/kotlin/kscience/kmath/structures/NDFieldTest.kt
index 1129a8a36..0f5f99b49 100644
--- a/kmath-core/src/commonTest/kotlin/kscience/kmath/structures/NDFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/kscience/kmath/structures/NDFieldTest.kt
@@ -1,5 +1,6 @@
 package kscience.kmath.structures
 
+import kscience.kmath.nd.NDField
 import kscience.kmath.operations.internal.FieldVerifier
 import kotlin.test.Test
 import kotlin.test.assertEquals
diff --git a/kmath-core/src/commonTest/kotlin/kscience/kmath/structures/NumberNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/kscience/kmath/structures/NumberNDFieldTest.kt
index 22a0d3629..baf3656fc 100644
--- a/kmath-core/src/commonTest/kotlin/kscience/kmath/structures/NumberNDFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/kscience/kmath/structures/NumberNDFieldTest.kt
@@ -1,7 +1,8 @@
 package kscience.kmath.structures
 
+import kscience.kmath.nd.NDField
+import kscience.kmath.nd.NDStructure
 import kscience.kmath.operations.Norm
-import kscience.kmath.operations.invoke
 import kscience.kmath.structures.NDElement.Companion.real2D
 import kotlin.math.abs
 import kotlin.math.pow
diff --git a/kmath-coroutines/src/jvmMain/kotlin/kscience/kmath/structures/LazyNDStructure.kt b/kmath-coroutines/src/jvmMain/kotlin/kscience/kmath/structures/LazyNDStructure.kt
index 7aa746797..cac692eb2 100644
--- a/kmath-coroutines/src/jvmMain/kotlin/kscience/kmath/structures/LazyNDStructure.kt
+++ b/kmath-coroutines/src/jvmMain/kotlin/kscience/kmath/structures/LazyNDStructure.kt
@@ -2,6 +2,8 @@ package kscience.kmath.structures
 
 import kotlinx.coroutines.*
 import kscience.kmath.coroutines.Math
+import kscience.kmath.nd.DefaultStrides
+import kscience.kmath.nd.NDStructure
 
 public class LazyNDStructure<T>(
     public val scope: CoroutineScope,
@@ -19,7 +21,7 @@ public class LazyNDStructure<T>(
 
     public override fun elements(): Sequence<Pair<IntArray, T>> {
         val strides = DefaultStrides(shape)
-        val res = runBlocking { strides.indices().toList().map { index -> index to await(index) } }
+        val res = runBlocking { kscience.kmath.nd.map { index -> index to await(index) } }
         return res.asSequence()
     }
 
diff --git a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
index 0244eae7f..52443da0e 100644
--- a/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
+++ b/kmath-dimensions/src/commonMain/kotlin/kscience/kmath/dimensions/Wrappers.kt
@@ -1,9 +1,8 @@
 package kscience.kmath.dimensions
 
 import kscience.kmath.linear.*
+import kscience.kmath.nd.Structure2D
 import kscience.kmath.operations.invoke
-import kscience.kmath.structures.Matrix
-import kscience.kmath.structures.Structure2D
 
 /**
  * A matrix with compile-time controlled dimension
diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt
index 82a5399fd..13e5f009b 100644
--- a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt
+++ b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrix.kt
@@ -2,8 +2,7 @@ package kscience.kmath.ejml
 
 import kscience.kmath.linear.*
 import kscience.kmath.misc.UnstableKMathAPI
-import kscience.kmath.structures.Matrix
-import kscience.kmath.structures.NDStructure
+import kscience.kmath.nd.NDStructure
 import kscience.kmath.structures.RealBuffer
 import org.ejml.dense.row.factory.DecompositionFactory_DDRM
 import org.ejml.simple.SimpleMatrix
diff --git a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
index 8184d0110..ebd4dd586 100644
--- a/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
+++ b/kmath-ejml/src/main/kotlin/kscience/kmath/ejml/EjmlMatrixContext.kt
@@ -3,10 +3,9 @@ package kscience.kmath.ejml
 import kscience.kmath.linear.InverseMatrixFeature
 import kscience.kmath.linear.MatrixContext
 import kscience.kmath.linear.Point
-import kscience.kmath.linear.origin
 import kscience.kmath.misc.UnstableKMathAPI
-import kscience.kmath.structures.Matrix
-import kscience.kmath.structures.getFeature
+import kscience.kmath.nd.Matrix
+import kscience.kmath.nd.getFeature
 import org.ejml.simple.SimpleMatrix
 
 /**
diff --git a/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt b/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt
index 455b52d9d..6af70b6f5 100644
--- a/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt
+++ b/kmath-ejml/src/test/kotlin/kscience/kmath/ejml/EjmlMatrixTest.kt
@@ -5,7 +5,7 @@ import kscience.kmath.linear.LupDecompositionFeature
 import kscience.kmath.linear.MatrixFeature
 import kscience.kmath.linear.plus
 import kscience.kmath.misc.UnstableKMathAPI
-import kscience.kmath.structures.getFeature
+import kscience.kmath.nd.getFeature
 import org.ejml.dense.row.factory.DecompositionFactory_DDRM
 import org.ejml.simple.SimpleMatrix
 import kotlin.random.Random
diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
index 274030aff..defaba126 100644
--- a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
+++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/RealMatrix.kt
@@ -1,15 +1,11 @@
 package kscience.kmath.real
 
 import kscience.kmath.linear.MatrixContext
-import kscience.kmath.linear.VirtualMatrix
-import kscience.kmath.linear.inverseWithLUP
 import kscience.kmath.linear.real
 import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.nd.Matrix
 import kscience.kmath.structures.Buffer
-import kscience.kmath.structures.Matrix
 import kscience.kmath.structures.RealBuffer
-import kscience.kmath.structures.asIterable
-import kotlin.math.pow
 
 /*
  *  Functions for convenient "numpy-like" operations with Double matrices.
diff --git a/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/realND.kt b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/realND.kt
new file mode 100644
index 000000000..f8757b132
--- /dev/null
+++ b/kmath-for-real/src/commonMain/kotlin/kscience/kmath/real/realND.kt
@@ -0,0 +1,31 @@
+package kscience.kmath.real
+
+import kscience.kmath.nd.NDBuffer
+import kscience.kmath.operations.RealField
+import kscience.kmath.structures.RealBuffer
+
+/**
+ * Map one [NDBuffer] using function without indices.
+ */
+public inline fun NDBuffer<Double>.mapInline(crossinline transform: RealField.(Double) -> Double): NDBuffer<Double> {
+    val array = DoubleArray(strides.linearSize) { offset -> RealField.transform(buffer[offset]) }
+    return NDBuffer(strides, RealBuffer(array))
+}
+
+/**
+ * Element by element application of any operation on elements to the whole array. Just like in numpy.
+ */
+public operator fun Function1<Double, Double>.invoke(ndElement: NDBuffer<Double>): NDBuffer<Double> =
+    ndElement.mapInline { this@invoke(it) }
+
+/* plus and minus */
+
+/**
+ * Summation operation for [NDBuffer] and single element
+ */
+public operator fun NDBuffer<Double>.plus(arg: Double): NDBuffer<Double> = mapInline { it + arg }
+
+/**
+ * Subtraction operation between [NDBuffer] and single element
+ */
+public operator fun NDBuffer<Double>.minus(arg: Double): NDBuffer<Double> = mapInline { it - arg }
\ No newline at end of file
diff --git a/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealMatrixTest.kt b/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealMatrixTest.kt
index a89f99b3c..68aa8cf52 100644
--- a/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealMatrixTest.kt
+++ b/kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/RealMatrixTest.kt
@@ -1,9 +1,7 @@
 package kaceince.kmath.real
 
-import kscience.kmath.linear.build
+import kscience.kmath.nd.Matrix
 import kscience.kmath.real.*
-import kscience.kmath.structures.Matrix
-import kscience.kmath.structures.contentEquals
 import kotlin.test.Test
 import kotlin.test.assertEquals
 import kotlin.test.assertTrue
diff --git a/kmath-functions/src/commonMain/kotlin/kscience/kmath/interpolation/XYPointSet.kt b/kmath-functions/src/commonMain/kotlin/kscience/kmath/interpolation/XYPointSet.kt
index 2abb7742c..bdece28e7 100644
--- a/kmath-functions/src/commonMain/kotlin/kscience/kmath/interpolation/XYPointSet.kt
+++ b/kmath-functions/src/commonMain/kotlin/kscience/kmath/interpolation/XYPointSet.kt
@@ -1,7 +1,7 @@
 package kscience.kmath.interpolation
 
+import kscience.kmath.nd.Structure2D
 import kscience.kmath.structures.Buffer
-import kscience.kmath.structures.Structure2D
 
 public interface XYPointSet<X, Y> {
     public val size: Int
diff --git a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt
index f95264ee1..42dbaabac 100644
--- a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt
+++ b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt
@@ -1,6 +1,8 @@
 package kscience.kmath.histogram
 
 import kscience.kmath.linear.Point
+import kscience.kmath.nd.DefaultStrides
+import kscience.kmath.nd.NDStructure
 import kscience.kmath.operations.SpaceOperations
 import kscience.kmath.operations.invoke
 import kscience.kmath.structures.*
@@ -15,7 +17,7 @@ public data class BinDef<T : Comparable<T>>(
         require(vector.size == center.size) { "Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}" }
         val upper = space { center + sizes / 2.0 }
         val lower = space { center - sizes / 2.0 }
-        return vector.asSequence().mapIndexed { i, value -> value in lower[i]..upper[i] }.all { it }
+        return kscience.kmath.nd.mapIndexed { i, value -> value in lower[i]..upper[i] }.all { it }
     }
 }
 
@@ -68,7 +70,7 @@ public class RealHistogram(
     public fun getValue(point: Buffer<out Double>): Long = getValue(getIndex(point))
 
     private fun getDef(index: IntArray): BinDef<Double> {
-        val center = index.mapIndexed { axis, i ->
+        val center = kscience.kmath.nd.mapIndexed { axis, i ->
             when (i) {
                 0 -> Double.NEGATIVE_INFINITY
                 strides.shape[axis] - 1 -> Double.POSITIVE_INFINITY
@@ -98,7 +100,7 @@ public class RealHistogram(
     }
 
     public override operator fun iterator(): Iterator<MultivariateBin<Double>> =
-        weights.elements().map { (index, value) -> MultivariateBin(getDef(index), value.sum()) }
+        kscience.kmath.nd.map { (index, value) -> MultivariateBin(getDef(index), value.sum()) }
             .iterator()
 
     /**
diff --git a/kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayAlgebra.kt b/kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayAlgebra.kt
index db2a44861..b958d1ccf 100644
--- a/kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayAlgebra.kt
+++ b/kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayAlgebra.kt
@@ -1,11 +1,9 @@
 package kscience.kmath.nd4j
 
 import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.nd.*
 import kscience.kmath.operations.*
-import kscience.kmath.structures.NDAlgebra
-import kscience.kmath.structures.NDField
-import kscience.kmath.structures.NDRing
-import kscience.kmath.structures.NDSpace
+import kscience.kmath.structures.*
 import org.nd4j.linalg.api.ndarray.INDArray
 import org.nd4j.linalg.factory.Nd4j
 
@@ -15,11 +13,22 @@ import org.nd4j.linalg.factory.Nd4j
  * @param T the type of ND-structure element.
  * @param C the type of the element context.
  */
-public interface Nd4jArrayAlgebra<T, C> : NDAlgebra<T, C, Nd4jArrayStructure<T>> {
+public interface Nd4jArrayAlgebra<T, C> : NDAlgebra<T, C> {
     /**
      * Wraps [INDArray] to [N].
      */
     public fun INDArray.wrap(): Nd4jArrayStructure<T>
+    public val NDStructure<T>.ndArray: INDArray
+        get() = when  {
+            !shape.contentEquals(this@Nd4jArrayAlgebra.shape) -> throw ShapeMismatchException(
+                this@Nd4jArrayAlgebra.shape,
+                shape
+            )
+            this is Nd4jArrayStructure -> ndArray //TODO check strides
+            else -> {
+                TODO()
+            }
+        }
 
     public override fun produce(initializer: C.(IntArray) -> T): Nd4jArrayStructure<T> {
         val struct = Nd4j.create(*shape)!!.wrap()
@@ -27,29 +36,26 @@ public interface Nd4jArrayAlgebra<T, C> : NDAlgebra<T, C, Nd4jArrayStructure<T>>
         return struct
     }
 
-    public override fun map(arg: Nd4jArrayStructure<T>, transform: C.(T) -> T): Nd4jArrayStructure<T> {
-        check(arg)
+    public override fun map(arg: NDStructure<T>, transform: C.(T) -> T): Nd4jArrayStructure<T> {
         val newStruct = arg.ndArray.dup().wrap()
         newStruct.elements().forEach { (idx, value) -> newStruct[idx] = elementContext.transform(value) }
         return newStruct
     }
 
     public override fun mapIndexed(
-        arg: Nd4jArrayStructure<T>,
+        arg: NDStructure<T>,
         transform: C.(index: IntArray, T) -> T,
     ): Nd4jArrayStructure<T> {
-        check(arg)
         val new = Nd4j.create(*shape).wrap()
         new.indicesIterator().forEach { idx -> new[idx] = elementContext.transform(idx, arg[idx]) }
         return new
     }
 
     public override fun combine(
-        a: Nd4jArrayStructure<T>,
-        b: Nd4jArrayStructure<T>,
+        a: NDStructure<T>,
+        b: NDStructure<T>,
         transform: C.(T, T) -> T,
     ): Nd4jArrayStructure<T> {
-        check(a, b)
         val new = Nd4j.create(*shape).wrap()
         new.indicesIterator().forEach { idx -> new[idx] = elementContext.transform(a[idx], b[idx]) }
         return new
@@ -62,38 +68,32 @@ public interface Nd4jArrayAlgebra<T, C> : NDAlgebra<T, C, Nd4jArrayStructure<T>>
  * @param T the type of the element contained in ND structure.
  * @param S the type of space of structure elements.
  */
-public interface Nd4jArraySpace<T, S : Space<T>> : NDSpace<T, S, Nd4jArrayStructure<T>>, Nd4jArrayAlgebra<T, S> {
+public interface Nd4jArraySpace<T, S : Space<T>> : NDSpace<T, S>, Nd4jArrayAlgebra<T, S> {
 
     public override val zero: Nd4jArrayStructure<T>
         get() = Nd4j.zeros(*shape).wrap()
 
-    public override fun add(a: Nd4jArrayStructure<T>, b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
-        check(a, b)
+    public override fun add(a: NDStructure<T>, b: NDStructure<T>): Nd4jArrayStructure<T> {
         return a.ndArray.add(b.ndArray).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<T>.minus(b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
-        check(this, b)
+    public override operator fun NDStructure<T>.minus(b: NDStructure<T>): Nd4jArrayStructure<T> {
         return ndArray.sub(b.ndArray).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<T>.unaryMinus(): Nd4jArrayStructure<T> {
-        check(this)
+    public override operator fun NDStructure<T>.unaryMinus(): Nd4jArrayStructure<T> {
         return ndArray.neg().wrap()
     }
 
-    public override fun multiply(a: Nd4jArrayStructure<T>, k: Number): Nd4jArrayStructure<T> {
-        check(a)
+    public override fun multiply(a: NDStructure<T>, k: Number): Nd4jArrayStructure<T> {
         return a.ndArray.mul(k).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<T>.div(k: Number): Nd4jArrayStructure<T> {
-        check(this)
+    public override operator fun NDStructure<T>.div(k: Number): Nd4jArrayStructure<T> {
         return ndArray.div(k).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<T>.times(k: Number): Nd4jArrayStructure<T> {
-        check(this)
+    public override operator fun NDStructure<T>.times(k: Number): Nd4jArrayStructure<T> {
         return ndArray.mul(k).wrap()
     }
 }
@@ -105,13 +105,12 @@ public interface Nd4jArraySpace<T, S : Space<T>> : NDSpace<T, S, Nd4jArrayStruct
  * @param R the type of ring of structure elements.
  */
 @OptIn(UnstableKMathAPI::class)
-public interface Nd4jArrayRing<T, R : Ring<T>> : NDRing<T, R, Nd4jArrayStructure<T>>, Nd4jArraySpace<T, R> {
+public interface Nd4jArrayRing<T, R : Ring<T>> : NDRing<T, R>, Nd4jArraySpace<T, R> {
 
     public override val one: Nd4jArrayStructure<T>
         get() = Nd4j.ones(*shape).wrap()
 
-    public override fun multiply(a: Nd4jArrayStructure<T>, b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
-        check(a, b)
+    public override fun multiply(a: NDStructure<T>, b: NDStructure<T>): Nd4jArrayStructure<T> {
         return a.ndArray.mul(b.ndArray).wrap()
     }
 //
@@ -168,17 +167,12 @@ public interface Nd4jArrayRing<T, R : Ring<T>> : NDRing<T, R, Nd4jArrayStructure
  * @param N the type of ND structure.
  * @param F the type field of structure elements.
  */
-public interface Nd4jArrayField<T, F : Field<T>> : NDField<T, F, Nd4jArrayStructure<T>>, Nd4jArrayRing<T, F> {
+public interface Nd4jArrayField<T, F : Field<T>> : NDField<T, F>, Nd4jArrayRing<T, F> {
 
-    public override fun divide(a: Nd4jArrayStructure<T>, b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
-        check(a, b)
-        return a.ndArray.div(b.ndArray).wrap()
-    }
+    public override fun divide(a: NDStructure<T>, b: NDStructure<T>): Nd4jArrayStructure<T> =
+        a.ndArray.div(b.ndArray).wrap()
 
-    public override operator fun Number.div(b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
-        check(b)
-        return b.ndArray.rdiv(this).wrap()
-    }
+    public override operator fun Number.div(b: NDStructure<T>): Nd4jArrayStructure<T> = b.ndArray.rdiv(this).wrap()
 
 
     public companion object {
@@ -219,35 +213,29 @@ public class RealNd4jArrayField(public override val shape: IntArray) : Nd4jArray
     public override val elementContext: RealField
         get() = RealField
 
-    public override fun INDArray.wrap(): Nd4jArrayStructure<Double> = check(asRealStructure())
+    public override fun INDArray.wrap(): Nd4jArrayStructure<Double> = checkShape(asRealStructure())
 
-    public override operator fun Nd4jArrayStructure<Double>.div(arg: Double): Nd4jArrayStructure<Double> {
-        check(this)
+    public override operator fun NDStructure<Double>.div(arg: Double): Nd4jArrayStructure<Double> {
         return ndArray.div(arg).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<Double>.plus(arg: Double): Nd4jArrayStructure<Double> {
-        check(this)
+    public override operator fun NDStructure<Double>.plus(arg: Double): Nd4jArrayStructure<Double> {
         return ndArray.add(arg).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<Double>.minus(arg: Double): Nd4jArrayStructure<Double> {
-        check(this)
+    public override operator fun NDStructure<Double>.minus(arg: Double): Nd4jArrayStructure<Double> {
         return ndArray.sub(arg).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<Double>.times(arg: Double): Nd4jArrayStructure<Double> {
-        check(this)
+    public override operator fun NDStructure<Double>.times(arg: Double): Nd4jArrayStructure<Double> {
         return ndArray.mul(arg).wrap()
     }
 
-    public override operator fun Double.div(arg: Nd4jArrayStructure<Double>): Nd4jArrayStructure<Double> {
-        check(arg)
+    public override operator fun Double.div(arg: NDStructure<Double>): Nd4jArrayStructure<Double> {
         return arg.ndArray.rdiv(this).wrap()
     }
 
-    public override operator fun Double.minus(arg: Nd4jArrayStructure<Double>): Nd4jArrayStructure<Double> {
-        check(arg)
+    public override operator fun Double.minus(arg: NDStructure<Double>): Nd4jArrayStructure<Double> {
         return arg.ndArray.rsub(this).wrap()
     }
 }
@@ -259,35 +247,29 @@ public class FloatNd4jArrayField(public override val shape: IntArray) : Nd4jArra
     public override val elementContext: FloatField
         get() = FloatField
 
-    public override fun INDArray.wrap(): Nd4jArrayStructure<Float> = check(asFloatStructure())
+    public override fun INDArray.wrap(): Nd4jArrayStructure<Float> = checkShape(asFloatStructure())
 
-    public override operator fun Nd4jArrayStructure<Float>.div(arg: Float): Nd4jArrayStructure<Float> {
-        check(this)
+    public override operator fun NDStructure<Float>.div(arg: Float): Nd4jArrayStructure<Float> {
         return ndArray.div(arg).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<Float>.plus(arg: Float): Nd4jArrayStructure<Float> {
-        check(this)
+    public override operator fun NDStructure<Float>.plus(arg: Float): Nd4jArrayStructure<Float> {
         return ndArray.add(arg).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<Float>.minus(arg: Float): Nd4jArrayStructure<Float> {
-        check(this)
+    public override operator fun NDStructure<Float>.minus(arg: Float): Nd4jArrayStructure<Float> {
         return ndArray.sub(arg).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<Float>.times(arg: Float): Nd4jArrayStructure<Float> {
-        check(this)
+    public override operator fun NDStructure<Float>.times(arg: Float): Nd4jArrayStructure<Float> {
         return ndArray.mul(arg).wrap()
     }
 
-    public override operator fun Float.div(arg: Nd4jArrayStructure<Float>): Nd4jArrayStructure<Float> {
-        check(arg)
+    public override operator fun Float.div(arg: NDStructure<Float>): Nd4jArrayStructure<Float> {
         return arg.ndArray.rdiv(this).wrap()
     }
 
-    public override operator fun Float.minus(arg: Nd4jArrayStructure<Float>): Nd4jArrayStructure<Float> {
-        check(arg)
+    public override operator fun Float.minus(arg: NDStructure<Float>): Nd4jArrayStructure<Float> {
         return arg.ndArray.rsub(this).wrap()
     }
 }
@@ -301,23 +283,19 @@ public class IntNd4jArrayRing(public override val shape: IntArray) : Nd4jArrayRi
 
     public override fun INDArray.wrap(): Nd4jArrayStructure<Int> = check(asIntStructure())
 
-    public override operator fun Nd4jArrayStructure<Int>.plus(arg: Int): Nd4jArrayStructure<Int> {
-        check(this)
+    public override operator fun NDStructure<Int>.plus(arg: Int): Nd4jArrayStructure<Int> {
         return ndArray.add(arg).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<Int>.minus(arg: Int): Nd4jArrayStructure<Int> {
-        check(this)
+    public override operator fun NDStructure<Int>.minus(arg: Int): Nd4jArrayStructure<Int> {
         return ndArray.sub(arg).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<Int>.times(arg: Int): Nd4jArrayStructure<Int> {
-        check(this)
+    public override operator fun NDStructure<Int>.times(arg: Int): Nd4jArrayStructure<Int> {
         return ndArray.mul(arg).wrap()
     }
 
-    public override operator fun Int.minus(arg: Nd4jArrayStructure<Int>): Nd4jArrayStructure<Int> {
-        check(arg)
+    public override operator fun Int.minus(arg: NDStructure<Int>): Nd4jArrayStructure<Int> {
         return arg.ndArray.rsub(this).wrap()
     }
 }
@@ -331,23 +309,19 @@ public class LongNd4jArrayRing(public override val shape: IntArray) : Nd4jArrayR
 
     public override fun INDArray.wrap(): Nd4jArrayStructure<Long> = check(asLongStructure())
 
-    public override operator fun Nd4jArrayStructure<Long>.plus(arg: Long): Nd4jArrayStructure<Long> {
-        check(this)
+    public override operator fun NDStructure<Long>.plus(arg: Long): Nd4jArrayStructure<Long> {
         return ndArray.add(arg).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<Long>.minus(arg: Long): Nd4jArrayStructure<Long> {
-        check(this)
+    public override operator fun NDStructure<Long>.minus(arg: Long): Nd4jArrayStructure<Long> {
         return ndArray.sub(arg).wrap()
     }
 
-    public override operator fun Nd4jArrayStructure<Long>.times(arg: Long): Nd4jArrayStructure<Long> {
-        check(this)
+    public override operator fun NDStructure<Long>.times(arg: Long): Nd4jArrayStructure<Long> {
         return ndArray.mul(arg).wrap()
     }
 
-    public override operator fun Long.minus(arg: Nd4jArrayStructure<Long>): Nd4jArrayStructure<Long> {
-        check(arg)
+    public override operator fun Long.minus(arg: NDStructure<Long>): Nd4jArrayStructure<Long> {
         return arg.ndArray.rsub(this).wrap()
     }
 }
diff --git a/kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayStructure.kt b/kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayStructure.kt
index d47a293c3..2a6b56602 100644
--- a/kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayStructure.kt
+++ b/kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayStructure.kt
@@ -1,7 +1,7 @@
 package kscience.kmath.nd4j
 
-import kscience.kmath.structures.MutableNDStructure
-import kscience.kmath.structures.NDStructure
+import kscience.kmath.nd.MutableNDStructure
+import kscience.kmath.nd.NDStructure
 import org.nd4j.linalg.api.ndarray.INDArray
 
 /**
diff --git a/kmath-nd4j/src/test/kotlin/kscience/kmath/nd4j/Nd4jArrayStructureTest.kt b/kmath-nd4j/src/test/kotlin/kscience/kmath/nd4j/Nd4jArrayStructureTest.kt
index 7e46211c1..2f1606061 100644
--- a/kmath-nd4j/src/test/kotlin/kscience/kmath/nd4j/Nd4jArrayStructureTest.kt
+++ b/kmath-nd4j/src/test/kotlin/kscience/kmath/nd4j/Nd4jArrayStructureTest.kt
@@ -1,6 +1,6 @@
 package kscience.kmath.nd4j
 
-import kscience.kmath.structures.get
+import kscience.kmath.nd.get
 import org.nd4j.linalg.factory.Nd4j
 import kotlin.test.Test
 import kotlin.test.assertEquals
diff --git a/kmath-viktor/src/main/kotlin/kscience/kmath/viktor/ViktorNDStructure.kt b/kmath-viktor/src/main/kotlin/kscience/kmath/viktor/ViktorNDStructure.kt
index 2471362cb..5dbc297b4 100644
--- a/kmath-viktor/src/main/kotlin/kscience/kmath/viktor/ViktorNDStructure.kt
+++ b/kmath-viktor/src/main/kotlin/kscience/kmath/viktor/ViktorNDStructure.kt
@@ -1,10 +1,10 @@
 package kscience.kmath.viktor
 
+import kscience.kmath.nd.DefaultStrides
+import kscience.kmath.nd.MutableNDStructure
+import kscience.kmath.nd.NDField
+import kscience.kmath.nd.Strides
 import kscience.kmath.operations.RealField
-import kscience.kmath.structures.DefaultStrides
-import kscience.kmath.structures.MutableNDStructure
-import kscience.kmath.structures.NDField
-import kscience.kmath.structures.Strides
 import org.jetbrains.bio.viktor.F64Array
 
 @Suppress("OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")