From 20a7c6c4f11f334a2bd9eab37aab04e790c26f9a Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Mon, 5 Nov 2018 16:11:53 +0300
Subject: [PATCH 01/13] MultivariateHistogram fix

---
 .../kotlin/scientifik/kmath/histogram/FastHistogram.kt     | 6 +++---
 .../kmath/histogram/MultivariateHistogramTest.kt           | 7 ++++---
 2 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
index 5f69f9d12..3852de316 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
@@ -17,7 +17,7 @@ class MultivariateBin(override val center: RealVector, val sizes: RealVector, va
         return vector.asSequence().mapIndexed { i, value -> value in (center[i] - sizes[i] / 2)..(center[i] + sizes[i] / 2) }.all { it }
     }
 
-    override val value: Number get() = counter.sum()
+    override val value get() = counter.sum()
     internal operator fun inc() = this.also { counter.increment() }
 
     override val dimension: Int get() = center.size
@@ -51,8 +51,8 @@ class FastHistogram(
             val center = indexArray.mapIndexed { axis, index ->
                 when (index) {
                     0 -> Double.NEGATIVE_INFINITY
-                    actualSizes[axis] -> Double.POSITIVE_INFINITY
-                    else -> lower[axis] + (index - 1) * binSize[axis]
+                    actualSizes[axis] - 1 -> Double.POSITIVE_INFINITY
+                    else -> lower[axis] + (index.toDouble() - 0.5) * binSize[axis]
                 }
             }.toVector()
             MultivariateBin(center, binSize)
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/histogram/MultivariateHistogramTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/histogram/MultivariateHistogramTest.kt
index 2fbbbc88d..842a12ae4 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/histogram/MultivariateHistogramTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/histogram/MultivariateHistogramTest.kt
@@ -14,10 +14,11 @@ class MultivariateHistogramTest {
                 (-1.0..1.0),
                 (-1.0..1.0)
         )
-        histogram.put(0.6, 0.6)
+        histogram.put(0.55, 0.55)
         val bin = histogram.find { it.value.toInt() > 0 }!!
-        assertTrue { bin.contains(Vector.ofReal(0.6, 0.6)) }
-        assertFalse { bin.contains(Vector.ofReal(-0.6, 0.6)) }
+        assertTrue { bin.contains(Vector.ofReal(0.55, 0.55)) }
+        assertTrue { bin.contains(Vector.ofReal(0.6, 0.5)) }
+        assertFalse { bin.contains(Vector.ofReal(-0.55, 0.55)) }
     }
 
     @Test

From be18014d5400e5c1b19fb5f37b2214e19e521662 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Tue, 6 Nov 2018 21:22:03 +0300
Subject: [PATCH 02/13] RealNDField

---
 .../kmath/structures/LazyStructure.kt         | 20 ++++++++
 .../scientifik/kmath/structures/NDField.kt    | 10 ++--
 .../kmath/structures/RealNDField.kt           | 46 +++++++++++++++++++
 .../expressions/FieldExpressionContextTest.kt | 11 +++++
 .../kmath/structures/RealNDFieldTest.kt       | 16 +++++++
 5 files changed, 99 insertions(+), 4 deletions(-)
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LazyStructure.kt
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LazyStructure.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LazyStructure.kt
new file mode 100644
index 000000000..0428535f9
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LazyStructure.kt
@@ -0,0 +1,20 @@
+package scientifik.kmath.structures
+
+//
+//class LazyStructureField<T: Any>(val field: Field<T>): Field<LazyStructure<T>>{
+//
+//}
+//
+//class LazyStructure<T : Any> : NDStructure<T> {
+//
+//    override val shape: IntArray
+//        get() = TODO("not implemented") //To change initializer of created properties use File | Settings | File Templates.
+//
+//    override fun get(index: IntArray): T {
+//        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
+//    }
+//
+//    override fun iterator(): Iterator<Pair<IntArray, T>> {
+//        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
+//    }
+//}
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
index 9461bee85..7b991200b 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
@@ -1,6 +1,5 @@
 package scientifik.kmath.structures
 
-import scientifik.kmath.operations.DoubleField
 import scientifik.kmath.operations.Field
 import scientifik.kmath.operations.FieldElement
 
@@ -107,14 +106,15 @@ abstract class NDField<T>(val shape: IntArray, val field: Field<T>) : Field<NDAr
 /**
  *  Immutable [NDStructure] coupled to the context. Emulates Python ndarray
  */
-data class NDArray<T>(override val context: NDField<T>, private val structure: NDStructure<T>) : FieldElement<NDArray<T>, NDField<T>>, NDStructure<T> by structure {
+class NDArray<T>(override val context: NDField<T>, private val structure: NDStructure<T>) : FieldElement<NDArray<T>, NDField<T>>, NDStructure<T> by structure {
 
     //TODO ensure structure is immutable
 
     override val self: NDArray<T>
         get() = this
 
-    fun transform(action: (IntArray, T) -> T): NDArray<T> = context.produce { action(it, get(*it)) }
+    inline fun transform(crossinline action: (IntArray, T) -> T): NDArray<T> = context.produce { action(it, get(*it)) }
+    inline fun transform(crossinline action: (T) -> T): NDArray<T> = context.produce { action(get(*it)) }
 }
 
 /**
@@ -173,7 +173,7 @@ object NDArrays {
      * Create a platform-optimized NDArray of doubles
      */
     fun realNDArray(shape: IntArray, initializer: (IntArray) -> Double = { 0.0 }): NDArray<Double> {
-        return GenericNDField(shape, DoubleField).produce(initializer)
+        return RealNDField(shape).produce(initializer)
     }
 
     fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDArray<Double> {
@@ -188,6 +188,8 @@ object NDArrays {
         return realNDArray(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
     }
 
+    inline fun produceReal(shape: IntArray, block: RealNDField.() -> RealNDArray) = RealNDField(shape).run(block)
+
 //    /**
 //     * Simple boxing NDField
 //     */
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
new file mode 100644
index 000000000..2ba5fd706
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
@@ -0,0 +1,46 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.DoubleField
+import scientifik.kmath.operations.ExponentialOperations
+import scientifik.kmath.operations.PowerOperations
+import scientifik.kmath.operations.TrigonometricOperations
+import kotlin.math.*
+
+typealias RealNDArray = NDArray<Double>
+
+
+class RealNDField(shape: IntArray) : NDField<Double>(shape, DoubleField),
+        TrigonometricOperations<RealNDArray>,
+        PowerOperations<RealNDArray>,
+        ExponentialOperations<RealNDArray> {
+
+    override fun produceStructure(initializer: (IntArray) -> Double): NDStructure<Double> {
+        return genericNdStructure(shape, initializer)
+    }
+
+    override fun power(arg: RealNDArray, pow: Double): RealNDArray {
+        return arg.transform { d -> d.pow(pow) }
+    }
+
+    override fun exp(arg: RealNDArray): RealNDArray {
+        return arg.transform { d -> exp(d) }
+    }
+
+    override fun ln(arg: RealNDArray): RealNDArray {
+        return arg.transform { d -> ln(d) }
+    }
+
+    override fun sin(arg: RealNDArray): RealNDArray {
+        return arg.transform { d -> sin(d) }
+    }
+
+    override fun cos(arg: RealNDArray): RealNDArray {
+        return arg.transform { d -> cos(d) }
+    }
+
+    fun abs(arg: RealNDArray): RealNDArray {
+        return arg.transform { d -> abs(d) }
+    }
+}
+
+
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
index 543e13d79..8e2845b9e 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
@@ -39,4 +39,15 @@ class FieldExpressionContextTest {
         val expression = FieldExpressionContext(DoubleField).expression()
         assertEquals(expression("x" to 1.0), 4.0)
     }
+
+    @Test
+    fun valueExpression() {
+        val expressionBuilder: FieldExpressionContext<Double>.()->Expression<Double> = {
+            val x = variable("x")
+            x * x + 2 * x + 1.0
+        }
+
+        val expression = FieldExpressionContext(DoubleField).expressionBuilder()
+        assertEquals(expression("x" to 1.0), 4.0)
+    }
 }
\ No newline at end of file
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt
index 329f13189..179fb5ade 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt
@@ -1,6 +1,8 @@
 package scientifik.kmath.structures
 
+import scientifik.kmath.structures.NDArrays.produceReal
 import scientifik.kmath.structures.NDArrays.real2DArray
+import kotlin.math.abs
 import kotlin.math.pow
 import kotlin.test.Test
 import kotlin.test.assertEquals
@@ -40,4 +42,18 @@ class RealNDFieldTest {
         val result = function(array1) + 1.0
         assertEquals(10.0, result[1,1])
     }
+
+    @Test
+    fun testLibraryFunction() {
+        val abs: (Double) -> Double = ::abs
+        val result = abs(array1)
+        assertEquals(10.0, result[1,1])
+    }
+
+    @Test
+    fun testAbs(){
+        val res = produceReal(array1.shape){
+           1 + abs(array1) + exp(array2)
+        }
+    }
 }

From 68bd0ae0af3b7f57a50bac8ba3dfe7069801d2e7 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Wed, 14 Nov 2018 23:26:38 +0300
Subject: [PATCH 03/13] Norm is a separate context

---
 build.gradle                                  |  8 ++-
 kmath-core/build.gradle                       |  5 --
 .../scientifik/kmath/linear/LinearAlgrebra.kt | 16 +++---
 .../scientifik/kmath/operations/Fields.kt     | 49 +++++++++--------
 .../kmath/operations/OptionalOperations.kt    |  8 ++-
 .../scientifik/kmath/structures/Buffers.kt    |  3 +-
 .../kmath/structures/ExtendedNDField.kt       | 39 ++++++++++++++
 .../scientifik/kmath/structures/NDField.kt    |  7 +--
 .../kmath/structures/RealNDField.kt           | 46 ----------------
 .../kmath/operations/RealFieldTest.kt         |  4 +-
 ...ealNDFieldTest.kt => NumberNDFieldTest.kt} | 23 +++++---
 kmath-io/build.gradle                         | 53 +++++++++++++++++++
 kmath-jmh/build.gradle                        |  5 --
 settings.gradle                               |  3 ++
 14 files changed, 168 insertions(+), 101 deletions(-)
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
 rename kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/{RealNDFieldTest.kt => NumberNDFieldTest.kt} (68%)
 create mode 100644 kmath-io/build.gradle

diff --git a/build.gradle b/build.gradle
index 51d476ed6..4bd2464a9 100644
--- a/build.gradle
+++ b/build.gradle
@@ -1,7 +1,8 @@
 buildscript {
-    ext.kotlin_version = '1.3.0'
+    ext.kotlin_version = '1.3.10'
 
     repositories {
+        mavenCentral()
         jcenter()
     }
 
@@ -17,6 +18,11 @@ allprojects {
 
     group = 'scientifik'
     version = '0.0.1-SNAPSHOT'
+
+    repositories {
+        mavenCentral()
+        jcenter()
+    }
 }
 
 if(file('artifactory.gradle').exists()){
diff --git a/kmath-core/build.gradle b/kmath-core/build.gradle
index 18be067e6..c63e518fc 100644
--- a/kmath-core/build.gradle
+++ b/kmath-core/build.gradle
@@ -2,11 +2,6 @@ plugins {
     id 'kotlin-multiplatform'
 }
 
-repositories {
-    maven { url = 'http://dl.bintray.com/kotlin/kotlin-eap' }
-    mavenCentral()
-}
-
 kotlin {
     targets {
         fromPreset(presets.jvm, 'jvm')
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
index e585b884d..39931b676 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
@@ -1,9 +1,6 @@
 package scientifik.kmath.linear
 
-import scientifik.kmath.operations.DoubleField
-import scientifik.kmath.operations.Field
-import scientifik.kmath.operations.Space
-import scientifik.kmath.operations.SpaceElement
+import scientifik.kmath.operations.*
 import scientifik.kmath.structures.*
 
 /**
@@ -70,8 +67,6 @@ abstract class MatrixSpace<T : Any>(val rows: Int, val columns: Int, val field:
         result = 31 * result + field.hashCode()
         return result
     }
-
-
 }
 
 infix fun <T : Any> Matrix<T>.dot(b: Matrix<T>): Matrix<T> = this.context.multiply(this, b)
@@ -191,11 +186,11 @@ interface Vector<T : Any> : SpaceElement<Vector<T>, VectorSpace<T>>, Buffer<T>,
 typealias NDFieldFactory<T> = (IntArray) -> NDField<T>
 
 internal fun <T : Any> genericNDFieldFactory(field: Field<T>): NDFieldFactory<T> = { index -> GenericNDField(index, field) }
-internal val realNDFieldFactory: NDFieldFactory<Double> = { index -> GenericNDField(index, DoubleField) }
+internal val realNDFieldFactory: NDFieldFactory<Double> = { index -> ExtendedNDField(index, DoubleField) }
 
 
 /**
- * NDArray-based implementation of vector space. By default uses slow [SimpleNDField], but could be overridden with custom [NDField] factory.
+ * NDArray-based implementation of vector space. By default uses slow [GenericNDField], but could be overridden with custom [NDField] factory.
  */
 class ArrayMatrixSpace<T : Any>(
         rows: Int,
@@ -318,3 +313,8 @@ fun <T : Any> Vector<T>.toMatrix(): Matrix<T> {
     return Matrix.of(size, 1, context.field) { i, _ -> get(i) }
 }
 
+object VectorL2Norm: Norm<Vector<out Number>, Double> {
+    override fun norm(arg: Vector<out Number>): Double {
+        return kotlin.math.sqrt(arg.sumByDouble { it.toDouble() })
+    }
+}
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
index eeba76222..5702803fc 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
@@ -2,10 +2,20 @@ package scientifik.kmath.operations
 
 import kotlin.math.pow
 
+/**
+ * Advanced Number-like field that implements basic operations
+ */
+interface ExtendedField<N : Any> :
+        Field<N>,
+        TrigonometricOperations<N>,
+        PowerOperations<N>,
+        ExponentialOperations<N>
+
+
 /**
  * Field for real values
  */
-object RealField : Field<Real>, TrigonometricOperations<Real>, PowerOperations<Real>, ExponentialOperations<Real> {
+object RealField : ExtendedField<Real>, Norm<Real, Real> {
     override val zero: Real = Real(0.0)
     override fun add(a: Real, b: Real): Real = Real(a.value + b.value)
     override val one: Real = Real(1.0)
@@ -13,33 +23,24 @@ object RealField : Field<Real>, TrigonometricOperations<Real>, PowerOperations<R
     override fun multiply(a: Real, k: Double): Real = Real(a.value * k)
     override fun divide(a: Real, b: Real): Real = Real(a.value / b.value)
 
-    override fun sin(arg: Real): Real  = Real(kotlin.math.sin(arg.value))
+    override fun sin(arg: Real): Real = Real(kotlin.math.sin(arg.value))
     override fun cos(arg: Real): Real = Real(kotlin.math.cos(arg.value))
 
-    override fun power(arg: Real, pow: Double): Real  = Real(arg.value.pow(pow))
+    override fun power(arg: Real, pow: Double): Real = Real(arg.value.pow(pow))
 
-    override fun exp(arg: Real): Real  = Real(kotlin.math.exp(arg.value))
+    override fun exp(arg: Real): Real = Real(kotlin.math.exp(arg.value))
 
-    override fun ln(arg: Real): Real  = Real(kotlin.math.ln(arg.value))
+    override fun ln(arg: Real): Real = Real(kotlin.math.ln(arg.value))
+
+    override fun norm(arg: Real): Real = Real(kotlin.math.abs(arg.value))
 }
 
 /**
  * Real field element wrapping double.
  *
- * TODO could be replaced by inline class in kotlin 1.3 if it would allow to avoid boxing
+ * TODO inline does not work due to compiler bug. Waiting for fix
  */
-data class Real(val value: Double) : Number(), FieldElement<Real, RealField> {
-    /*
-     * The class uses composition instead of inheritance since Double is final
-     */
-
-    override fun toByte(): Byte = value.toByte()
-    override fun toChar(): Char = value.toChar()
-    override fun toDouble(): Double = value
-    override fun toFloat(): Float = value.toFloat()
-    override fun toInt(): Int = value.toInt()
-    override fun toLong(): Long = value.toLong()
-    override fun toShort(): Short = value.toShort()
+class Real(val value: Double) : FieldElement<Real, RealField> {
 
     //values are dynamically calculated to save memory
     override val self
@@ -53,19 +54,21 @@ data class Real(val value: Double) : Number(), FieldElement<Real, RealField> {
 /**
  * A field for double without boxing. Does not produce appropriate field element
  */
-object DoubleField : Field<Double>, TrigonometricOperations<Double>, PowerOperations<Double>, ExponentialOperations<Double> {
+object DoubleField : ExtendedField<Double>, Norm<Double, Double> {
     override val zero: Double = 0.0
     override fun add(a: Double, b: Double): Double = a + b
     override fun multiply(a: Double, @Suppress("PARAMETER_NAME_CHANGED_ON_OVERRIDE") b: Double): Double = a * b
     override val one: Double = 1.0
     override fun divide(a: Double, b: Double): Double = a / b
 
-    override fun sin(arg: Double): Double  = kotlin.math.sin(arg)
+    override fun sin(arg: Double): Double = kotlin.math.sin(arg)
     override fun cos(arg: Double): Double = kotlin.math.cos(arg)
 
-    override fun power(arg: Double, pow: Double): Double  = arg.pow(pow)
+    override fun power(arg: Double, pow: Double): Double = arg.pow(pow)
 
-    override fun exp(arg: Double): Double  =kotlin.math.exp(arg)
+    override fun exp(arg: Double): Double = kotlin.math.exp(arg)
 
-    override fun ln(arg: Double): Double  = kotlin.math.ln(arg)
+    override fun ln(arg: Double): Double = kotlin.math.ln(arg)
+
+    override fun norm(arg: Double): Double = kotlin.math.abs(arg)
 }
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt
index 7628cb4fc..9e5c8a801 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt
@@ -45,4 +45,10 @@ interface ExponentialOperations<T> {
 }
 
 fun <T : MathElement<T, out ExponentialOperations<T>>> exp(arg: T): T = arg.context.exp(arg)
-fun <T : MathElement<T, out ExponentialOperations<T>>> ln(arg: T): T = arg.context.ln(arg)
\ No newline at end of file
+fun <T : MathElement<T, out ExponentialOperations<T>>> ln(arg: T): T = arg.context.ln(arg)
+
+interface Norm<in T, out R> {
+    fun norm(arg: T): R
+}
+
+fun <T : MathElement<T, out Norm<T, R>>, R> norm(arg: T): R = arg.context.norm(arg)
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
index 04983ef85..1be3514f2 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
@@ -42,6 +42,7 @@ inline class ListBuffer<T>(private val list: MutableList<T>) : MutableBuffer<T>
 }
 
 class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
+    //Can't inline because array invariant
     override val size: Int
         get() = array.size
 
@@ -56,7 +57,7 @@ class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
     override fun copy(): MutableBuffer<T> = ArrayBuffer(array.copyOf())
 }
 
-class DoubleBuffer(private val array: DoubleArray) : MutableBuffer<Double> {
+inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer<Double> {
     override val size: Int
         get() = array.size
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt
new file mode 100644
index 000000000..6392aeec6
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt
@@ -0,0 +1,39 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.*
+
+
+/**
+ * NDField that supports [ExtendedField] operations on its elements
+ */
+class ExtendedNDField<N: Any>(shape: IntArray, override val field: ExtendedField<N>) : NDField<N>(shape, field),
+        TrigonometricOperations<NDArray<N>>,
+        PowerOperations<NDArray<N>>,
+        ExponentialOperations<NDArray<N>> {
+
+    override fun produceStructure(initializer: (IntArray) -> N): NDStructure<N> {
+        return genericNdStructure(shape, initializer)
+    }
+
+    override fun power(arg: NDArray<N>, pow: Double): NDArray<N> {
+        return arg.transform { d -> with(field){power(d,pow)} }
+    }
+
+    override fun exp(arg: NDArray<N>): NDArray<N> {
+        return arg.transform { d -> with(field){exp(d)} }
+    }
+
+    override fun ln(arg: NDArray<N>): NDArray<N> {
+        return arg.transform { d -> with(field){ln(d)} }
+    }
+
+    override fun sin(arg: NDArray<N>): NDArray<N> {
+        return arg.transform { d -> with(field){sin(d)} }
+    }
+
+    override fun cos(arg: NDArray<N>): NDArray<N> {
+        return arg.transform { d -> with(field){cos(d)} }
+    }
+}
+
+
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
index 7b991200b..74bffa7ad 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
@@ -1,5 +1,6 @@
 package scientifik.kmath.structures
 
+import scientifik.kmath.operations.DoubleField
 import scientifik.kmath.operations.Field
 import scientifik.kmath.operations.FieldElement
 
@@ -14,7 +15,7 @@ class ShapeMismatchException(val expected: IntArray, val actual: IntArray) : Run
  * @param field - operations field defined on individual array element
  * @param T the type of the element contained in NDArray
  */
-abstract class NDField<T>(val shape: IntArray, val field: Field<T>) : Field<NDArray<T>> {
+abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field<NDArray<T>> {
 
     abstract fun produceStructure(initializer: (IntArray) -> T): NDStructure<T>
 
@@ -173,7 +174,7 @@ object NDArrays {
      * Create a platform-optimized NDArray of doubles
      */
     fun realNDArray(shape: IntArray, initializer: (IntArray) -> Double = { 0.0 }): NDArray<Double> {
-        return RealNDField(shape).produce(initializer)
+        return ExtendedNDField(shape, DoubleField).produce(initializer)
     }
 
     fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDArray<Double> {
@@ -188,7 +189,7 @@ object NDArrays {
         return realNDArray(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
     }
 
-    inline fun produceReal(shape: IntArray, block: RealNDField.() -> RealNDArray) = RealNDField(shape).run(block)
+    inline fun produceReal(shape: IntArray, block: ExtendedNDField<Double>.() -> NDArray<Double>) = ExtendedNDField(shape, DoubleField).run(block)
 
 //    /**
 //     * Simple boxing NDField
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
deleted file mode 100644
index 2ba5fd706..000000000
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/RealNDField.kt
+++ /dev/null
@@ -1,46 +0,0 @@
-package scientifik.kmath.structures
-
-import scientifik.kmath.operations.DoubleField
-import scientifik.kmath.operations.ExponentialOperations
-import scientifik.kmath.operations.PowerOperations
-import scientifik.kmath.operations.TrigonometricOperations
-import kotlin.math.*
-
-typealias RealNDArray = NDArray<Double>
-
-
-class RealNDField(shape: IntArray) : NDField<Double>(shape, DoubleField),
-        TrigonometricOperations<RealNDArray>,
-        PowerOperations<RealNDArray>,
-        ExponentialOperations<RealNDArray> {
-
-    override fun produceStructure(initializer: (IntArray) -> Double): NDStructure<Double> {
-        return genericNdStructure(shape, initializer)
-    }
-
-    override fun power(arg: RealNDArray, pow: Double): RealNDArray {
-        return arg.transform { d -> d.pow(pow) }
-    }
-
-    override fun exp(arg: RealNDArray): RealNDArray {
-        return arg.transform { d -> exp(d) }
-    }
-
-    override fun ln(arg: RealNDArray): RealNDArray {
-        return arg.transform { d -> ln(d) }
-    }
-
-    override fun sin(arg: RealNDArray): RealNDArray {
-        return arg.transform { d -> sin(d) }
-    }
-
-    override fun cos(arg: RealNDArray): RealNDArray {
-        return arg.transform { d -> cos(d) }
-    }
-
-    fun abs(arg: RealNDArray): RealNDArray {
-        return arg.transform { d -> abs(d) }
-    }
-}
-
-
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt
index b07158cc3..551bf0ee4 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt
@@ -7,8 +7,8 @@ class RealFieldTest {
     @Test
     fun testSqrt() {
         val sqrt = with(RealField) {
-            sqrt(25 * one)
+            sqrt(  25 * one)
         }
-        assertEquals(5.0, sqrt.toDouble())
+        assertEquals(5.0, sqrt.value)
     }
 }
\ No newline at end of file
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
similarity index 68%
rename from kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt
rename to kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
index 179fb5ade..844535fa8 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
@@ -1,5 +1,8 @@
 package scientifik.kmath.structures
 
+import scientifik.kmath.linear.Vector
+import scientifik.kmath.linear.VectorL2Norm
+import scientifik.kmath.operations.Norm
 import scientifik.kmath.structures.NDArrays.produceReal
 import scientifik.kmath.structures.NDArrays.real2DArray
 import kotlin.math.abs
@@ -7,7 +10,7 @@ import kotlin.math.pow
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
-class RealNDFieldTest {
+class NumberNDFieldTest {
     val array1 = real2DArray(3, 3) { i, j -> (i + j).toDouble() }
     val array2 = real2DArray(3, 3) { i, j -> (i - j).toDouble() }
 
@@ -46,14 +49,22 @@ class RealNDFieldTest {
     @Test
     fun testLibraryFunction() {
         val abs: (Double) -> Double = ::abs
-        val result = abs(array1)
-        assertEquals(10.0, result[1,1])
+        val result = abs(array2)
+        assertEquals(2.0, result[0,2])
+    }
+
+    object L2Norm: Norm<NDArray<out Number>, Double> {
+        override fun norm(arg: NDArray<out Number>): Double {
+            return kotlin.math.sqrt(arg.sumByDouble { it.second.toDouble() })
+        }
     }
 
     @Test
-    fun testAbs(){
-        val res = produceReal(array1.shape){
-           1 + abs(array1) + exp(array2)
+    fun testInternalContext(){
+        produceReal(array1.shape){
+            with(L2Norm) {
+                1 + norm(array1) + exp(array2)
+            }
         }
     }
 }
diff --git a/kmath-io/build.gradle b/kmath-io/build.gradle
new file mode 100644
index 000000000..d0d93262e
--- /dev/null
+++ b/kmath-io/build.gradle
@@ -0,0 +1,53 @@
+plugins {
+    id 'kotlin-multiplatform'
+}
+
+kotlin {
+    targets {
+        fromPreset(presets.jvm, 'jvm')
+        fromPreset(presets.js, 'js')
+        // For ARM, preset should be changed to presets.iosArm32 or presets.iosArm64
+        // For Linux, preset should be changed to e.g. presets.linuxX64
+        // For MacOS, preset should be changed to e.g. presets.macosX64
+        //fromPreset(presets.mingwX64, 'mingw')
+    }
+    sourceSets {
+        commonMain {
+            dependencies {
+                api project(":kmath-core")
+                implementation 'org.jetbrains.kotlin:kotlin-stdlib-common'
+            }
+        }
+        commonTest {
+            dependencies {
+                implementation 'org.jetbrains.kotlin:kotlin-test-common'
+                implementation 'org.jetbrains.kotlin:kotlin-test-annotations-common'
+            }
+        }
+        jvmMain {
+            dependencies {
+                implementation 'org.jetbrains.kotlin:kotlin-stdlib-jdk8'
+            }
+        }
+        jvmTest {
+            dependencies {
+                implementation 'org.jetbrains.kotlin:kotlin-test'
+                implementation 'org.jetbrains.kotlin:kotlin-test-junit'
+            }
+        }
+        jsMain {
+            dependencies {
+                implementation 'org.jetbrains.kotlin:kotlin-stdlib-js'
+            }
+        }
+        jsTest {
+            dependencies {
+                implementation 'org.jetbrains.kotlin:kotlin-test-js'
+            }
+        }
+//        mingwMain {
+//        }
+//        mingwTest {
+//        }
+    }
+}
diff --git a/kmath-jmh/build.gradle b/kmath-jmh/build.gradle
index a16bf3022..97337cb40 100644
--- a/kmath-jmh/build.gradle
+++ b/kmath-jmh/build.gradle
@@ -4,11 +4,6 @@ plugins {
     id "me.champeau.gradle.jmh" version "0.4.7"
 }
 
-repositories {
-    maven { url = 'http://dl.bintray.com/kotlin/kotlin-eap' }
-    mavenCentral()
-}
-
 dependencies {
     implementation project(':kmath-core')
     jmh 'org.jetbrains.kotlin:kotlin-stdlib-jdk8'
diff --git a/settings.gradle b/settings.gradle
index df4ccc99e..ddcbef5dd 100644
--- a/settings.gradle
+++ b/settings.gradle
@@ -1,6 +1,7 @@
 pluginManagement {
     repositories {
         mavenCentral()
+        jcenter()
         maven { url = 'https://plugins.gradle.org/m2/' }
     }
 }
@@ -9,5 +10,7 @@ enableFeaturePreview('GRADLE_METADATA')
 
 rootProject.name = 'kmath'
 include ':kmath-core'
+include ':kmath-io'
+
 include ':kmath-jmh'
 

From dcf8f4c6fd38d7e4c72d5dfe09b8556a43bab9f6 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Thu, 15 Nov 2018 11:07:51 +0300
Subject: [PATCH 04/13] inline class bug

---
 .../commonMain/kotlin/scientifik/kmath/operations/Fields.kt   | 4 ++--
 .../kotlin/scientifik/kmath/operations/RealFieldTest.kt       | 2 ++
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
index 5702803fc..3e4f47869 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
@@ -38,9 +38,9 @@ object RealField : ExtendedField<Real>, Norm<Real, Real> {
 /**
  * Real field element wrapping double.
  *
- * TODO inline does not work due to compiler bug. Waiting for fix
+ * TODO inline does not work due to compiler bug. Waiting for fix for KT-27586
  */
-class Real(val value: Double) : FieldElement<Real, RealField> {
+inline class Real(val value: Double) : FieldElement<Real, RealField> {
 
     //values are dynamically calculated to save memory
     override val self
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt
index 551bf0ee4..0d33204df 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt
@@ -6,6 +6,8 @@ import kotlin.test.assertEquals
 class RealFieldTest {
     @Test
     fun testSqrt() {
+
+        //fails because KT-27586
         val sqrt = with(RealField) {
             sqrt(  25 * one)
         }

From 62c7099f8aea8882c3c5c69568f62f4b50c1f400 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Tue, 20 Nov 2018 22:09:07 +0300
Subject: [PATCH 05/13] Minor change to histogram API. Project build to kts.

---
 build.gradle                                  | 24 -----------
 build.gradle.kts                              | 28 +++++++++++++
 .../kmath/histogram/FastHistogram.kt          |  8 ++--
 .../scientifik/kmath/histogram/Histogram.kt   | 20 +++++----
 .../kmath/histogram/PhantomHistogram.kt       | 42 +++++++++++++++++++
 .../scientifik/kmath/structures/Buffers.kt    | 17 +++++++-
 .../kmath/histogram/UnivariateHistogram.kt    |  7 +++-
 settings.gradle                               | 13 ------
 settings.gradle.kts                           | 14 +++++++
 9 files changed, 123 insertions(+), 50 deletions(-)
 delete mode 100644 build.gradle
 create mode 100644 build.gradle.kts
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
 delete mode 100644 settings.gradle
 create mode 100644 settings.gradle.kts

diff --git a/build.gradle b/build.gradle
deleted file mode 100644
index 51d476ed6..000000000
--- a/build.gradle
+++ /dev/null
@@ -1,24 +0,0 @@
-buildscript {
-    ext.kotlin_version = '1.3.0'
-
-    repositories {
-        jcenter()
-    }
-
-    dependencies {
-        classpath "org.jetbrains.kotlin:kotlin-gradle-plugin:$kotlin_version"
-        classpath "org.jfrog.buildinfo:build-info-extractor-gradle:4+"
-    }
-}
-
-allprojects {
-    apply plugin: 'maven-publish'
-    apply plugin: "com.jfrog.artifactory"
-
-    group = 'scientifik'
-    version = '0.0.1-SNAPSHOT'
-}
-
-if(file('artifactory.gradle').exists()){
-    apply from: 'artifactory.gradle'
-}
\ No newline at end of file
diff --git a/build.gradle.kts b/build.gradle.kts
new file mode 100644
index 000000000..2d7b9f119
--- /dev/null
+++ b/build.gradle.kts
@@ -0,0 +1,28 @@
+buildscript {
+    val kotlin_version = "1.3.10"
+
+    repositories {
+        jcenter()
+    }
+
+    dependencies {
+        classpath("org.jetbrains.kotlin:kotlin-gradle-plugin:$kotlin_version")
+        classpath("org.jfrog.buildinfo:build-info-extractor-gradle:4+")
+    }
+}
+
+plugins {
+    id("com.jfrog.artifactory") version "4.8.1"  apply false
+}
+
+allprojects {
+    apply(plugin = "maven-publish")
+    apply(plugin = "com.jfrog.artifactory")
+
+    group = "scientifik"
+    version = "0.0.1-SNAPSHOT"
+}
+
+if(file("artifactory.gradle").exists()){
+    apply(from = "artifactory.gradle")
+}
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
index 3852de316..621ea90c0 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
@@ -24,13 +24,14 @@ class MultivariateBin(override val center: RealVector, val sizes: RealVector, va
 }
 
 /**
- * Uniform multivariate histogram with fixed borders. Based on NDStructure implementation with complexity of m for bin search, where m is the number of dimensions
+ * Uniform multivariate histogram with fixed borders. Based on NDStructure implementation with complexity of m for bin search, where m is the number of dimensions.
+ * The histogram is optimized for speed, but have large size in memory
  */
 class FastHistogram(
         private val lower: RealVector,
         private val upper: RealVector,
         private val binNums: IntArray = IntArray(lower.size) { 20 }
-) : Histogram<Double, MultivariateBin> {
+) : MutableHistogram<Double, MultivariateBin> {
 
     init {
         // argument checks
@@ -76,7 +77,8 @@ class FastHistogram(
         return bins[index]
     }
 
-    override fun put(point: Buffer<out Double>) {
+    override fun put(point: Buffer<out Double>, weight: Double) {
+        if (weight != 1.0) TODO("Implement weighting")
         this[point]?.inc() ?: error("Could not find appropriate bin (should not be possible)")
     }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
index 1b5c231b4..e143bd90d 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
@@ -4,12 +4,14 @@ import scientifik.kmath.operations.Space
 import scientifik.kmath.structures.ArrayBuffer
 import scientifik.kmath.structures.Buffer
 
+typealias Point<T> = Buffer<T>
+
 /**
  * A simple geometric domain
  * TODO move to geometry module
  */
 interface Domain<T: Any> {
-    operator fun contains(vector: Buffer<out T>): Boolean
+    operator fun contains(vector: Point<out T>): Boolean
     val dimension: Int
 }
 
@@ -21,7 +23,7 @@ interface Bin<T: Any> : Domain<T> {
      * The value of this bin
      */
     val value: Number
-    val center: Buffer<T>
+    val center: Point<T>
 }
 
 interface Histogram<T: Any, out B : Bin<T>> : Iterable<B> {
@@ -29,27 +31,31 @@ interface Histogram<T: Any, out B : Bin<T>> : Iterable<B> {
     /**
      * Find existing bin, corresponding to given coordinates
      */
-    operator fun get(point: Buffer<out T>): B?
+    operator fun get(point: Point<out T>): B?
 
     /**
      * Dimension of the histogram
      */
     val dimension: Int
 
+}
+
+interface MutableHistogram<T: Any, out B : Bin<T>>: Histogram<T,B>{
+
     /**
      * Increment appropriate bin
      */
-    fun put(point: Buffer<out T>)
+    fun put(point: Point<out T>, weight: Double = 1.0)
 }
 
-fun <T: Any> Histogram<T,*>.put(vararg point: T) = put(ArrayBuffer(point))
+fun <T: Any> MutableHistogram<T,*>.put(vararg point: T) = put(ArrayBuffer(point))
 
-fun <T: Any> Histogram<T,*>.fill(sequence: Iterable<Buffer<T>>) = sequence.forEach { put(it) }
+fun <T: Any> MutableHistogram<T,*>.fill(sequence: Iterable<Point<T>>) = sequence.forEach { put(it) }
 
 /**
  * Pass a sequence builder into histogram
  */
-fun <T: Any> Histogram<T, *>.fill(buider: suspend SequenceScope<Buffer<T>>.() -> Unit) = fill(sequence(buider).asIterable())
+fun <T: Any> MutableHistogram<T, *>.fill(buider: suspend SequenceScope<Point<T>>.() -> Unit) = fill(sequence(buider).asIterable())
 
 /**
  * A space to perform arithmetic operations on histograms
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
new file mode 100644
index 000000000..469856daa
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
@@ -0,0 +1,42 @@
+package scientifik.kmath.histogram
+
+import scientifik.kmath.linear.RealVector
+import scientifik.kmath.structures.NDStructure
+
+class BinTemplate(val center: RealVector, val sizes: RealVector) {
+    fun contains(vector: Point<out Double>): Boolean {
+        if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}")
+        return vector.asSequence().mapIndexed { i, value -> value in (center[i] - sizes[i] / 2)..(center[i] + sizes[i] / 2) }.all { it }
+    }
+}
+
+class PhantomBin(val template: BinTemplate, override val value: Number) : Bin<Double> {
+
+    override fun contains(vector: Point<out Double>): Boolean = template.contains(vector)
+
+    override val dimension: Int
+        get() = template.center.size
+
+    override val center: Point<Double>
+        get() = template.center
+
+}
+
+class PhantomHistogram(
+        val bins: Map<BinTemplate, IntArray>,
+        val data: NDStructure<Double>
+) : Histogram<Double, PhantomBin> {
+
+    override val dimension: Int
+        get() = data.dimension
+
+    override fun iterator(): Iterator<PhantomBin> {
+        return bins.asSequence().map {entry-> PhantomBin(entry.key,data[entry.value]) }.iterator()
+    }
+
+    override fun get(point: Point<out Double>): PhantomBin? {
+        val template = bins.keys.find { it.contains(point) }
+        return template?.let { PhantomBin(it, data[bins[it]!!]) }
+    }
+
+}
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
index 04983ef85..3e2712f0c 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
@@ -56,7 +56,7 @@ class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
     override fun copy(): MutableBuffer<T> = ArrayBuffer(array.copyOf())
 }
 
-class DoubleBuffer(private val array: DoubleArray) : MutableBuffer<Double> {
+inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer<Double> {
     override val size: Int
         get() = array.size
 
@@ -71,6 +71,21 @@ class DoubleBuffer(private val array: DoubleArray) : MutableBuffer<Double> {
     override fun copy(): MutableBuffer<Double> = DoubleBuffer(array.copyOf())
 }
 
+inline class IntBuffer(private val array: IntArray): MutableBuffer<Int>{
+    override val size: Int
+        get() = array.size
+
+    override fun get(index: Int): Int = array[index]
+
+    override fun set(index: Int, value: Int) {
+        array[index] = value
+    }
+
+    override fun iterator(): Iterator<Int>  = array.iterator()
+
+    override fun copy(): MutableBuffer<Int> = IntBuffer(array.copyOf())
+}
+
 inline fun <reified T : Any> buffer(size: Int, noinline initializer: (Int) -> T): Buffer<T> {
     return ArrayBuffer(Array(size, initializer))
 }
diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt
index 2382c379c..26ad520e6 100644
--- a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt
+++ b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt
@@ -26,7 +26,7 @@ class UnivariateBin(val position: Double, val size: Double, val counter: LongCou
 /**
  * Univariate histogram with log(n) bin search speed
  */
-class UnivariateHistogram private constructor(private val factory: (Double) -> UnivariateBin) : Histogram<Double,UnivariateBin> {
+class UnivariateHistogram private constructor(private val factory: (Double) -> UnivariateBin) : MutableHistogram<Double,UnivariateBin> {
 
     private val bins: TreeMap<Double, UnivariateBin> = TreeMap()
 
@@ -58,7 +58,10 @@ class UnivariateHistogram private constructor(private val factory: (Double) -> U
         (get(value) ?: createBin(value)).inc()
     }
 
-    override fun put(point: Buffer<out Double>) = put(point[0])
+    override fun put(point: Buffer<out Double>, weight: Double) {
+        if (weight != 1.0) TODO("Implement weighting")
+        put(point[0])
+    }
 
     companion object {
         fun uniform(binSize: Double, start: Double = 0.0): UnivariateHistogram {
diff --git a/settings.gradle b/settings.gradle
deleted file mode 100644
index df4ccc99e..000000000
--- a/settings.gradle
+++ /dev/null
@@ -1,13 +0,0 @@
-pluginManagement {
-    repositories {
-        mavenCentral()
-        maven { url = 'https://plugins.gradle.org/m2/' }
-    }
-}
-
-enableFeaturePreview('GRADLE_METADATA')
-
-rootProject.name = 'kmath'
-include ':kmath-core'
-include ':kmath-jmh'
-
diff --git a/settings.gradle.kts b/settings.gradle.kts
new file mode 100644
index 000000000..0f90d4c0b
--- /dev/null
+++ b/settings.gradle.kts
@@ -0,0 +1,14 @@
+pluginManagement {
+    repositories {
+        mavenCentral()
+        maven("https://plugins.gradle.org/m2/")
+    }
+}
+
+enableFeaturePreview("GRADLE_METADATA")
+
+rootProject.name = "kmath"
+include(
+        ":kmath-core",
+        ":kmath-jmh"
+)

From 5cd6301f452c0605fa84b671a0bf55655a52bdea Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Wed, 21 Nov 2018 09:55:39 +0300
Subject: [PATCH 06/13] NDArray renamed to NDElement

---
 build.gradle.kts                              |  6 +-
 .../kmath/histogram/PhantomHistogram.kt       |  4 ++
 .../scientifik/kmath/linear/LinearAlgrebra.kt | 14 ++--
 .../kmath/structures/ExtendedNDField.kt       | 21 +++---
 .../scientifik/kmath/structures/NDField.kt    | 66 +++++++++----------
 .../kmath/structures/NumberNDFieldTest.kt     |  6 +-
 6 files changed, 62 insertions(+), 55 deletions(-)

diff --git a/build.gradle.kts b/build.gradle.kts
index 2d7b9f119..8f8bcaa34 100644
--- a/build.gradle.kts
+++ b/build.gradle.kts
@@ -1,12 +1,14 @@
 buildscript {
-    val kotlin_version = "1.3.10"
+    extra["kotlinVersion"] = "1.3.10"
+
+    val kotlinVersion: String by extra
 
     repositories {
         jcenter()
     }
 
     dependencies {
-        classpath("org.jetbrains.kotlin:kotlin-gradle-plugin:$kotlin_version")
+        classpath("org.jetbrains.kotlin:kotlin-gradle-plugin:$kotlinVersion")
         classpath("org.jfrog.buildinfo:build-info-extractor-gradle:4+")
     }
 }
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
index 469856daa..0bd5de0e2 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
@@ -22,6 +22,10 @@ class PhantomBin(val template: BinTemplate, override val value: Number) : Bin<Do
 
 }
 
+/**
+ * Immutable histogram with explicit structure for content and additional external bin description.
+ * Bin search is slow, but full histogram algebra is supported.
+ */
 class PhantomHistogram(
         val bins: Map<BinTemplate, IntArray>,
         val data: NDStructure<Double>
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
index 39931b676..62852ff5d 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
@@ -225,7 +225,7 @@ class ArrayVectorSpace<T : Any>(
 /**
  * Member of [ArrayMatrixSpace] which wraps 2-D array
  */
-class ArrayMatrix<T : Any> internal constructor(override val context: ArrayMatrixSpace<T>, val array: NDArray<T>) : Matrix<T> {
+class ArrayMatrix<T : Any> internal constructor(override val context: ArrayMatrixSpace<T>, val element: NDElement<T>) : Matrix<T> {
 
     constructor(context: ArrayMatrixSpace<T>, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) })
 
@@ -234,32 +234,32 @@ class ArrayMatrix<T : Any> internal constructor(override val context: ArrayMatri
     override val columns: Int get() = context.columns
 
     override fun get(i: Int, j: Int): T {
-        return array[i, j]
+        return element[i, j]
     }
 
     override val self: ArrayMatrix<T> get() = this
 }
 
 
-class ArrayVector<T : Any> internal constructor(override val context: ArrayVectorSpace<T>, val array: NDArray<T>) : Vector<T> {
+class ArrayVector<T : Any> internal constructor(override val context: ArrayVectorSpace<T>, val element: NDElement<T>) : Vector<T> {
 
     constructor(context: ArrayVectorSpace<T>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
 
     init {
-        if (context.size != array.shape[0]) {
+        if (context.size != element.shape[0]) {
             error("Array dimension mismatch")
         }
     }
 
     override fun get(index: Int): T {
-        return array[index]
+        return element[index]
     }
 
     override val self: ArrayVector<T> get() = this
 
-    override fun iterator(): Iterator<T> = (0 until size).map { array[it] }.iterator()
+    override fun iterator(): Iterator<T> = (0 until size).map { element[it] }.iterator()
 
-    override fun copy(): ArrayVector<T> = ArrayVector(context, array)
+    override fun copy(): ArrayVector<T> = ArrayVector(context, element)
 
     override fun toString(): String = this.joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() }
 }
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt
index 6392aeec6..08751d934 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt
@@ -1,37 +1,40 @@
 package scientifik.kmath.structures
 
-import scientifik.kmath.operations.*
+import scientifik.kmath.operations.ExponentialOperations
+import scientifik.kmath.operations.ExtendedField
+import scientifik.kmath.operations.PowerOperations
+import scientifik.kmath.operations.TrigonometricOperations
 
 
 /**
  * NDField that supports [ExtendedField] operations on its elements
  */
 class ExtendedNDField<N: Any>(shape: IntArray, override val field: ExtendedField<N>) : NDField<N>(shape, field),
-        TrigonometricOperations<NDArray<N>>,
-        PowerOperations<NDArray<N>>,
-        ExponentialOperations<NDArray<N>> {
+        TrigonometricOperations<NDElement<N>>,
+        PowerOperations<NDElement<N>>,
+        ExponentialOperations<NDElement<N>> {
 
     override fun produceStructure(initializer: (IntArray) -> N): NDStructure<N> {
         return genericNdStructure(shape, initializer)
     }
 
-    override fun power(arg: NDArray<N>, pow: Double): NDArray<N> {
+    override fun power(arg: NDElement<N>, pow: Double): NDElement<N> {
         return arg.transform { d -> with(field){power(d,pow)} }
     }
 
-    override fun exp(arg: NDArray<N>): NDArray<N> {
+    override fun exp(arg: NDElement<N>): NDElement<N> {
         return arg.transform { d -> with(field){exp(d)} }
     }
 
-    override fun ln(arg: NDArray<N>): NDArray<N> {
+    override fun ln(arg: NDElement<N>): NDElement<N> {
         return arg.transform { d -> with(field){ln(d)} }
     }
 
-    override fun sin(arg: NDArray<N>): NDArray<N> {
+    override fun sin(arg: NDElement<N>): NDElement<N> {
         return arg.transform { d -> with(field){sin(d)} }
     }
 
-    override fun cos(arg: NDArray<N>): NDArray<N> {
+    override fun cos(arg: NDElement<N>): NDElement<N> {
         return arg.transform { d -> with(field){cos(d)} }
     }
 }
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
index 74bffa7ad..6d3ce05a5 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
@@ -15,7 +15,7 @@ class ShapeMismatchException(val expected: IntArray, val actual: IntArray) : Run
  * @param field - operations field defined on individual array element
  * @param T the type of the element contained in NDArray
  */
-abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field<NDArray<T>> {
+abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field<NDElement<T>> {
 
     abstract fun produceStructure(initializer: (IntArray) -> T): NDStructure<T>
 
@@ -23,17 +23,17 @@ abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field
      * Create new instance of NDArray using field shape and given initializer
      * The producer takes list of indices as argument and returns contained value
      */
-    fun produce(initializer: (IntArray) -> T): NDArray<T> = NDArray(this, produceStructure(initializer))
+    fun produce(initializer: (IntArray) -> T): NDElement<T> = NDElement(this, produceStructure(initializer))
 
-    override val zero: NDArray<T> by lazy {
+    override val zero: NDElement<T> by lazy {
         produce { this.field.zero }
     }
 
     /**
      * Check the shape of given NDArray and throw exception if it does not coincide with shape of the field
      */
-    private fun checkShape(vararg arrays: NDArray<T>) {
-        arrays.forEach {
+    private fun checkShape(vararg elements: NDElement<T>) {
+        elements.forEach {
             if (!shape.contentEquals(it.shape)) {
                 throw ShapeMismatchException(shape, it.shape)
             }
@@ -43,7 +43,7 @@ abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field
     /**
      * Element-by-element addition
      */
-    override fun add(a: NDArray<T>, b: NDArray<T>): NDArray<T> {
+    override fun add(a: NDElement<T>, b: NDElement<T>): NDElement<T> {
         checkShape(a, b)
         return produce { with(field) { a[it] + b[it] } }
     }
@@ -51,18 +51,18 @@ abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field
     /**
      * Multiply all elements by cinstant
      */
-    override fun multiply(a: NDArray<T>, k: Double): NDArray<T> {
+    override fun multiply(a: NDElement<T>, k: Double): NDElement<T> {
         checkShape(a)
         return produce { with(field) { a[it] * k } }
     }
 
-    override val one: NDArray<T>
+    override val one: NDElement<T>
         get() = produce { this.field.one }
 
     /**
      * Element-by-element multiplication
      */
-    override fun multiply(a: NDArray<T>, b: NDArray<T>): NDArray<T> {
+    override fun multiply(a: NDElement<T>, b: NDElement<T>): NDElement<T> {
         checkShape(a)
         return produce { with(field) { a[it] * b[it] } }
     }
@@ -70,7 +70,7 @@ abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field
     /**
      * Element-by-element division
      */
-    override fun divide(a: NDArray<T>, b: NDArray<T>): NDArray<T> {
+    override fun divide(a: NDElement<T>, b: NDElement<T>): NDElement<T> {
         checkShape(a)
         return produce { with(field) { a[it] / b[it] } }
     }
@@ -78,12 +78,12 @@ abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field
     /**
      * Reverse sum operation
      */
-    operator fun <T> T.plus(arg: NDArray<T>): NDArray<T> = arg + this
+    operator fun <T> T.plus(arg: NDElement<T>): NDElement<T> = arg + this
 
     /**
      * Reverse minus operation
      */
-    operator fun <T> T.minus(arg: NDArray<T>): NDArray<T> = arg.transform { _, value ->
+    operator fun <T> T.minus(arg: NDElement<T>): NDElement<T> = arg.transform { _, value ->
         with(arg.context.field) {
             this@minus - value
         }
@@ -92,12 +92,12 @@ abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field
     /**
      * Reverse product operation
      */
-    operator fun <T> T.times(arg: NDArray<T>): NDArray<T> = arg * this
+    operator fun <T> T.times(arg: NDElement<T>): NDElement<T> = arg * this
 
     /**
      * Reverse division operation
      */
-    operator fun <T> T.div(arg: NDArray<T>): NDArray<T> = arg.transform { _, value ->
+    operator fun <T> T.div(arg: NDElement<T>): NDElement<T> = arg.transform { _, value ->
         with(arg.context.field) {
             this@div / value
         }
@@ -107,37 +107,37 @@ abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field
 /**
  *  Immutable [NDStructure] coupled to the context. Emulates Python ndarray
  */
-class NDArray<T>(override val context: NDField<T>, private val structure: NDStructure<T>) : FieldElement<NDArray<T>, NDField<T>>, NDStructure<T> by structure {
+class NDElement<T>(override val context: NDField<T>, private val structure: NDStructure<T>) : FieldElement<NDElement<T>, NDField<T>>, NDStructure<T> by structure {
 
     //TODO ensure structure is immutable
 
-    override val self: NDArray<T>
+    override val self: NDElement<T>
         get() = this
 
-    inline fun transform(crossinline action: (IntArray, T) -> T): NDArray<T> = context.produce { action(it, get(*it)) }
-    inline fun transform(crossinline action: (T) -> T): NDArray<T> = context.produce { action(get(*it)) }
+    inline fun transform(crossinline action: (IntArray, T) -> T): NDElement<T> = context.produce { action(it, get(*it)) }
+    inline fun transform(crossinline action: (T) -> T): NDElement<T> = context.produce { action(get(*it)) }
 }
 
 /**
  * Element by element application of any operation on elements to the whole array. Just like in numpy
  */
-operator fun <T> Function1<T, T>.invoke(ndArray: NDArray<T>): NDArray<T> = ndArray.transform { _, value -> this(value) }
+operator fun <T> Function1<T, T>.invoke(ndElement: NDElement<T>): NDElement<T> = ndElement.transform { _, value -> this(value) }
 
 /* plus and minus */
 
 /**
- * Summation operation for [NDArray] and single element
+ * Summation operation for [NDElement] and single element
  */
-operator fun <T> NDArray<T>.plus(arg: T): NDArray<T> = transform { _, value ->
+operator fun <T> NDElement<T>.plus(arg: T): NDElement<T> = transform { _, value ->
     with(context.field) {
         arg + value
     }
 }
 
 /**
- * Subtraction operation between [NDArray] and single element
+ * Subtraction operation between [NDElement] and single element
  */
-operator fun <T> NDArray<T>.minus(arg: T): NDArray<T> = transform { _, value ->
+operator fun <T> NDElement<T>.minus(arg: T): NDElement<T> = transform { _, value ->
     with(context.field) {
         arg - value
     }
@@ -146,18 +146,18 @@ operator fun <T> NDArray<T>.minus(arg: T): NDArray<T> = transform { _, value ->
 /* prod and div */
 
 /**
- * Product operation for [NDArray] and single element
+ * Product operation for [NDElement] and single element
  */
-operator fun <T> NDArray<T>.times(arg: T): NDArray<T> = transform { _, value ->
+operator fun <T> NDElement<T>.times(arg: T): NDElement<T> = transform { _, value ->
     with(context.field) {
         arg * value
     }
 }
 
 /**
- * Division operation between [NDArray] and single element
+ * Division operation between [NDElement] and single element
  */
-operator fun <T> NDArray<T>.div(arg: T): NDArray<T> = transform { _, value ->
+operator fun <T> NDElement<T>.div(arg: T): NDElement<T> = transform { _, value ->
     with(context.field) {
         arg / value
     }
@@ -173,23 +173,23 @@ object NDArrays {
     /**
      * Create a platform-optimized NDArray of doubles
      */
-    fun realNDArray(shape: IntArray, initializer: (IntArray) -> Double = { 0.0 }): NDArray<Double> {
+    fun realNDArray(shape: IntArray, initializer: (IntArray) -> Double = { 0.0 }): NDElement<Double> {
         return ExtendedNDField(shape, DoubleField).produce(initializer)
     }
 
-    fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDArray<Double> {
+    fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDElement<Double> {
         return realNDArray(intArrayOf(dim)) { initializer(it[0]) }
     }
 
-    fun real2DArray(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }): NDArray<Double> {
+    fun real2DArray(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }): NDElement<Double> {
         return realNDArray(intArrayOf(dim1, dim2)) { initializer(it[0], it[1]) }
     }
 
-    fun real3DArray(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }): NDArray<Double> {
+    fun real3DArray(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }): NDElement<Double> {
         return realNDArray(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
     }
 
-    inline fun produceReal(shape: IntArray, block: ExtendedNDField<Double>.() -> NDArray<Double>) = ExtendedNDField(shape, DoubleField).run(block)
+    inline fun produceReal(shape: IntArray, block: ExtendedNDField<Double>.() -> NDElement<Double>) = ExtendedNDField(shape, DoubleField).run(block)
 
 //    /**
 //     * Simple boxing NDField
@@ -199,7 +199,7 @@ object NDArrays {
     /**
      * Simple boxing NDArray
      */
-    fun <T : Any> create(field: Field<T>, shape: IntArray, initializer: (IntArray) -> T): NDArray<T> {
+    fun <T : Any> create(field: Field<T>, shape: IntArray, initializer: (IntArray) -> T): NDElement<T> {
         return GenericNDField(shape, field).produce { initializer(it) }
     }
 }
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
index 844535fa8..cd58a63b4 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
@@ -1,7 +1,5 @@
 package scientifik.kmath.structures
 
-import scientifik.kmath.linear.Vector
-import scientifik.kmath.linear.VectorL2Norm
 import scientifik.kmath.operations.Norm
 import scientifik.kmath.structures.NDArrays.produceReal
 import scientifik.kmath.structures.NDArrays.real2DArray
@@ -53,8 +51,8 @@ class NumberNDFieldTest {
         assertEquals(2.0, result[0,2])
     }
 
-    object L2Norm: Norm<NDArray<out Number>, Double> {
-        override fun norm(arg: NDArray<out Number>): Double {
+    object L2Norm: Norm<NDElement<out Number>, Double> {
+        override fun norm(arg: NDElement<out Number>): Double {
             return kotlin.math.sqrt(arg.sumByDouble { it.second.toDouble() })
         }
     }

From 8f5276bdcbe51658729f7828fa8ca7ba0b08f6c0 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Fri, 23 Nov 2018 14:08:28 +0300
Subject: [PATCH 07/13] Broaden scope of VectorSpace to include Points

---
 .../scientifik/kmath/histogram/Histogram.kt   |  13 +-
 .../kmath/histogram/PhantomHistogram.kt       |  41 ++++--
 .../scientifik/kmath/linear/LinearAlgrebra.kt | 100 +--------------
 .../kotlin/scientifik/kmath/linear/Vector.kt  | 118 ++++++++++++++++++
 .../scientifik/kmath/structures/Buffers.kt    |   2 +-
 5 files changed, 151 insertions(+), 123 deletions(-)
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
index e143bd90d..f605f6df3 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
@@ -1,6 +1,5 @@
 package scientifik.kmath.histogram
 
-import scientifik.kmath.operations.Space
 import scientifik.kmath.structures.ArrayBuffer
 import scientifik.kmath.structures.Buffer
 
@@ -55,14 +54,4 @@ fun <T: Any> MutableHistogram<T,*>.fill(sequence: Iterable<Point<T>>) = sequence
 /**
  * Pass a sequence builder into histogram
  */
-fun <T: Any> MutableHistogram<T, *>.fill(buider: suspend SequenceScope<Point<T>>.() -> Unit) = fill(sequence(buider).asIterable())
-
-/**
- * A space to perform arithmetic operations on histograms
- */
-interface HistogramSpace<T: Any, B : Bin<T>, H : Histogram<T,B>> : Space<H> {
-    /**
-     * Rules for performing operations on bins
-     */
-    val binSpace: Space<Bin<T>>
-}
\ No newline at end of file
+fun <T: Any> MutableHistogram<T, *>.fill(buider: suspend SequenceScope<Point<T>>.() -> Unit) = fill(sequence(buider).asIterable())
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
index 0bd5de0e2..041eccfe2 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
@@ -1,23 +1,38 @@
 package scientifik.kmath.histogram
 
-import scientifik.kmath.linear.RealVector
+import scientifik.kmath.linear.Vector
+import scientifik.kmath.operations.Space
 import scientifik.kmath.structures.NDStructure
 
-class BinTemplate(val center: RealVector, val sizes: RealVector) {
-    fun contains(vector: Point<out Double>): Boolean {
+class BinTemplate<T : Comparable<T>>(val center: Vector<T>, val sizes: Vector<T>) {
+    fun contains(vector: Point<out T>): Boolean {
         if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}")
-        return vector.asSequence().mapIndexed { i, value -> value in (center[i] - sizes[i] / 2)..(center[i] + sizes[i] / 2) }.all { it }
+        val upper = center + sizes/2.0
+        val lower = center - sizes/2.0
+        return vector.asSequence().mapIndexed { i, value ->
+            value in lower[i]..upper[i]
+        }.all { it }
     }
 }
 
-class PhantomBin(val template: BinTemplate, override val value: Number) : Bin<Double> {
+/**
+ * A space to perform arithmetic operations on histograms
+ */
+interface HistogramSpace<T : Any, B : Bin<T>, H : Histogram<T, B>> : Space<H> {
+    /**
+     * Rules for performing operations on bins
+     */
+    val binSpace: Space<Bin<T>>
+}
 
-    override fun contains(vector: Point<out Double>): Boolean = template.contains(vector)
+class PhantomBin<T : Comparable<T>>(val template: BinTemplate<T>, override val value: Number) : Bin<T> {
+
+    override fun contains(vector: Point<out T>): Boolean = template.contains(vector)
 
     override val dimension: Int
         get() = template.center.size
 
-    override val center: Point<Double>
+    override val center: Point<T>
         get() = template.center
 
 }
@@ -26,19 +41,19 @@ class PhantomBin(val template: BinTemplate, override val value: Number) : Bin<Do
  * Immutable histogram with explicit structure for content and additional external bin description.
  * Bin search is slow, but full histogram algebra is supported.
  */
-class PhantomHistogram(
-        val bins: Map<BinTemplate, IntArray>,
+class PhantomHistogram<T : Comparable<T>>(
+        val bins: Map<BinTemplate<T>, IntArray>,
         val data: NDStructure<Double>
-) : Histogram<Double, PhantomBin> {
+) : Histogram<T, PhantomBin<T>> {
 
     override val dimension: Int
         get() = data.dimension
 
-    override fun iterator(): Iterator<PhantomBin> {
-        return bins.asSequence().map {entry-> PhantomBin(entry.key,data[entry.value]) }.iterator()
+    override fun iterator(): Iterator<PhantomBin<T>> {
+        return bins.asSequence().map { entry -> PhantomBin(entry.key, data[entry.value]) }.iterator()
     }
 
-    override fun get(point: Point<out Double>): PhantomBin? {
+    override fun get(point: Point<out T>): PhantomBin<T>? {
         val template = bins.keys.find { it.contains(point) }
         return template?.let { PhantomBin(it, data[bins[it]!!]) }
     }
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
index 62852ff5d..9e8692b42 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
@@ -1,7 +1,9 @@
 package scientifik.kmath.linear
 
 import scientifik.kmath.operations.*
-import scientifik.kmath.structures.*
+import scientifik.kmath.structures.ExtendedNDField
+import scientifik.kmath.structures.GenericNDField
+import scientifik.kmath.structures.NDField
 
 /**
  * The space for linear elements. Supports scalar product alongside with standard linear operations.
@@ -139,49 +141,7 @@ interface Matrix<T : Any> : SpaceElement<Matrix<T>, MatrixSpace<T>> {
 }
 
 
-/**
- * A linear space for vectors
- */
-abstract class VectorSpace<T : Any>(val size: Int, val field: Field<T>) : Space<Vector<T>> {
 
-    abstract fun produce(initializer: (Int) -> T): Vector<T>
-
-    override val zero: Vector<T> by lazy { produce { field.zero } }
-
-    override fun add(a: Vector<T>, b: Vector<T>): Vector<T> = produce { with(field) { a[it] + b[it] } }
-
-    override fun multiply(a: Vector<T>, k: Double): Vector<T> = produce { with(field) { a[it] * k } }
-}
-
-
-interface Vector<T : Any> : SpaceElement<Vector<T>, VectorSpace<T>>, Buffer<T>, Iterable<T> {
-    override val size: Int get() = context.size
-
-    companion object {
-        /**
-         * Create vector with custom field
-         */
-        fun <T : Any> of(size: Int, field: Field<T>, initializer: (Int) -> T) =
-                ArrayVector(ArrayVectorSpace(size, field), initializer)
-
-        /**
-         * Create vector of [Double]
-         */
-        fun ofReal(size: Int, initializer: (Int) -> Double) =
-                ArrayVector(ArrayVectorSpace(size, DoubleField, realNDFieldFactory), initializer)
-
-        fun ofReal(vararg point: Double) = point.toVector()
-
-        fun equals(v1: Vector<*>, v2: Vector<*>): Boolean {
-            if (v1 === v2) return true
-            if (v1.context != v2.context) return false
-            for (i in 0 until v2.size) {
-                if (v1[i] != v2[i]) return false
-            }
-            return true
-        }
-    }
-}
 
 typealias NDFieldFactory<T> = (IntArray) -> NDField<T>
 
@@ -210,61 +170,7 @@ class ArrayMatrixSpace<T : Any>(
     }
 }
 
-class ArrayVectorSpace<T : Any>(
-        size: Int,
-        field: Field<T>,
-        val ndFactory: NDFieldFactory<T> = genericNDFieldFactory(field)
-) : VectorSpace<T>(size, field) {
-    val ndField by lazy {
-        ndFactory(intArrayOf(size))
-    }
 
-    override fun produce(initializer: (Int) -> T): Vector<T> = ArrayVector(this, initializer)
-}
-
-/**
- * Member of [ArrayMatrixSpace] which wraps 2-D array
- */
-class ArrayMatrix<T : Any> internal constructor(override val context: ArrayMatrixSpace<T>, val element: NDElement<T>) : Matrix<T> {
-
-    constructor(context: ArrayMatrixSpace<T>, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) })
-
-    override val rows: Int get() = context.rows
-
-    override val columns: Int get() = context.columns
-
-    override fun get(i: Int, j: Int): T {
-        return element[i, j]
-    }
-
-    override val self: ArrayMatrix<T> get() = this
-}
-
-
-class ArrayVector<T : Any> internal constructor(override val context: ArrayVectorSpace<T>, val element: NDElement<T>) : Vector<T> {
-
-    constructor(context: ArrayVectorSpace<T>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
-
-    init {
-        if (context.size != element.shape[0]) {
-            error("Array dimension mismatch")
-        }
-    }
-
-    override fun get(index: Int): T {
-        return element[index]
-    }
-
-    override val self: ArrayVector<T> get() = this
-
-    override fun iterator(): Iterator<T> = (0 until size).map { element[it] }.iterator()
-
-    override fun copy(): ArrayVector<T> = ArrayVector(context, element)
-
-    override fun toString(): String = this.joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() }
-}
-
-typealias RealVector = Vector<Double>
 
 /**
  * A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
new file mode 100644
index 000000000..70480542f
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
@@ -0,0 +1,118 @@
+package scientifik.kmath.linear
+
+import scientifik.kmath.histogram.Point
+import scientifik.kmath.operations.DoubleField
+import scientifik.kmath.operations.Field
+import scientifik.kmath.operations.Space
+import scientifik.kmath.operations.SpaceElement
+import scientifik.kmath.structures.NDElement
+import scientifik.kmath.structures.get
+
+/**
+ * A linear space for vectors.
+ * Could be used on any point-like structure
+ */
+abstract class VectorSpace<T : Any>(val size: Int, val field: Field<T>) : Space<Point<T>> {
+
+    abstract fun produce(initializer: (Int) -> T): Vector<T>
+
+    override val zero: Vector<T> by lazy { produce { field.zero } }
+
+    override fun add(a: Point<T>, b: Point<T>): Vector<T> = produce { with(field) { a[it] + b[it] } }
+
+    override fun multiply(a: Point<T>, k: Double): Vector<T> = produce { with(field) { a[it] * k } }
+}
+
+
+/**
+ * A point coupled to the linear space
+ */
+interface Vector<T : Any> : SpaceElement<Point<T>, VectorSpace<T>>, Point<T>, Iterable<T> {
+    override val size: Int get() = context.size
+
+    override operator fun plus(b: Point<T>): Vector<T> = context.add(self, b)
+    override operator fun minus(b: Point<T>): Vector<T> = context.add(self, context.multiply(b, -1.0))
+    override operator fun times(k: Number): Vector<T> = context.multiply(self, k.toDouble())
+    override operator fun div(k: Number): Vector<T> = context.multiply(self, 1.0 / k.toDouble())
+
+    companion object {
+        /**
+         * Create vector with custom field
+         */
+        fun <T : Any> of(size: Int, field: Field<T>, initializer: (Int) -> T) =
+                ArrayVector(ArrayVectorSpace(size, field), initializer)
+
+        /**
+         * Create vector of [Double]
+         */
+        fun ofReal(size: Int, initializer: (Int) -> Double) =
+                ArrayVector(ArrayVectorSpace(size, DoubleField, realNDFieldFactory), initializer)
+
+        fun ofReal(vararg point: Double) = point.toVector()
+
+        fun equals(v1: Vector<*>, v2: Vector<*>): Boolean {
+            if (v1 === v2) return true
+            if (v1.context != v2.context) return false
+            for (i in 0 until v2.size) {
+                if (v1[i] != v2[i]) return false
+            }
+            return true
+        }
+    }
+}
+
+class ArrayVectorSpace<T : Any>(
+        size: Int,
+        field: Field<T>,
+        val ndFactory: NDFieldFactory<T> = genericNDFieldFactory(field)
+) : VectorSpace<T>(size, field) {
+    val ndField by lazy {
+        ndFactory(intArrayOf(size))
+    }
+
+    override fun produce(initializer: (Int) -> T): Vector<T> = ArrayVector(this, initializer)
+}
+
+/**
+ * Member of [ArrayMatrixSpace] which wraps 2-D array
+ */
+class ArrayMatrix<T : Any> internal constructor(override val context: ArrayMatrixSpace<T>, val element: NDElement<T>) : Matrix<T> {
+
+    constructor(context: ArrayMatrixSpace<T>, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) })
+
+    override val rows: Int get() = context.rows
+
+    override val columns: Int get() = context.columns
+
+    override fun get(i: Int, j: Int): T {
+        return element[i, j]
+    }
+
+    override val self: ArrayMatrix<T> get() = this
+}
+
+
+class ArrayVector<T : Any> internal constructor(override val context: VectorSpace<T>, val element: NDElement<T>) : Vector<T> {
+
+    constructor(context: ArrayVectorSpace<T>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
+
+    init {
+        if (context.size != element.shape[0]) {
+            error("Array dimension mismatch")
+        }
+    }
+
+    override fun get(index: Int): T {
+        return element[index]
+    }
+
+    override val self: ArrayVector<T> get() = this
+
+    override fun iterator(): Iterator<T> = (0 until size).map { element[it] }.iterator()
+
+    override fun copy(): ArrayVector<T> = ArrayVector(context, element)
+
+    override fun toString(): String = this.joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() }
+}
+
+typealias RealVector = Vector<Double>
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
index b68072332..052f9caa4 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
@@ -2,7 +2,7 @@ package scientifik.kmath.structures
 
 
 /**
- * A generic linear buffer for both primitives and objects
+ * A generic random access structure for both primitives and objects
  */
 interface Buffer<T> : Iterable<T> {
 

From 29da5bad2842a08480e5e1403dd97e8d7055a2d8 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Fri, 23 Nov 2018 15:46:14 +0300
Subject: [PATCH 08/13] Added phantom histograms

---
 .../scientifik/kmath/histogram/FastHistogram.kt   | 15 +++++++++++++++
 .../kmath/histogram/PhantomHistogram.kt           |  5 +++--
 2 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
index 621ea90c0..37c4eed2a 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
@@ -84,6 +84,21 @@ class FastHistogram(
 
     override fun iterator(): Iterator<MultivariateBin> = bins.asSequence().map { it.second }.iterator()
 
+    /**
+     * Convert this histogram into NDStructure containing bin values but not bin descriptions
+     */
+    fun asND(): NDStructure<Number> {
+        return ndStructure(this.bins.shape) { bins[it].value }
+    }
+
+    /**
+     * Create a phantom lightweight immutable copy of this histogram
+     */
+    fun asPhantom(): PhantomHistogram<Double> {
+        val binTemplates = bins.associate { (index, bin) -> BinTemplate<Double>(bin.center, bin.sizes) to index }
+        return PhantomHistogram(binTemplates, asND())
+    }
+
     companion object {
 
         /**
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
index 041eccfe2..9ea17c730 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
@@ -4,7 +4,7 @@ import scientifik.kmath.linear.Vector
 import scientifik.kmath.operations.Space
 import scientifik.kmath.structures.NDStructure
 
-class BinTemplate<T : Comparable<T>>(val center: Vector<T>, val sizes: Vector<T>) {
+data class BinTemplate<T : Comparable<T>>(val center: Vector<T>, val sizes: Vector<T>) {
     fun contains(vector: Point<out T>): Boolean {
         if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}")
         val upper = center + sizes/2.0
@@ -40,10 +40,11 @@ class PhantomBin<T : Comparable<T>>(val template: BinTemplate<T>, override val v
 /**
  * Immutable histogram with explicit structure for content and additional external bin description.
  * Bin search is slow, but full histogram algebra is supported.
+ * @param bins map a template into structure index
  */
 class PhantomHistogram<T : Comparable<T>>(
         val bins: Map<BinTemplate<T>, IntArray>,
-        val data: NDStructure<Double>
+        val data: NDStructure<Number>
 ) : Histogram<T, PhantomBin<T>> {
 
     override val dimension: Int

From 5e3e0eb09d5c087f7239c70a16a5633b9e77db61 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Fri, 23 Nov 2018 18:42:40 +0300
Subject: [PATCH 09/13] All fields in Structures, Vectors and Matrices are now
 generic

---
 .../scientifik/kmath/histogram/Histogram.kt   |   3 +
 .../kmath/histogram/PhantomHistogram.kt       |   6 +-
 .../kmath/linear/LUDecomposition.kt           |  18 +-
 .../scientifik/kmath/linear/LinearAlgrebra.kt | 194 ++----------------
 .../kotlin/scientifik/kmath/linear/Matrix.kt  | 171 +++++++++++++++
 .../kotlin/scientifik/kmath/linear/Vector.kt  |  49 +++--
 .../kmath/structures/ExtendedNDField.kt       |  32 +--
 .../scientifik/kmath/structures/NDField.kt    | 113 +++++-----
 .../kmath/structures/NumberNDFieldTest.kt     |  14 +-
 9 files changed, 306 insertions(+), 294 deletions(-)
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
index f605f6df3..8e429087a 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
@@ -2,6 +2,7 @@ package scientifik.kmath.histogram
 
 import scientifik.kmath.structures.ArrayBuffer
 import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.DoubleBuffer
 
 typealias Point<T> = Buffer<T>
 
@@ -49,6 +50,8 @@ interface MutableHistogram<T: Any, out B : Bin<T>>: Histogram<T,B>{
 
 fun <T: Any> MutableHistogram<T,*>.put(vararg point: T) = put(ArrayBuffer(point))
 
+fun MutableHistogram<Double,*>.put(vararg point: Number) = put(DoubleBuffer(point.map { it.toDouble() }.toDoubleArray()))
+
 fun <T: Any> MutableHistogram<T,*>.fill(sequence: Iterable<Point<T>>) = sequence.forEach { put(it) }
 
 /**
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
index 9ea17c730..25a1bf33a 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
@@ -4,11 +4,11 @@ import scientifik.kmath.linear.Vector
 import scientifik.kmath.operations.Space
 import scientifik.kmath.structures.NDStructure
 
-data class BinTemplate<T : Comparable<T>>(val center: Vector<T>, val sizes: Vector<T>) {
+data class BinTemplate<T : Comparable<T>>(val center: Vector<T, *>, val sizes: Vector<T, *>) {
     fun contains(vector: Point<out T>): Boolean {
         if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}")
-        val upper = center + sizes/2.0
-        val lower = center - sizes/2.0
+        val upper = center + sizes / 2.0
+        val lower = center - sizes / 2.0
         return vector.asSequence().mapIndexed { i, value ->
             value in lower[i]..upper[i]
         }.all { it }
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
index 0b71405c6..223e63419 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
@@ -1,5 +1,7 @@
 package scientifik.kmath.linear
 
+import scientifik.kmath.operations.DoubleField
+import scientifik.kmath.operations.Field
 import scientifik.kmath.structures.MutableNDStructure
 import scientifik.kmath.structures.NDStructure
 import scientifik.kmath.structures.genericNdStructure
@@ -9,7 +11,7 @@ import kotlin.math.absoluteValue
 /**
  * Implementation based on Apache common-maths LU-decomposition
  */
-abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
+abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matrix<T, F>) {
 
     private val field get() = matrix.context.field
     /** Entries of LU decomposition.  */
@@ -31,7 +33,7 @@ abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
      * L is a lower-triangular matrix
      * @return the L matrix (or null if decomposed matrix is singular)
      */
-    val l: Matrix<out T> by lazy {
+    val l: Matrix<out T, F> by lazy {
         matrix.context.produce { i, j ->
             when {
                 j < i -> lu[i, j]
@@ -48,7 +50,7 @@ abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
      * U is an upper-triangular matrix
      * @return the U matrix (or null if decomposed matrix is singular)
      */
-    val u: Matrix<out T> by lazy {
+    val u: Matrix<out T, F> by lazy {
         matrix.context.produce { i, j ->
             if (j >= i) lu[i, j] else field.zero
         }
@@ -64,7 +66,7 @@ abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
      * @return the P rows permutation matrix (or null if decomposed matrix is singular)
      * @see .getPivot
      */
-    val p: Matrix<out T> by lazy {
+    val p: Matrix<out T, F> by lazy {
         matrix.context.produce { i, j ->
             //TODO ineffective. Need sparse matrix for that
             if (j == pivot[i]) field.one else field.zero
@@ -181,7 +183,7 @@ abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
 
 }
 
-class RealLUDecomposition(matrix: Matrix<Double>, private val singularityThreshold: Double = DEFAULT_TOO_SMALL) : LUDecomposition<Double>(matrix) {
+class RealLUDecomposition(matrix: RealMatrix, private val singularityThreshold: Double = DEFAULT_TOO_SMALL) : LUDecomposition<Double, DoubleField>(matrix) {
     override fun isSingular(value: Double): Boolean {
         return value.absoluteValue < singularityThreshold
     }
@@ -194,12 +196,12 @@ class RealLUDecomposition(matrix: Matrix<Double>, private val singularityThresho
 
 
 /** Specialized solver.  */
-object RealLUSolver : LinearSolver<Double> {
+object RealLUSolver : LinearSolver<Double, DoubleField> {
 
 
-    fun decompose(mat: Matrix<Double>, threshold: Double = 1e-11): RealLUDecomposition = RealLUDecomposition(mat, threshold)
+    fun decompose(mat: Matrix<Double, DoubleField>, threshold: Double = 1e-11): RealLUDecomposition = RealLUDecomposition(mat, threshold)
 
-    override fun solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> {
+    override fun solve(a: RealMatrix, b: RealMatrix): RealMatrix {
         val decomposition = decompose(a, a.context.field.zero)
 
         if (b.rows != a.rows) {
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
index 9e8692b42..acd27dbd9 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
@@ -1,184 +1,17 @@
 package scientifik.kmath.linear
 
-import scientifik.kmath.operations.*
-import scientifik.kmath.structures.ExtendedNDField
-import scientifik.kmath.structures.GenericNDField
-import scientifik.kmath.structures.NDField
-
-/**
- * The space for linear elements. Supports scalar product alongside with standard linear operations.
- * @param T type of individual element of the vector or matrix
- * @param V the type of vector space element
- */
-abstract class MatrixSpace<T : Any>(val rows: Int, val columns: Int, val field: Field<T>) : Space<Matrix<T>> {
-
-    /**
-     * Produce the element of this space
-     */
-    abstract fun produce(initializer: (Int, Int) -> T): Matrix<T>
-
-    /**
-     * Produce new matrix space with given dimensions. The space produced could be raised from cache since [MatrixSpace] does not have mutable elements
-     */
-    abstract fun produceSpace(rows: Int, columns: Int): MatrixSpace<T>
-
-    override val zero: Matrix<T> by lazy {
-        produce { _, _ -> field.zero }
-    }
-
-//    val one: Matrix<T> by lazy {
-//        produce { i, j -> if (i == j) field.one else field.zero }
-//    }
-
-    override fun add(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
-        return produce { i, j -> with(field) { a[i, j] + b[i, j] } }
-    }
-
-    override fun multiply(a: Matrix<T>, k: Double): Matrix<T> {
-        //TODO it is possible to implement scalable linear elements which normed values and adjustable scale to save memory and processing poser
-        return produce { i, j -> with(field) { a[i, j] * k } }
-    }
-
-    /**
-     * Dot product. Throws exception on dimension mismatch
-     */
-    fun multiply(a: Matrix<T>, b: Matrix<T>): Matrix<T> {
-        if (a.rows != b.columns) {
-            //TODO replace by specific exception
-            error("Dimension mismatch in linear structure dot product: [${a.rows},${a.columns}]*[${b.rows},${b.columns}]")
-        }
-        return produceSpace(a.rows, b.columns).produce { i, j ->
-            (0 until a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) }
-        }
-    }
-
-    override fun equals(other: Any?): Boolean {
-        if (this === other) return true
-        if (other !is MatrixSpace<*>) return false
-
-        if (rows != other.rows) return false
-        if (columns != other.columns) return false
-        if (field != other.field) return false
-
-        return true
-    }
-
-    override fun hashCode(): Int {
-        var result = rows
-        result = 31 * result + columns
-        result = 31 * result + field.hashCode()
-        return result
-    }
-}
-
-infix fun <T : Any> Matrix<T>.dot(b: Matrix<T>): Matrix<T> = this.context.multiply(this, b)
-
-/**
- * A matrix-like structure
- */
-interface Matrix<T : Any> : SpaceElement<Matrix<T>, MatrixSpace<T>> {
-    /**
-     * Number of rows
-     */
-    val rows: Int
-    /**
-     * Number of columns
-     */
-    val columns: Int
-
-    /**
-     * Get element in row [i] and column [j]. Throws error in case of call ounside structure dimensions
-     */
-    operator fun get(i: Int, j: Int): T
-
-    override val self: Matrix<T>
-        get() = this
-
-    fun transpose(): Matrix<T> {
-        return object : Matrix<T> {
-            override val context: MatrixSpace<T> = this@Matrix.context
-            override val rows: Int = this@Matrix.columns
-            override val columns: Int = this@Matrix.rows
-            override fun get(i: Int, j: Int): T = this@Matrix[j, i]
-        }
-    }
-
-    companion object {
-
-        /**
-         * Create [ArrayMatrix] with custom field
-         */
-        fun <T : Any> of(rows: Int, columns: Int, field: Field<T>, initializer: (Int, Int) -> T) =
-                ArrayMatrix(ArrayMatrixSpace(rows, columns, field), initializer)
-
-        /**
-         * Create [ArrayMatrix] of doubles. The implementation in general should be faster than generic one due to boxing.
-         */
-        fun ofReal(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
-                ArrayMatrix(ArrayMatrixSpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
-
-        /**
-         * Create a diagonal value matrix. By default value equals [Field.one].
-         */
-        fun <T : Any> diagonal(rows: Int, columns: Int, field: Field<T>, values: (Int) -> T = { field.one }): Matrix<T> {
-            return of(rows, columns, field) { i, j -> if (i == j) values(i) else field.zero }
-        }
-
-        /**
-         * Equality check on two generic matrices
-         */
-        fun equals(mat1: Matrix<*>, mat2: Matrix<*>): Boolean {
-            if (mat1 === mat2) return true
-            if (mat1.context != mat2.context) return false
-            for (i in 0 until mat1.rows) {
-                for (j in 0 until mat2.columns) {
-                    if (mat1[i, j] != mat2[i, j]) return false
-                }
-            }
-            return true
-        }
-    }
-}
-
-
-
-
-typealias NDFieldFactory<T> = (IntArray) -> NDField<T>
-
-internal fun <T : Any> genericNDFieldFactory(field: Field<T>): NDFieldFactory<T> = { index -> GenericNDField(index, field) }
-internal val realNDFieldFactory: NDFieldFactory<Double> = { index -> ExtendedNDField(index, DoubleField) }
-
-
-/**
- * NDArray-based implementation of vector space. By default uses slow [GenericNDField], but could be overridden with custom [NDField] factory.
- */
-class ArrayMatrixSpace<T : Any>(
-        rows: Int,
-        columns: Int,
-        field: Field<T>,
-        val ndFactory: NDFieldFactory<T> = genericNDFieldFactory(field)
-) : MatrixSpace<T>(rows, columns, field) {
-
-    val ndField by lazy {
-        ndFactory(intArrayOf(rows, columns))
-    }
-
-    override fun produce(initializer: (Int, Int) -> T): Matrix<T> = ArrayMatrix(this, initializer)
-
-    override fun produceSpace(rows: Int, columns: Int): ArrayMatrixSpace<T> {
-        return ArrayMatrixSpace(rows, columns, field, ndFactory)
-    }
-}
-
+import scientifik.kmath.operations.DoubleField
+import scientifik.kmath.operations.Field
+import scientifik.kmath.operations.Norm
 
 
 /**
  * A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors
  */
-interface LinearSolver<T : Any> {
-    fun solve(a: Matrix<T>, b: Matrix<T>): Matrix<T>
-    fun solve(a: Matrix<T>, b: Vector<T>): Vector<T> = solve(a, b.toMatrix()).toVector()
-    fun inverse(a: Matrix<T>): Matrix<T> = solve(a, Matrix.diagonal(a.rows, a.columns, a.context.field))
+interface LinearSolver<T : Any, F : Field<T>> {
+    fun solve(a: Matrix<T, F>, b: Matrix<T, F>): Matrix<T, F>
+    fun solve(a: Matrix<T, F>, b: Vector<T, F>): Vector<T, F> = solve(a, b.toMatrix()).toVector()
+    fun inverse(a: Matrix<T, F>): Matrix<T, F> = solve(a, Matrix.diagonal(a.rows, a.columns, a.context.field))
 }
 
 /**
@@ -192,7 +25,7 @@ fun List<Double>.toVector() = Vector.ofReal(this.size) { this[it] }
 /**
  * Convert matrix to vector if it is possible
  */
-fun <T : Any> Matrix<T>.toVector(): Vector<T> {
+fun <T : Any, F : Field<T>> Matrix<T, F>.toVector(): Vector<T, F> {
     return when {
         this.columns == 1 -> {
 //            if (this is ArrayMatrix) {
@@ -208,7 +41,7 @@ fun <T : Any> Matrix<T>.toVector(): Vector<T> {
     }
 }
 
-fun <T : Any> Vector<T>.toMatrix(): Matrix<T> {
+fun <T : Any, F : Field<T>> Vector<T, F>.toMatrix(): Matrix<T, F> {
 //    return if (this is ArrayVector) {
 //        //Reuse existing underlying array
 //        ArrayMatrix(ArrayMatrixSpace(size, 1, context.field, context.ndFactory), array)
@@ -216,11 +49,14 @@ fun <T : Any> Vector<T>.toMatrix(): Matrix<T> {
 //        //Generic vector
 //        matrix(size, 1, context.field) { i, j -> get(i) }
 //    }
-    return Matrix.of(size, 1, context.field) { i, _ -> get(i) }
+    return Matrix.of(size, 1, context.space) { i, _ -> get(i) }
 }
 
-object VectorL2Norm: Norm<Vector<out Number>, Double> {
-    override fun norm(arg: Vector<out Number>): Double {
+object VectorL2Norm : Norm<Vector<out Number, *>, Double> {
+    override fun norm(arg: Vector<out Number, *>): Double {
         return kotlin.math.sqrt(arg.sumByDouble { it.toDouble() })
     }
 }
+
+typealias RealVector = Vector<Double, DoubleField>
+typealias RealMatrix = Matrix<Double, DoubleField>
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
new file mode 100644
index 000000000..f7ae04849
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
@@ -0,0 +1,171 @@
+package scientifik.kmath.linear
+
+import scientifik.kmath.operations.*
+import scientifik.kmath.structures.ExtendedNDField
+import scientifik.kmath.structures.GenericNDField
+import scientifik.kmath.structures.NDField
+
+/**
+ * The space for linear elements. Supports scalar product alongside with standard linear operations.
+ * @param T type of individual element of the vector or matrix
+ * @param V the type of vector space element
+ */
+abstract class MatrixSpace<T : Any, F : Ring<T>>(val rows: Int, val columns: Int, val field: F) : Space<Matrix<T, F>> {
+
+    /**
+     * Produce the element of this space
+     */
+    abstract fun produce(initializer: (Int, Int) -> T): Matrix<T, F>
+
+    /**
+     * Produce new matrix space with given dimensions. The space produced could be raised from cache since [MatrixSpace] does not have mutable elements
+     */
+    abstract fun produceSpace(rows: Int, columns: Int): MatrixSpace<T, F>
+
+    override val zero: Matrix<T, F> by lazy {
+        produce { _, _ -> field.zero }
+    }
+
+//    val one: Matrix<T> by lazy {
+//        produce { i, j -> if (i == j) field.one else field.zero }
+//    }
+
+    override fun add(a: Matrix<T, F>, b: Matrix<T, F>): Matrix<T, F> {
+        return produce { i, j -> with(field) { a[i, j] + b[i, j] } }
+    }
+
+    override fun multiply(a: Matrix<T, F>, k: Double): Matrix<T, F> {
+        //TODO it is possible to implement scalable linear elements which normed values and adjustable scale to save memory and processing poser
+        return produce { i, j -> with(field) { a[i, j] * k } }
+    }
+
+    /**
+     * Dot product. Throws exception on dimension mismatch
+     */
+    fun multiply(a: Matrix<T, F>, b: Matrix<T, F>): Matrix<T, F> {
+        if (a.rows != b.columns) {
+            //TODO replace by specific exception
+            error("Dimension mismatch in linear structure dot product: [${a.rows},${a.columns}]*[${b.rows},${b.columns}]")
+        }
+        return produceSpace(a.rows, b.columns).produce { i, j ->
+            (0 until a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) }
+        }
+    }
+
+    override fun equals(other: Any?): Boolean {
+        if (this === other) return true
+        if (other !is MatrixSpace<*,*>) return false
+
+        if (rows != other.rows) return false
+        if (columns != other.columns) return false
+        if (field != other.field) return false
+
+        return true
+    }
+
+    override fun hashCode(): Int {
+        var result = rows
+        result = 31 * result + columns
+        result = 31 * result + field.hashCode()
+        return result
+    }
+}
+
+infix fun <T : Any, F : Field<T>> Matrix<T, F>.dot(b: Matrix<T, F>): Matrix<T, F> = this.context.multiply(this, b)
+
+/**
+ * A matrix-like structure
+ */
+interface Matrix<T : Any, F: Ring<T>> : SpaceElement<Matrix<T, F>, MatrixSpace<T, F>> {
+    /**
+     * Number of rows
+     */
+    val rows: Int
+    /**
+     * Number of columns
+     */
+    val columns: Int
+
+    /**
+     * Get element in row [i] and column [j]. Throws error in case of call ounside structure dimensions
+     */
+    operator fun get(i: Int, j: Int): T
+
+    override val self: Matrix<T, F>
+        get() = this
+
+    fun transpose(): Matrix<T, F> {
+        return object : Matrix<T, F> {
+            override val context: MatrixSpace<T, F> = this@Matrix.context
+            override val rows: Int = this@Matrix.columns
+            override val columns: Int = this@Matrix.rows
+            override fun get(i: Int, j: Int): T = this@Matrix[j, i]
+        }
+    }
+
+    companion object {
+
+        /**
+         * Create [ArrayMatrix] with custom field
+         */
+        fun <T : Any, F: Field<T>> of(rows: Int, columns: Int, field: F, initializer: (Int, Int) -> T) =
+                ArrayMatrix(ArrayMatrixSpace(rows, columns, field), initializer)
+
+        /**
+         * Create [ArrayMatrix] of doubles. The implementation in general should be faster than generic one due to boxing.
+         */
+        fun ofReal(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
+                ArrayMatrix(ArrayMatrixSpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
+
+        /**
+         * Create a diagonal value matrix. By default value equals [Field.one].
+         */
+        fun <T : Any, F: Field<T>> diagonal(rows: Int, columns: Int, field: F, values: (Int) -> T = { field.one }): Matrix<T, F> {
+            return of(rows, columns, field) { i, j -> if (i == j) values(i) else field.zero }
+        }
+
+        /**
+         * Equality check on two generic matrices
+         */
+        fun equals(mat1: Matrix<*, *>, mat2: Matrix<*, *>): Boolean {
+            if (mat1 === mat2) return true
+            if (mat1.context != mat2.context) return false
+            for (i in 0 until mat1.rows) {
+                for (j in 0 until mat2.columns) {
+                    if (mat1[i, j] != mat2[i, j]) return false
+                }
+            }
+            return true
+        }
+    }
+}
+
+
+
+
+typealias NDFieldFactory<T, F> = (IntArray) -> NDField<T, F>
+
+internal fun <T : Any, F : Field<T>> genericNDFieldFactory(field: F): NDFieldFactory<T, F> = { index -> GenericNDField(index, field) }
+internal val realNDFieldFactory: NDFieldFactory<Double, DoubleField> = { index -> ExtendedNDField(index, DoubleField) }
+
+
+/**
+ * NDArray-based implementation of vector space. By default uses slow [GenericNDField], but could be overridden with custom [NDField] factory.
+ */
+class ArrayMatrixSpace<T : Any, F : Field<T>>(
+        rows: Int,
+        columns: Int,
+        field: F,
+        val ndFactory: NDFieldFactory<T, F> = genericNDFieldFactory(field)
+) : MatrixSpace<T, F>(rows, columns, field) {
+
+    val ndField by lazy {
+        ndFactory(intArrayOf(rows, columns))
+    }
+
+    override fun produce(initializer: (Int, Int) -> T): Matrix<T, F> = ArrayMatrix(this, initializer)
+
+    override fun produceSpace(rows: Int, columns: Int): ArrayMatrixSpace<T, F> {
+        return ArrayMatrixSpace(rows, columns, field, ndFactory)
+    }
+}
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
index 70480542f..3e0f91ebe 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
@@ -12,34 +12,34 @@ import scientifik.kmath.structures.get
  * A linear space for vectors.
  * Could be used on any point-like structure
  */
-abstract class VectorSpace<T : Any>(val size: Int, val field: Field<T>) : Space<Point<T>> {
+abstract class VectorSpace<T : Any, S : Space<T>>(val size: Int, val space: S) : Space<Point<T>> {
 
-    abstract fun produce(initializer: (Int) -> T): Vector<T>
+    abstract fun produce(initializer: (Int) -> T): Vector<T, S>
 
-    override val zero: Vector<T> by lazy { produce { field.zero } }
+    override val zero: Vector<T, S> by lazy { produce { space.zero } }
 
-    override fun add(a: Point<T>, b: Point<T>): Vector<T> = produce { with(field) { a[it] + b[it] } }
+    override fun add(a: Point<T>, b: Point<T>): Vector<T, S> = produce { with(space) { a[it] + b[it] } }
 
-    override fun multiply(a: Point<T>, k: Double): Vector<T> = produce { with(field) { a[it] * k } }
+    override fun multiply(a: Point<T>, k: Double): Vector<T, S> = produce { with(space) { a[it] * k } }
 }
 
 
 /**
  * A point coupled to the linear space
  */
-interface Vector<T : Any> : SpaceElement<Point<T>, VectorSpace<T>>, Point<T>, Iterable<T> {
+interface Vector<T : Any, S : Space<T>> : SpaceElement<Point<T>, VectorSpace<T, S>>, Point<T>, Iterable<T> {
     override val size: Int get() = context.size
 
-    override operator fun plus(b: Point<T>): Vector<T> = context.add(self, b)
-    override operator fun minus(b: Point<T>): Vector<T> = context.add(self, context.multiply(b, -1.0))
-    override operator fun times(k: Number): Vector<T> = context.multiply(self, k.toDouble())
-    override operator fun div(k: Number): Vector<T> = context.multiply(self, 1.0 / k.toDouble())
+    override operator fun plus(b: Point<T>): Vector<T, S> = context.add(self, b)
+    override operator fun minus(b: Point<T>): Vector<T, S> = context.add(self, context.multiply(b, -1.0))
+    override operator fun times(k: Number): Vector<T, S> = context.multiply(self, k.toDouble())
+    override operator fun div(k: Number): Vector<T, S> = context.multiply(self, 1.0 / k.toDouble())
 
     companion object {
         /**
          * Create vector with custom field
          */
-        fun <T : Any> of(size: Int, field: Field<T>, initializer: (Int) -> T) =
+        fun <T : Any, F: Field<T>> of(size: Int, field: F, initializer: (Int) -> T) =
                 ArrayVector(ArrayVectorSpace(size, field), initializer)
 
         /**
@@ -50,7 +50,7 @@ interface Vector<T : Any> : SpaceElement<Point<T>, VectorSpace<T>>, Point<T>, It
 
         fun ofReal(vararg point: Double) = point.toVector()
 
-        fun equals(v1: Vector<*>, v2: Vector<*>): Boolean {
+        fun equals(v1: Vector<*,*>, v2: Vector<*,*>): Boolean {
             if (v1 === v2) return true
             if (v1.context != v2.context) return false
             for (i in 0 until v2.size) {
@@ -61,24 +61,24 @@ interface Vector<T : Any> : SpaceElement<Point<T>, VectorSpace<T>>, Point<T>, It
     }
 }
 
-class ArrayVectorSpace<T : Any>(
+class ArrayVectorSpace<T : Any, F : Field<T>>(
         size: Int,
-        field: Field<T>,
-        val ndFactory: NDFieldFactory<T> = genericNDFieldFactory(field)
-) : VectorSpace<T>(size, field) {
+        field: F,
+        val ndFactory: NDFieldFactory<T, F> = genericNDFieldFactory(field)
+) : VectorSpace<T, F>(size, field) {
     val ndField by lazy {
         ndFactory(intArrayOf(size))
     }
 
-    override fun produce(initializer: (Int) -> T): Vector<T> = ArrayVector(this, initializer)
+    override fun produce(initializer: (Int) -> T): Vector<T, F> = ArrayVector(this, initializer)
 }
 
 /**
  * Member of [ArrayMatrixSpace] which wraps 2-D array
  */
-class ArrayMatrix<T : Any> internal constructor(override val context: ArrayMatrixSpace<T>, val element: NDElement<T>) : Matrix<T> {
+class ArrayMatrix<T : Any, F : Field<T>> internal constructor(override val context: ArrayMatrixSpace<T, F>, val element: NDElement<T, F>) : Matrix<T, F> {
 
-    constructor(context: ArrayMatrixSpace<T>, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) })
+    constructor(context: ArrayMatrixSpace<T, F>, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) })
 
     override val rows: Int get() = context.rows
 
@@ -88,13 +88,13 @@ class ArrayMatrix<T : Any> internal constructor(override val context: ArrayMatri
         return element[i, j]
     }
 
-    override val self: ArrayMatrix<T> get() = this
+    override val self: ArrayMatrix<T, F> get() = this
 }
 
 
-class ArrayVector<T : Any> internal constructor(override val context: VectorSpace<T>, val element: NDElement<T>) : Vector<T> {
+class ArrayVector<T : Any, F : Field<T>> internal constructor(override val context: VectorSpace<T, F>, val element: NDElement<T, F>) : Vector<T, F> {
 
-    constructor(context: ArrayVectorSpace<T>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
+    constructor(context: ArrayVectorSpace<T, F>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
 
     init {
         if (context.size != element.shape[0]) {
@@ -106,13 +106,12 @@ class ArrayVector<T : Any> internal constructor(override val context: VectorSpac
         return element[index]
     }
 
-    override val self: ArrayVector<T> get() = this
+    override val self: ArrayVector<T, F> get() = this
 
     override fun iterator(): Iterator<T> = (0 until size).map { element[it] }.iterator()
 
-    override fun copy(): ArrayVector<T> = ArrayVector(context, element)
+    override fun copy(): ArrayVector<T, F> = ArrayVector(context, element)
 
     override fun toString(): String = this.joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() }
 }
 
-typealias RealVector = Vector<Double>
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt
index 08751d934..93ffad5a3 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/ExtendedNDField.kt
@@ -9,33 +9,33 @@ import scientifik.kmath.operations.TrigonometricOperations
 /**
  * NDField that supports [ExtendedField] operations on its elements
  */
-class ExtendedNDField<N: Any>(shape: IntArray, override val field: ExtendedField<N>) : NDField<N>(shape, field),
-        TrigonometricOperations<NDElement<N>>,
-        PowerOperations<NDElement<N>>,
-        ExponentialOperations<NDElement<N>> {
+class ExtendedNDField<N : Any, F : ExtendedField<N>>(shape: IntArray, field: F) : NDField<N, F>(shape, field),
+        TrigonometricOperations<NDElement<N, F>>,
+        PowerOperations<NDElement<N, F>>,
+        ExponentialOperations<NDElement<N, F>> {
 
-    override fun produceStructure(initializer: (IntArray) -> N): NDStructure<N> {
-        return genericNdStructure(shape, initializer)
+    override fun produceStructure(initializer: F.(IntArray) -> N): NDStructure<N> {
+        return genericNdStructure(shape) { field.initializer(it) }
     }
 
-    override fun power(arg: NDElement<N>, pow: Double): NDElement<N> {
-        return arg.transform { d -> with(field){power(d,pow)} }
+    override fun power(arg: NDElement<N, F>, pow: Double): NDElement<N, F> {
+        return arg.transform { d -> with(field) { power(d, pow) } }
     }
 
-    override fun exp(arg: NDElement<N>): NDElement<N> {
-        return arg.transform { d -> with(field){exp(d)} }
+    override fun exp(arg: NDElement<N, F>): NDElement<N, F> {
+        return arg.transform { d -> with(field) { exp(d) } }
     }
 
-    override fun ln(arg: NDElement<N>): NDElement<N> {
-        return arg.transform { d -> with(field){ln(d)} }
+    override fun ln(arg: NDElement<N, F>): NDElement<N, F> {
+        return arg.transform { d -> with(field) { ln(d) } }
     }
 
-    override fun sin(arg: NDElement<N>): NDElement<N> {
-        return arg.transform { d -> with(field){sin(d)} }
+    override fun sin(arg: NDElement<N, F>): NDElement<N, F> {
+        return arg.transform { d -> with(field) { sin(d) } }
     }
 
-    override fun cos(arg: NDElement<N>): NDElement<N> {
-        return arg.transform { d -> with(field){cos(d)} }
+    override fun cos(arg: NDElement<N, F>): NDElement<N, F> {
+        return arg.transform { d -> with(field) { cos(d) } }
     }
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
index 6d3ce05a5..d8e9cc31b 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
@@ -15,24 +15,24 @@ class ShapeMismatchException(val expected: IntArray, val actual: IntArray) : Run
  * @param field - operations field defined on individual array element
  * @param T the type of the element contained in NDArray
  */
-abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field<NDElement<T>> {
+abstract class NDField<T, F : Field<T>>(val shape: IntArray, val field: F) : Field<NDElement<T, F>> {
 
-    abstract fun produceStructure(initializer: (IntArray) -> T): NDStructure<T>
+    abstract fun produceStructure(initializer: F.(IntArray) -> T): NDStructure<T>
 
     /**
      * Create new instance of NDArray using field shape and given initializer
      * The producer takes list of indices as argument and returns contained value
      */
-    fun produce(initializer: (IntArray) -> T): NDElement<T> = NDElement(this, produceStructure(initializer))
+    fun produce(initializer: F.(IntArray) -> T): NDElement<T, F> = NDElement(this, produceStructure(initializer))
 
-    override val zero: NDElement<T> by lazy {
-        produce { this.field.zero }
+    override val zero: NDElement<T, F> by lazy {
+        produce { zero }
     }
 
     /**
      * Check the shape of given NDArray and throw exception if it does not coincide with shape of the field
      */
-    private fun checkShape(vararg elements: NDElement<T>) {
+    private fun checkShape(vararg elements: NDElement<T, F>) {
         elements.forEach {
             if (!shape.contentEquals(it.shape)) {
                 throw ShapeMismatchException(shape, it.shape)
@@ -43,7 +43,7 @@ abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field
     /**
      * Element-by-element addition
      */
-    override fun add(a: NDElement<T>, b: NDElement<T>): NDElement<T> {
+    override fun add(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
         checkShape(a, b)
         return produce { with(field) { a[it] + b[it] } }
     }
@@ -51,18 +51,18 @@ abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field
     /**
      * Multiply all elements by cinstant
      */
-    override fun multiply(a: NDElement<T>, k: Double): NDElement<T> {
+    override fun multiply(a: NDElement<T, F>, k: Double): NDElement<T, F> {
         checkShape(a)
         return produce { with(field) { a[it] * k } }
     }
 
-    override val one: NDElement<T>
-        get() = produce { this.field.one }
+    override val one: NDElement<T, F>
+        get() = produce { one }
 
     /**
      * Element-by-element multiplication
      */
-    override fun multiply(a: NDElement<T>, b: NDElement<T>): NDElement<T> {
+    override fun multiply(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
         checkShape(a)
         return produce { with(field) { a[it] * b[it] } }
     }
@@ -70,65 +70,65 @@ abstract class NDField<T>(val shape: IntArray, open val field: Field<T>) : Field
     /**
      * Element-by-element division
      */
-    override fun divide(a: NDElement<T>, b: NDElement<T>): NDElement<T> {
+    override fun divide(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
         checkShape(a)
         return produce { with(field) { a[it] / b[it] } }
     }
 
-    /**
-     * Reverse sum operation
-     */
-    operator fun <T> T.plus(arg: NDElement<T>): NDElement<T> = arg + this
-
-    /**
-     * Reverse minus operation
-     */
-    operator fun <T> T.minus(arg: NDElement<T>): NDElement<T> = arg.transform { _, value ->
-        with(arg.context.field) {
-            this@minus - value
-        }
-    }
-
-    /**
-     * Reverse product operation
-     */
-    operator fun <T> T.times(arg: NDElement<T>): NDElement<T> = arg * this
-
-    /**
-     * Reverse division operation
-     */
-    operator fun <T> T.div(arg: NDElement<T>): NDElement<T> = arg.transform { _, value ->
-        with(arg.context.field) {
-            this@div / value
-        }
-    }
+//    /**
+//     * Reverse sum operation
+//     */
+//    operator fun T.plus(arg: NDElement<T, F>): NDElement<T, F> = arg + this
+//
+//    /**
+//     * Reverse minus operation
+//     */
+//    operator fun T.minus(arg: NDElement<T, F>): NDElement<T, F> = arg.transform { _, value ->
+//        with(arg.context.field) {
+//            this@minus - value
+//        }
+//    }
+//
+//    /**
+//     * Reverse product operation
+//     */
+//    operator fun T.times(arg: NDElement<T, F>): NDElement<T, F> = arg * this
+//
+//    /**
+//     * Reverse division operation
+//     */
+//    operator fun T.div(arg: NDElement<T, F>): NDElement<T, F> = arg.transform { _, value ->
+//        with(arg.context.field) {
+//            this@div / value
+//        }
+//    }
 }
 
 /**
  *  Immutable [NDStructure] coupled to the context. Emulates Python ndarray
  */
-class NDElement<T>(override val context: NDField<T>, private val structure: NDStructure<T>) : FieldElement<NDElement<T>, NDField<T>>, NDStructure<T> by structure {
+class NDElement<T, F : Field<T>>(override val context: NDField<T, F>, private val structure: NDStructure<T>) : FieldElement<NDElement<T, F>, NDField<T, F>>, NDStructure<T> by structure {
 
     //TODO ensure structure is immutable
 
-    override val self: NDElement<T>
+    override val self: NDElement<T, F>
         get() = this
 
-    inline fun transform(crossinline action: (IntArray, T) -> T): NDElement<T> = context.produce { action(it, get(*it)) }
-    inline fun transform(crossinline action: (T) -> T): NDElement<T> = context.produce { action(get(*it)) }
+    inline fun transform(crossinline action: (IntArray, T) -> T): NDElement<T, F> = context.produce { action(it, get(*it)) }
+    inline fun transform(crossinline action: (T) -> T): NDElement<T, F> = context.produce { action(get(*it)) }
 }
 
 /**
  * Element by element application of any operation on elements to the whole array. Just like in numpy
  */
-operator fun <T> Function1<T, T>.invoke(ndElement: NDElement<T>): NDElement<T> = ndElement.transform { _, value -> this(value) }
+operator fun <T, F : Field<T>> Function1<T, T>.invoke(ndElement: NDElement<T, F>): NDElement<T, F> = ndElement.transform { _, value -> this(value) }
 
 /* plus and minus */
 
 /**
  * Summation operation for [NDElement] and single element
  */
-operator fun <T> NDElement<T>.plus(arg: T): NDElement<T> = transform { _, value ->
+operator fun <T, F : Field<T>> NDElement<T, F>.plus(arg: T): NDElement<T, F> = transform { _, value ->
     with(context.field) {
         arg + value
     }
@@ -137,7 +137,7 @@ operator fun <T> NDElement<T>.plus(arg: T): NDElement<T> = transform { _, value
 /**
  * Subtraction operation between [NDElement] and single element
  */
-operator fun <T> NDElement<T>.minus(arg: T): NDElement<T> = transform { _, value ->
+operator fun <T, F : Field<T>> NDElement<T, F>.minus(arg: T): NDElement<T, F> = transform { _, value ->
     with(context.field) {
         arg - value
     }
@@ -148,7 +148,7 @@ operator fun <T> NDElement<T>.minus(arg: T): NDElement<T> = transform { _, value
 /**
  * Product operation for [NDElement] and single element
  */
-operator fun <T> NDElement<T>.times(arg: T): NDElement<T> = transform { _, value ->
+operator fun <T, F : Field<T>> NDElement<T, F>.times(arg: T): NDElement<T, F> = transform { _, value ->
     with(context.field) {
         arg * value
     }
@@ -157,14 +157,14 @@ operator fun <T> NDElement<T>.times(arg: T): NDElement<T> = transform { _, value
 /**
  * Division operation between [NDElement] and single element
  */
-operator fun <T> NDElement<T>.div(arg: T): NDElement<T> = transform { _, value ->
+operator fun <T, F : Field<T>> NDElement<T, F>.div(arg: T): NDElement<T, F> = transform { _, value ->
     with(context.field) {
         arg / value
     }
 }
 
-class GenericNDField<T : Any>(shape: IntArray, field: Field<T>) : NDField<T>(shape, field) {
-    override fun produceStructure(initializer: (IntArray) -> T): NDStructure<T> = genericNdStructure(shape, initializer)
+class GenericNDField<T : Any, F : Field<T>>(shape: IntArray, field: F) : NDField<T, F>(shape, field) {
+    override fun produceStructure(initializer: F.(IntArray) -> T): NDStructure<T> = genericNdStructure(shape) { field.initializer(it) }
 }
 
 //typealias NDFieldFactory<T> = (IntArray)->NDField<T>
@@ -173,23 +173,24 @@ object NDArrays {
     /**
      * Create a platform-optimized NDArray of doubles
      */
-    fun realNDArray(shape: IntArray, initializer: (IntArray) -> Double = { 0.0 }): NDElement<Double> {
+    fun realNDArray(shape: IntArray, initializer: DoubleField.(IntArray) -> Double = { 0.0 }): NDElement<Double, DoubleField> {
         return ExtendedNDField(shape, DoubleField).produce(initializer)
     }
 
-    fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDElement<Double> {
+    fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDElement<Double, DoubleField> {
         return realNDArray(intArrayOf(dim)) { initializer(it[0]) }
     }
 
-    fun real2DArray(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }): NDElement<Double> {
+    fun real2DArray(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }): NDElement<Double, DoubleField> {
         return realNDArray(intArrayOf(dim1, dim2)) { initializer(it[0], it[1]) }
     }
 
-    fun real3DArray(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }): NDElement<Double> {
+    fun real3DArray(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }): NDElement<Double, DoubleField> {
         return realNDArray(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
     }
 
-    inline fun produceReal(shape: IntArray, block: ExtendedNDField<Double>.() -> NDElement<Double>) = ExtendedNDField(shape, DoubleField).run(block)
+    inline fun produceReal(shape: IntArray, block: ExtendedNDField<Double, DoubleField>.() -> NDElement<Double, DoubleField>) =
+            ExtendedNDField(shape, DoubleField).run(block)
 
 //    /**
 //     * Simple boxing NDField
@@ -199,7 +200,7 @@ object NDArrays {
     /**
      * Simple boxing NDArray
      */
-    fun <T : Any> create(field: Field<T>, shape: IntArray, initializer: (IntArray) -> T): NDElement<T> {
+    fun <T : Any, F : Field<T>> create(field: F, shape: IntArray, initializer: (IntArray) -> T): NDElement<T, F> {
         return GenericNDField(shape, field).produce { initializer(it) }
     }
 }
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
index cd58a63b4..990adb0c7 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt
@@ -32,7 +32,7 @@ class NumberNDFieldTest {
         for (i in 0..2) {
             for (j in 0..2) {
                 val expected = (i * 10 + j).toDouble()
-                assertEquals(expected, array[i, j],"Error at index [$i, $j]")
+                assertEquals(expected, array[i, j], "Error at index [$i, $j]")
             }
         }
     }
@@ -41,25 +41,25 @@ class NumberNDFieldTest {
     fun testExternalFunction() {
         val function: (Double) -> Double = { x -> x.pow(2) + 2 * x + 1 }
         val result = function(array1) + 1.0
-        assertEquals(10.0, result[1,1])
+        assertEquals(10.0, result[1, 1])
     }
 
     @Test
     fun testLibraryFunction() {
         val abs: (Double) -> Double = ::abs
         val result = abs(array2)
-        assertEquals(2.0, result[0,2])
+        assertEquals(2.0, result[0, 2])
     }
 
-    object L2Norm: Norm<NDElement<out Number>, Double> {
-        override fun norm(arg: NDElement<out Number>): Double {
+    object L2Norm : Norm<NDElement<out Number, *>, Double> {
+        override fun norm(arg: NDElement<out Number, *>): Double {
             return kotlin.math.sqrt(arg.sumByDouble { it.second.toDouble() })
         }
     }
 
     @Test
-    fun testInternalContext(){
-        produceReal(array1.shape){
+    fun testInternalContext() {
+        produceReal(array1.shape) {
             with(L2Norm) {
                 1 + norm(array1) + exp(array2)
             }

From 15cc4a22e273f8fbff7a88b71ad09ac15d89b223 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Fri, 23 Nov 2018 21:44:42 +0300
Subject: [PATCH 10/13] Fixed some performance issues with histogram

---
 .../kmath/histogram/FastHistogram.kt          | 56 +++++++++++--------
 .../kotlin/scientifik/kmath/linear/Matrix.kt  | 22 +++++++-
 .../kotlin/scientifik/kmath/linear/Vector.kt  | 16 ------
 .../scientifik/kmath/structures/Buffers.kt    | 23 +++++++-
 .../kmath/structures/NDStructure.kt           |  4 +-
 5 files changed, 74 insertions(+), 47 deletions(-)

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
index 37c4eed2a..f9388f668 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
@@ -1,24 +1,31 @@
 package scientifik.kmath.histogram
 
-import scientifik.kmath.linear.RealVector
 import scientifik.kmath.linear.toVector
 import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.ListBuffer
 import scientifik.kmath.structures.NDStructure
 import scientifik.kmath.structures.ndStructure
 import kotlin.math.floor
 
-class MultivariateBin(override val center: RealVector, val sizes: RealVector, val counter: LongCounter = LongCounter()) : Bin<Double> {
+typealias RealPoint = Point<Double>
+
+private operator fun RealPoint.minus(other: RealPoint) = ListBuffer((0 until size).map { get(it) - other[it] })
+
+private inline fun <T> Buffer<out Double>.mapIndexed(crossinline mapper: (Int, Double) -> T): Sequence<T> = (0 until size).asSequence().map { mapper(it, get(it)) }
+
+
+class MultivariateBin(override val center: RealPoint, val sizes: RealPoint, var counter: Long = 0) : Bin<Double> {
     init {
         if (center.size != sizes.size) error("Dimension mismatch in bin creation. Expected ${center.size}, but found ${sizes.size}")
     }
 
     override fun contains(vector: Buffer<out Double>): Boolean {
         if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}")
-        return vector.asSequence().mapIndexed { i, value -> value in (center[i] - sizes[i] / 2)..(center[i] + sizes[i] / 2) }.all { it }
+        return vector.mapIndexed { i, value -> value in (center[i] - sizes[i] / 2)..(center[i] + sizes[i] / 2) }.all { it }
     }
 
-    override val value get() = counter.sum()
-    internal operator fun inc() = this.also { counter.increment() }
+    override val value get() = counter
+    internal operator fun inc() = this.also { counter++ }
 
     override val dimension: Int get() = center.size
 }
@@ -28,8 +35,8 @@ class MultivariateBin(override val center: RealVector, val sizes: RealVector, va
  * The histogram is optimized for speed, but have large size in memory
  */
 class FastHistogram(
-        private val lower: RealVector,
-        private val upper: RealVector,
+        private val lower: RealPoint,
+        private val upper: RealPoint,
         private val binNums: IntArray = IntArray(lower.size) { 20 }
 ) : MutableHistogram<Double, MultivariateBin> {
 
@@ -37,25 +44,27 @@ class FastHistogram(
         // argument checks
         if (lower.size != upper.size) error("Dimension mismatch in histogram lower and upper limits.")
         if (lower.size != binNums.size) error("Dimension mismatch in bin count.")
-        if ((upper - lower).any { it <= 0 }) error("Range for one of axis is not strictly positive")
+        if ((upper - lower).asSequence().any { it <= 0 }) error("Range for one of axis is not strictly positive")
     }
 
 
     override val dimension: Int get() = lower.size
 
     //TODO optimize binSize performance if needed
-    private val binSize = (upper - lower).mapIndexed { index, value -> value / binNums[index] }.toVector()
+    private val binSize: RealPoint = ListBuffer((upper - lower).mapIndexed { index, value -> value / binNums[index] }.toList())
 
     private val bins: NDStructure<MultivariateBin> by lazy {
         val actualSizes = IntArray(binNums.size) { binNums[it] + 2 }
         ndStructure(actualSizes) { indexArray ->
-            val center = indexArray.mapIndexed { axis, index ->
-                when (index) {
-                    0 -> Double.NEGATIVE_INFINITY
-                    actualSizes[axis] - 1 -> Double.POSITIVE_INFINITY
-                    else -> lower[axis] + (index.toDouble() - 0.5) * binSize[axis]
-                }
-            }.toVector()
+            val center = ListBuffer(
+                    indexArray.mapIndexed { axis, index ->
+                        when (index) {
+                            0 -> Double.NEGATIVE_INFINITY
+                            actualSizes[axis] - 1 -> Double.POSITIVE_INFINITY
+                            else -> lower[axis] + (index.toDouble() - 0.5) * binSize[axis]
+                        }
+                    }
+            )
             MultivariateBin(center, binSize)
         }
     }
@@ -91,13 +100,14 @@ class FastHistogram(
         return ndStructure(this.bins.shape) { bins[it].value }
     }
 
-    /**
-     * Create a phantom lightweight immutable copy of this histogram
-     */
-    fun asPhantom(): PhantomHistogram<Double> {
-        val binTemplates = bins.associate { (index, bin) -> BinTemplate<Double>(bin.center, bin.sizes) to index }
-        return PhantomHistogram(binTemplates, asND())
-    }
+//    /**
+//     * Create a phantom lightweight immutable copy of this histogram
+//     */
+//    fun asPhantom(): PhantomHistogram<Double> {
+//        val center =
+//        val binTemplates = bins.associate { (index, bin) -> BinTemplate<Double>(bin.center, bin.sizes) to index }
+//        return PhantomHistogram(binTemplates, asND())
+//    }
 
     companion object {
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
index f7ae04849..243eb2d6c 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt
@@ -1,9 +1,7 @@
 package scientifik.kmath.linear
 
 import scientifik.kmath.operations.*
-import scientifik.kmath.structures.ExtendedNDField
-import scientifik.kmath.structures.GenericNDField
-import scientifik.kmath.structures.NDField
+import scientifik.kmath.structures.*
 
 /**
  * The space for linear elements. Supports scalar product alongside with standard linear operations.
@@ -169,3 +167,21 @@ class ArrayMatrixSpace<T : Any, F : Field<T>>(
         return ArrayMatrixSpace(rows, columns, field, ndFactory)
     }
 }
+
+/**
+ * Member of [ArrayMatrixSpace] which wraps 2-D array
+ */
+class ArrayMatrix<T : Any, F : Field<T>> internal constructor(override val context: ArrayMatrixSpace<T, F>, val element: NDElement<T, F>) : Matrix<T, F> {
+
+    constructor(context: ArrayMatrixSpace<T, F>, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) })
+
+    override val rows: Int get() = context.rows
+
+    override val columns: Int get() = context.columns
+
+    override fun get(i: Int, j: Int): T {
+        return element[i, j]
+    }
+
+    override val self: ArrayMatrix<T, F> get() = this
+}
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
index 3e0f91ebe..38bc54ab6 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
@@ -73,23 +73,7 @@ class ArrayVectorSpace<T : Any, F : Field<T>>(
     override fun produce(initializer: (Int) -> T): Vector<T, F> = ArrayVector(this, initializer)
 }
 
-/**
- * Member of [ArrayMatrixSpace] which wraps 2-D array
- */
-class ArrayMatrix<T : Any, F : Field<T>> internal constructor(override val context: ArrayMatrixSpace<T, F>, val element: NDElement<T, F>) : Matrix<T, F> {
 
-    constructor(context: ArrayMatrixSpace<T, F>, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) })
-
-    override val rows: Int get() = context.rows
-
-    override val columns: Int get() = context.columns
-
-    override fun get(i: Int, j: Int): T {
-        return element[i, j]
-    }
-
-    override val self: ArrayMatrix<T, F> get() = this
-}
 
 
 class ArrayVector<T : Any, F : Field<T>> internal constructor(override val context: VectorSpace<T, F>, val element: NDElement<T, F>) : Vector<T, F> {
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
index 052f9caa4..5b80fa135 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
@@ -4,12 +4,16 @@ package scientifik.kmath.structures
 /**
  * A generic random access structure for both primitives and objects
  */
-interface Buffer<T> : Iterable<T> {
+interface Buffer<T> {
 
     val size: Int
 
     operator fun get(index: Int): T
 
+    operator fun iterator(): Iterator<T>
+
+    fun asSequence(): Sequence<T> = iterator().asSequence()
+
     /**
      * A shallow copy of the buffer
      */
@@ -25,7 +29,20 @@ interface MutableBuffer<T> : Buffer<T> {
     override fun copy(): MutableBuffer<T>
 }
 
-inline class ListBuffer<T>(private val list: MutableList<T>) : MutableBuffer<T> {
+
+inline class ListBuffer<T>(private val list: List<T>) : Buffer<T> {
+
+    override val size: Int
+        get() = list.size
+
+    override fun get(index: Int): T = list[index]
+
+    override fun iterator(): Iterator<T>  = list.iterator()
+
+    override fun copy(): ListBuffer<T> = ListBuffer(ArrayList(list))
+}
+
+inline class MutableListBuffer<T>(private val list: MutableList<T>) : MutableBuffer<T> {
 
     override val size: Int
         get() = list.size
@@ -38,7 +55,7 @@ inline class ListBuffer<T>(private val list: MutableList<T>) : MutableBuffer<T>
 
     override fun iterator(): Iterator<T>  = list.iterator()
 
-    override fun copy(): MutableBuffer<T> = ListBuffer(ArrayList(list))
+    override fun copy(): MutableBuffer<T> = MutableListBuffer(ArrayList(list))
 }
 
 class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
index a6880362d..ed05613c5 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
@@ -76,7 +76,7 @@ class DefaultStrides(override val shape: IntArray) : Strides {
     override fun offset(index: IntArray): Int {
         return index.mapIndexed { i, value ->
             if (value < 0 || value >= shape[i]) {
-                throw RuntimeException("Index $value out of shape bounds: (0,${shape[i]})")
+                throw RuntimeException("Index $value out of shape bounds: (0,${this.shape[i]})")
             }
             value * strides[i]
         }.sum()
@@ -169,6 +169,6 @@ fun <T> genericNdStructure(shape: IntArray, initializer: (IntArray) -> T): Mutab
             yield(initializer(it))
         }
     }
-    val buffer = ListBuffer<T>(sequence.toMutableList())
+    val buffer = MutableListBuffer<T>(sequence.toMutableList())
     return MutableBufferNDStructure(strides, buffer)
 }

From d799668083af8cbd32be58bde37b8b78d6eb8d55 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sat, 24 Nov 2018 20:55:35 +0300
Subject: [PATCH 11/13] Faster implementation for FastHistogram

---
 .../kmath/histogram/FastHistogram.kt          | 107 ++++++++++--------
 .../scientifik/kmath/histogram/Histogram.kt   |   2 +
 .../kmath/histogram/PhantomHistogram.kt       |   7 +-
 .../kotlin/scientifik/kmath/linear/Vector.kt  |  14 ++-
 .../scientifik/kmath/structures/Buffers.kt    |   4 +-
 5 files changed, 77 insertions(+), 57 deletions(-)

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
index f9388f668..37f46b1d7 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
@@ -1,44 +1,48 @@
 package scientifik.kmath.histogram
 
 import scientifik.kmath.linear.toVector
-import scientifik.kmath.structures.Buffer
-import scientifik.kmath.structures.ListBuffer
-import scientifik.kmath.structures.NDStructure
-import scientifik.kmath.structures.ndStructure
+import scientifik.kmath.structures.*
 import kotlin.math.floor
 
-typealias RealPoint = Point<Double>
-
 private operator fun RealPoint.minus(other: RealPoint) = ListBuffer((0 until size).map { get(it) - other[it] })
 
 private inline fun <T> Buffer<out Double>.mapIndexed(crossinline mapper: (Int, Double) -> T): Sequence<T> = (0 until size).asSequence().map { mapper(it, get(it)) }
 
 
-class MultivariateBin(override val center: RealPoint, val sizes: RealPoint, var counter: Long = 0) : Bin<Double> {
-    init {
-        if (center.size != sizes.size) error("Dimension mismatch in bin creation. Expected ${center.size}, but found ${sizes.size}")
-    }
-
-    override fun contains(vector: Buffer<out Double>): Boolean {
-        if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}")
-        return vector.mapIndexed { i, value -> value in (center[i] - sizes[i] / 2)..(center[i] + sizes[i] / 2) }.all { it }
-    }
-
-    override val value get() = counter
-    internal operator fun inc() = this.also { counter++ }
-
-    override val dimension: Int get() = center.size
-}
+//class MultivariateBin(override val center: RealPoint, val sizes: RealPoint, var counter: Long = 0) : Bin<Double> {
+//    init {
+//        if (center.size != sizes.size) error("Dimension mismatch in bin creation. Expected ${center.size}, but found ${sizes.size}")
+//    }
+//
+//    override fun contains(vector: Buffer<out Double>): Boolean {
+//        if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}")
+//        return vector.mapIndexed { i, value -> value in (center[i] - sizes[i] / 2)..(center[i] + sizes[i] / 2) }.all { it }
+//    }
+//
+//    override val value get() = counter
+//    internal operator fun inc() = this.also { counter++ }
+//
+//    override val dimension: Int get() = center.size
+//}
 
 /**
  * Uniform multivariate histogram with fixed borders. Based on NDStructure implementation with complexity of m for bin search, where m is the number of dimensions.
- * The histogram is optimized for speed, but have large size in memory
  */
 class FastHistogram(
         private val lower: RealPoint,
         private val upper: RealPoint,
         private val binNums: IntArray = IntArray(lower.size) { 20 }
-) : MutableHistogram<Double, MultivariateBin> {
+) : MutableHistogram<Double, PhantomBin<Double>> {
+
+
+    private val strides = DefaultStrides(IntArray(binNums.size) { binNums[it] + 2 })
+
+    private val values: NDStructure<LongCounter> = ndStructure(strides) { LongCounter() }
+
+    //private val weight: NDStructure<DoubleCounter?> = ndStructure(strides){null}
+
+    //TODO optimize binSize performance if needed
+    private val binSize: RealPoint = ListBuffer((upper - lower).mapIndexed { index, value -> value / binNums[index] }.toList())
 
     init {
         // argument checks
@@ -50,24 +54,6 @@ class FastHistogram(
 
     override val dimension: Int get() = lower.size
 
-    //TODO optimize binSize performance if needed
-    private val binSize: RealPoint = ListBuffer((upper - lower).mapIndexed { index, value -> value / binNums[index] }.toList())
-
-    private val bins: NDStructure<MultivariateBin> by lazy {
-        val actualSizes = IntArray(binNums.size) { binNums[it] + 2 }
-        ndStructure(actualSizes) { indexArray ->
-            val center = ListBuffer(
-                    indexArray.mapIndexed { axis, index ->
-                        when (index) {
-                            0 -> Double.NEGATIVE_INFINITY
-                            actualSizes[axis] - 1 -> Double.POSITIVE_INFINITY
-                            else -> lower[axis] + (index.toDouble() - 0.5) * binSize[axis]
-                        }
-                    }
-            )
-            MultivariateBin(center, binSize)
-        }
-    }
 
     /**
      * Get internal [NDStructure] bin index for given axis
@@ -80,24 +66,51 @@ class FastHistogram(
         }
     }
 
+    private fun getIndex(point: Buffer<out Double>): IntArray = IntArray(dimension) { getIndex(it, point[it]) }
 
-    override fun get(point: Buffer<out Double>): MultivariateBin? {
-        val index = IntArray(dimension) { getIndex(it, point[it]) }
-        return bins[index]
+    private fun getValue(index: IntArray): Long {
+        return values[index].sum()
+    }
+
+    fun getValue(point: Buffer<out Double>): Long {
+        return getValue(getIndex(point))
+    }
+
+    private fun getTemplate(index: IntArray): BinTemplate<Double> {
+        val center = index.mapIndexed { axis, i ->
+            when (i) {
+                0 -> Double.NEGATIVE_INFINITY
+                strides.shape[axis] - 1 -> Double.POSITIVE_INFINITY
+                else -> lower[axis] + (i.toDouble() - 0.5) * binSize[axis]
+            }
+        }.toVector()
+        return BinTemplate(center, binSize)
+    }
+
+    fun getTemplate(point: Buffer<out Double>): BinTemplate<Double> {
+        return getTemplate(getIndex(point))
+    }
+
+    override fun get(point: Buffer<out Double>): PhantomBin<Double>? {
+        val index = getIndex(point)
+        return PhantomBin(getTemplate(index), getValue(index))
     }
 
     override fun put(point: Buffer<out Double>, weight: Double) {
         if (weight != 1.0) TODO("Implement weighting")
-        this[point]?.inc() ?: error("Could not find appropriate bin (should not be possible)")
+        val index = getIndex(point)
+        values[index].increment()
     }
 
-    override fun iterator(): Iterator<MultivariateBin> = bins.asSequence().map { it.second }.iterator()
+    override fun iterator(): Iterator<PhantomBin<Double>> = values.asSequence().map { (index, value) ->
+        PhantomBin(getTemplate(index), value.sum())
+    }.iterator()
 
     /**
      * Convert this histogram into NDStructure containing bin values but not bin descriptions
      */
     fun asND(): NDStructure<Number> {
-        return ndStructure(this.bins.shape) { bins[it].value }
+        return ndStructure(this.values.shape) { values[it].sum() }
     }
 
 //    /**
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
index 8e429087a..86110f140 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
@@ -6,6 +6,8 @@ import scientifik.kmath.structures.DoubleBuffer
 
 typealias Point<T> = Buffer<T>
 
+typealias RealPoint = Buffer<Double>
+
 /**
  * A simple geometric domain
  * TODO move to geometry module
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
index 25a1bf33a..ffffb0d7d 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/PhantomHistogram.kt
@@ -3,12 +3,13 @@ package scientifik.kmath.histogram
 import scientifik.kmath.linear.Vector
 import scientifik.kmath.operations.Space
 import scientifik.kmath.structures.NDStructure
+import scientifik.kmath.structures.asSequence
 
-data class BinTemplate<T : Comparable<T>>(val center: Vector<T, *>, val sizes: Vector<T, *>) {
+data class BinTemplate<T : Comparable<T>>(val center: Vector<T, *>, val sizes: Point<T>) {
     fun contains(vector: Point<out T>): Boolean {
         if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}")
-        val upper = center + sizes / 2.0
-        val lower = center - sizes / 2.0
+        val upper = center.context.run { center + sizes / 2.0}
+        val lower =  center.context.run {center - sizes / 2.0}
         return vector.asSequence().mapIndexed { i, value ->
             value in lower[i]..upper[i]
         }.all { it }
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
index 38bc54ab6..755d58d9b 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
@@ -39,18 +39,24 @@ interface Vector<T : Any, S : Space<T>> : SpaceElement<Point<T>, VectorSpace<T,
         /**
          * Create vector with custom field
          */
-        fun <T : Any, F: Field<T>> of(size: Int, field: F, initializer: (Int) -> T) =
+        fun <T : Any, F : Field<T>> of(size: Int, field: F, initializer: (Int) -> T) =
                 ArrayVector(ArrayVectorSpace(size, field), initializer)
 
+        private val realSpaceCache = HashMap<Int, ArrayVectorSpace<Double, DoubleField>>()
+
+        private fun getRealSpace(size: Int): ArrayVectorSpace<Double, DoubleField> {
+            return realSpaceCache.getOrPut(size){ArrayVectorSpace(size, DoubleField, realNDFieldFactory)}
+        }
+
         /**
          * Create vector of [Double]
          */
         fun ofReal(size: Int, initializer: (Int) -> Double) =
-                ArrayVector(ArrayVectorSpace(size, DoubleField, realNDFieldFactory), initializer)
+                ArrayVector(getRealSpace(size), initializer)
 
         fun ofReal(vararg point: Double) = point.toVector()
 
-        fun equals(v1: Vector<*,*>, v2: Vector<*,*>): Boolean {
+        fun equals(v1: Vector<*, *>, v2: Vector<*, *>): Boolean {
             if (v1 === v2) return true
             if (v1.context != v2.context) return false
             for (i in 0 until v2.size) {
@@ -74,8 +80,6 @@ class ArrayVectorSpace<T : Any, F : Field<T>>(
 }
 
 
-
-
 class ArrayVector<T : Any, F : Field<T>> internal constructor(override val context: VectorSpace<T, F>, val element: NDElement<T, F>) : Vector<T, F> {
 
     constructor(context: ArrayVectorSpace<T, F>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
index 5b80fa135..a4bb76060 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
@@ -12,14 +12,14 @@ interface Buffer<T> {
 
     operator fun iterator(): Iterator<T>
 
-    fun asSequence(): Sequence<T> = iterator().asSequence()
-
     /**
      * A shallow copy of the buffer
      */
     fun copy(): Buffer<T>
 }
 
+fun <T> Buffer<T>.asSequence(): Sequence<T> = iterator().asSequence()
+
 interface MutableBuffer<T> : Buffer<T> {
     operator fun set(index: Int, value: T)
 

From bdd9dccd4f37c27f3c33ecf563616f03d97eb8da Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sat, 24 Nov 2018 21:49:42 +0300
Subject: [PATCH 12/13] Final fixes for FastHistogram and NDStructure
 performance

---
 .../kmath/histogram/FastHistogram.kt          | 38 +++++--------------
 .../scientifik/kmath/histogram/Histogram.kt   |  1 +
 .../kmath/structures/NDStructure.kt           | 18 ++++-----
 3 files changed, 20 insertions(+), 37 deletions(-)

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
index 37f46b1d7..8b30e7ebf 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
@@ -8,23 +8,6 @@ private operator fun RealPoint.minus(other: RealPoint) = ListBuffer((0 until siz
 
 private inline fun <T> Buffer<out Double>.mapIndexed(crossinline mapper: (Int, Double) -> T): Sequence<T> = (0 until size).asSequence().map { mapper(it, get(it)) }
 
-
-//class MultivariateBin(override val center: RealPoint, val sizes: RealPoint, var counter: Long = 0) : Bin<Double> {
-//    init {
-//        if (center.size != sizes.size) error("Dimension mismatch in bin creation. Expected ${center.size}, but found ${sizes.size}")
-//    }
-//
-//    override fun contains(vector: Buffer<out Double>): Boolean {
-//        if (vector.size != center.size) error("Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}")
-//        return vector.mapIndexed { i, value -> value in (center[i] - sizes[i] / 2)..(center[i] + sizes[i] / 2) }.all { it }
-//    }
-//
-//    override val value get() = counter
-//    internal operator fun inc() = this.also { counter++ }
-//
-//    override val dimension: Int get() = center.size
-//}
-
 /**
  * Uniform multivariate histogram with fixed borders. Based on NDStructure implementation with complexity of m for bin search, where m is the number of dimensions.
  */
@@ -109,18 +92,17 @@ class FastHistogram(
     /**
      * Convert this histogram into NDStructure containing bin values but not bin descriptions
      */
-    fun asND(): NDStructure<Number> {
+    fun asNDStructure(): NDStructure<Number> {
         return ndStructure(this.values.shape) { values[it].sum() }
     }
 
-//    /**
-//     * Create a phantom lightweight immutable copy of this histogram
-//     */
-//    fun asPhantom(): PhantomHistogram<Double> {
-//        val center =
-//        val binTemplates = bins.associate { (index, bin) -> BinTemplate<Double>(bin.center, bin.sizes) to index }
-//        return PhantomHistogram(binTemplates, asND())
-//    }
+    /**
+     * Create a phantom lightweight immutable copy of this histogram
+     */
+    fun asPhantomHistogram(): PhantomHistogram<Double> {
+        val binTemplates = values.associate { (index, _) -> getTemplate(index) to index }
+        return PhantomHistogram(binTemplates, asNDStructure())
+    }
 
     companion object {
 
@@ -148,8 +130,8 @@ class FastHistogram(
          */
         fun fromRanges(vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>): FastHistogram {
             return FastHistogram(
-                    ranges.map { it.first.start }.toVector(),
-                    ranges.map { it.first.endInclusive }.toVector(),
+                    ListBuffer(ranges.map { it.first.start }),
+                    ListBuffer(ranges.map { it.first.endInclusive }),
                     ranges.map { it.second }.toIntArray()
             )
         }
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
index 86110f140..08214142e 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/Histogram.kt
@@ -53,6 +53,7 @@ interface MutableHistogram<T: Any, out B : Bin<T>>: Histogram<T,B>{
 fun <T: Any> MutableHistogram<T,*>.put(vararg point: T) = put(ArrayBuffer(point))
 
 fun MutableHistogram<Double,*>.put(vararg point: Number) = put(DoubleBuffer(point.map { it.toDouble() }.toDoubleArray()))
+fun MutableHistogram<Double,*>.put(vararg point: Double) = put(DoubleBuffer(point))
 
 fun <T: Any> MutableHistogram<T,*>.fill(sequence: Iterable<Point<T>>) = sequence.forEach { put(it) }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
index ed05613c5..ec985eb8b 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
@@ -83,15 +83,15 @@ class DefaultStrides(override val shape: IntArray) : Strides {
     }
 
     override fun index(offset: Int): IntArray {
-        return sequence {
-            var current = offset
-            var strideIndex = strides.size - 2
-            while (strideIndex >= 0) {
-                yield(current / strides[strideIndex])
-                current %= strides[strideIndex]
-                strideIndex--
-            }
-        }.toList().reversed().toIntArray()
+        val res = IntArray(shape.size)
+        var current = offset
+        var strideIndex = strides.size - 2
+        while (strideIndex >= 0) {
+            res[ strideIndex] = (current / strides[strideIndex])
+            current %= strides[strideIndex]
+            strideIndex--
+        }
+        return res
     }
 
     override val linearSize: Int

From cddea1869d0586f483caa90dff8429c2f3103354 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Tue, 11 Dec 2018 17:25:55 +0300
Subject: [PATCH 13/13] Lazy structures, structure is no longer iterable,
 removed copy on read-only buffers.

---
 README.md                                     |  4 +-
 build.gradle.kts                              |  9 ++-
 kmath-core/build.gradle                       |  8 +-
 .../kmath/histogram/FastHistogram.kt          |  4 +-
 .../kotlin/scientifik/kmath/linear/Vector.kt  |  2 -
 .../scientifik/kmath/operations/Complex.kt    |  4 +
 .../scientifik/kmath/operations/Fields.kt     | 15 ++++
 .../scientifik/kmath/structures/Buffers.kt    | 60 ++++++++++----
 .../scientifik/kmath/structures/NDField.kt    | 33 ++++----
 .../kmath/structures/NDStructure.kt           | 32 ++++----
 .../scientifik/kmath/structures/BufferSpec.kt | 55 +++++++++++++
 .../kmath/structures/ComplexBufferSpec.kt     | 25 ++++++
 .../kmath/structures/ObjectBuffer.kt          | 28 +++++++
 .../kmath/structures/RealBufferSpec.kt        | 23 ++++++
 .../kmath/structures/ComplexBufferSpecTest.kt | 17 ++++
 kmath-coroutines/build.gradle                 | 42 ++++++++++
 .../kmath/structures/CoroutinesExtra.kt       | 11 +++
 .../kmath/structures/LazyNDField.kt           | 79 +++++++++++++++++++
 .../kmath/structures/LazyNDFieldTest.kt       | 20 +++++
 .../kmath/structures/_CoroutinesExtra.kt      |  6 ++
 kmath-io/build.gradle                         |  4 +-
 kmath-jmh/build.gradle                        |  4 +-
 .../kmath/structures/ArrayBenchmark.kt        |  6 +-
 .../kmath/structures/BufferBenchmark.kt       | 38 +++++++++
 settings.gradle.kts                           |  4 +-
 25 files changed, 465 insertions(+), 68 deletions(-)
 create mode 100644 kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/BufferSpec.kt
 create mode 100644 kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt
 create mode 100644 kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ObjectBuffer.kt
 create mode 100644 kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/RealBufferSpec.kt
 create mode 100644 kmath-core/src/jvmTest/kotlin/scientifik/kmath/structures/ComplexBufferSpecTest.kt
 create mode 100644 kmath-coroutines/build.gradle
 create mode 100644 kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/CoroutinesExtra.kt
 create mode 100644 kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/LazyNDField.kt
 create mode 100644 kmath-coroutines/src/commonTest/kotlin/scientifik/kmath/structures/LazyNDFieldTest.kt
 create mode 100644 kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/_CoroutinesExtra.kt
 create mode 100644 kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/BufferBenchmark.kt

diff --git a/README.md b/README.md
index 607cda1bf..8e07b546a 100644
--- a/README.md
+++ b/README.md
@@ -58,11 +58,11 @@ repositories {
 
 Then use regular dependency like
 ```groovy
-compile(group: 'scientifik', name: 'kmath-core-jvm', version: '0.0.1-SNAPSHOT')
+compile(group: 'scientifik', name: 'kmath-core', version: '0.0.1-SNAPSHOT')
 ```
 or in kotlin
 ```kotlin
-compile(group = "scientifik", name = "kmath-core-jvm", version = "0.0.1-SNAPSHOT")
+compile(group = "scientifik", name = "kmath-core", version = "0.0.1-SNAPSHOT")
 ```
 
 Work builds could be obtained with [![](https://jitpack.io/v/altavir/kmath.svg)](https://jitpack.io/#altavir/kmath). 
diff --git a/build.gradle.kts b/build.gradle.kts
index 8f8bcaa34..5b95be71a 100644
--- a/build.gradle.kts
+++ b/build.gradle.kts
@@ -1,7 +1,11 @@
 buildscript {
-    extra["kotlinVersion"] = "1.3.10"
+    extra["kotlinVersion"] = "1.3.11"
+    extra["ioVersion"] = "0.1.2-dev-2"
+    extra["coroutinesVersion"] = "1.0.1"
 
     val kotlinVersion: String by extra
+    val ioVersion: String by extra
+    val coroutinesVersion: String by extra
 
     repositories {
         jcenter()
@@ -15,6 +19,7 @@ buildscript {
 
 plugins {
     id("com.jfrog.artifactory") version "4.8.1"  apply false
+//    id("org.jetbrains.kotlin.multiplatform") apply false
 }
 
 allprojects {
@@ -22,7 +27,7 @@ allprojects {
     apply(plugin = "com.jfrog.artifactory")
 
     group = "scientifik"
-    version = "0.0.1-SNAPSHOT"
+    version = "0.0.2-dev-1"
 }
 
 if(file("artifactory.gradle").exists()){
diff --git a/kmath-core/build.gradle b/kmath-core/build.gradle
index c63e518fc..d51ccdb65 100644
--- a/kmath-core/build.gradle
+++ b/kmath-core/build.gradle
@@ -1,5 +1,5 @@
 plugins {
-    id 'kotlin-multiplatform'
+    id "org.jetbrains.kotlin.multiplatform"
 }
 
 kotlin {
@@ -14,7 +14,7 @@ kotlin {
     sourceSets {
         commonMain {
             dependencies {
-                implementation 'org.jetbrains.kotlin:kotlin-stdlib-common'
+                api 'org.jetbrains.kotlin:kotlin-stdlib-common'
             }
         }
         commonTest {
@@ -25,7 +25,7 @@ kotlin {
         }
         jvmMain {
             dependencies {
-                implementation 'org.jetbrains.kotlin:kotlin-stdlib-jdk8'
+                api 'org.jetbrains.kotlin:kotlin-stdlib-jdk8'
             }
         }
         jvmTest {
@@ -36,7 +36,7 @@ kotlin {
         }
         jsMain {
             dependencies {
-                implementation 'org.jetbrains.kotlin:kotlin-stdlib-js'
+                api 'org.jetbrains.kotlin:kotlin-stdlib-js'
             }
         }
         jsTest {
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
index 8b30e7ebf..e48979430 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt
@@ -85,7 +85,7 @@ class FastHistogram(
         values[index].increment()
     }
 
-    override fun iterator(): Iterator<PhantomBin<Double>> = values.asSequence().map { (index, value) ->
+    override fun iterator(): Iterator<PhantomBin<Double>> = values.elements().map { (index, value) ->
         PhantomBin(getTemplate(index), value.sum())
     }.iterator()
 
@@ -100,7 +100,7 @@ class FastHistogram(
      * Create a phantom lightweight immutable copy of this histogram
      */
     fun asPhantomHistogram(): PhantomHistogram<Double> {
-        val binTemplates = values.associate { (index, _) -> getTemplate(index) to index }
+        val binTemplates = values.elements().associate { (index, _) -> getTemplate(index) to index }
         return PhantomHistogram(binTemplates, asNDStructure())
     }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
index 755d58d9b..5ab9907de 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt
@@ -98,8 +98,6 @@ class ArrayVector<T : Any, F : Field<T>> internal constructor(override val conte
 
     override fun iterator(): Iterator<T> = (0 until size).map { element[it] }.iterator()
 
-    override fun copy(): ArrayVector<T, F> = ArrayVector(context, element)
-
     override fun toString(): String = this.joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() }
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt
index dd3be256e..8d3c01d28 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt
@@ -38,4 +38,8 @@ data class Complex(val re: Double, val im: Double) : FieldElement<Complex, Compl
     val abs: Double
         get() = kotlin.math.sqrt(square)
 
+    companion object {
+
+    }
+
 }
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
index 3e4f47869..d993d37b4 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
@@ -49,6 +49,9 @@ inline class Real(val value: Double) : FieldElement<Real, RealField> {
     override val context
         get() = RealField
 
+    companion object {
+
+    }
 }
 
 /**
@@ -71,4 +74,16 @@ object DoubleField : ExtendedField<Double>, Norm<Double, Double> {
     override fun ln(arg: Double): Double = kotlin.math.ln(arg)
 
     override fun norm(arg: Double): Double = kotlin.math.abs(arg)
+}
+
+/**
+ * A field for double without boxing. Does not produce appropriate field element
+ */
+object IntField : Field<Int>{
+    override val zero: Int = 0
+    override fun add(a: Int, b: Int): Int = a + b
+    override fun multiply(a: Int, b: Int): Int = a * b
+    override fun multiply(a: Int, k: Double): Int = (k*a).toInt()
+    override val one: Int = 1
+    override fun divide(a: Int, b: Int): Int = a / b
 }
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
index a4bb76060..69e8163c6 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt
@@ -11,11 +11,6 @@ interface Buffer<T> {
     operator fun get(index: Int): T
 
     operator fun iterator(): Iterator<T>
-
-    /**
-     * A shallow copy of the buffer
-     */
-    fun copy(): Buffer<T>
 }
 
 fun <T> Buffer<T>.asSequence(): Sequence<T> = iterator().asSequence()
@@ -26,7 +21,7 @@ interface MutableBuffer<T> : Buffer<T> {
     /**
      * A shallow copy of the buffer
      */
-    override fun copy(): MutableBuffer<T>
+    fun copy(): MutableBuffer<T>
 }
 
 
@@ -37,9 +32,7 @@ inline class ListBuffer<T>(private val list: List<T>) : Buffer<T> {
 
     override fun get(index: Int): T = list[index]
 
-    override fun iterator(): Iterator<T>  = list.iterator()
-
-    override fun copy(): ListBuffer<T> = ListBuffer(ArrayList(list))
+    override fun iterator(): Iterator<T> = list.iterator()
 }
 
 inline class MutableListBuffer<T>(private val list: MutableList<T>) : MutableBuffer<T> {
@@ -53,7 +46,7 @@ inline class MutableListBuffer<T>(private val list: MutableList<T>) : MutableBuf
         list[index] = value
     }
 
-    override fun iterator(): Iterator<T>  = list.iterator()
+    override fun iterator(): Iterator<T> = list.iterator()
 
     override fun copy(): MutableBuffer<T> = MutableListBuffer(ArrayList(list))
 }
@@ -69,7 +62,7 @@ class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
         array[index] = value
     }
 
-    override fun iterator(): Iterator<T>  = array.iterator()
+    override fun iterator(): Iterator<T> = array.iterator()
 
     override fun copy(): MutableBuffer<T> = ArrayBuffer(array.copyOf())
 }
@@ -84,12 +77,12 @@ inline class DoubleBuffer(private val array: DoubleArray) : MutableBuffer<Double
         array[index] = value
     }
 
-    override fun iterator(): Iterator<Double>  = array.iterator()
+    override fun iterator(): Iterator<Double> = array.iterator()
 
     override fun copy(): MutableBuffer<Double> = DoubleBuffer(array.copyOf())
 }
 
-inline class IntBuffer(private val array: IntArray): MutableBuffer<Int>{
+inline class IntBuffer(private val array: IntArray) : MutableBuffer<Int> {
     override val size: Int
         get() = array.size
 
@@ -99,15 +92,48 @@ inline class IntBuffer(private val array: IntArray): MutableBuffer<Int>{
         array[index] = value
     }
 
-    override fun iterator(): Iterator<Int>  = array.iterator()
+    override fun iterator(): Iterator<Int> = array.iterator()
 
     override fun copy(): MutableBuffer<Int> = IntBuffer(array.copyOf())
 }
 
-inline fun <reified T : Any> buffer(size: Int, noinline initializer: (Int) -> T): Buffer<T> {
-    return ArrayBuffer(Array(size, initializer))
+inline class ReadOnlyBuffer<T>(private val buffer: MutableBuffer<T>) : Buffer<T> {
+    override val size: Int get() = buffer.size
+
+    override fun get(index: Int): T = buffer.get(index)
+
+    override fun iterator(): Iterator<T> = buffer.iterator()
 }
 
+/**
+ * Convert this buffer to read-only buffer
+ */
+fun <T> Buffer<T>.asReadOnly(): Buffer<T> = if (this is MutableBuffer) {
+    ReadOnlyBuffer(this)
+} else {
+    this
+}
+
+/**
+ * Create most appropriate immutable buffer for given type avoiding boxing wherever possible
+ */
+@Suppress("UNCHECKED_CAST")
+inline fun <reified T : Any> buffer(size: Int, noinline initializer: (Int) -> T): Buffer<T> {
+    return when (T::class) {
+        Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as Buffer<T>
+        Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as Buffer<T>
+        else -> ArrayBuffer(Array(size, initializer))
+    }
+}
+
+/**
+ * Create most appropriate mutable buffer for given type avoiding boxing wherever possible
+ */
+@Suppress("UNCHECKED_CAST")
 inline fun <reified T : Any> mutableBuffer(size: Int, noinline initializer: (Int) -> T): MutableBuffer<T> {
-    return ArrayBuffer(Array(size, initializer))
+    return when (T::class) {
+        Double::class -> DoubleBuffer(DoubleArray(size) { initializer(it) as Double }) as MutableBuffer<T>
+        Int::class -> IntBuffer(IntArray(size) { initializer(it) as Int }) as MutableBuffer<T>
+        else -> ArrayBuffer(Array(size, initializer))
+    }
 }
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
index d8e9cc31b..cddac5801 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDField.kt
@@ -23,7 +23,7 @@ abstract class NDField<T, F : Field<T>>(val shape: IntArray, val field: F) : Fie
      * Create new instance of NDArray using field shape and given initializer
      * The producer takes list of indices as argument and returns contained value
      */
-    fun produce(initializer: F.(IntArray) -> T): NDElement<T, F> = NDElement(this, produceStructure(initializer))
+    fun produce(initializer: F.(IntArray) -> T): NDElement<T, F> = NDStructureElement(this, produceStructure(initializer))
 
     override val zero: NDElement<T, F> by lazy {
         produce { zero }
@@ -83,7 +83,7 @@ abstract class NDField<T, F : Field<T>>(val shape: IntArray, val field: F) : Fie
 //    /**
 //     * Reverse minus operation
 //     */
-//    operator fun T.minus(arg: NDElement<T, F>): NDElement<T, F> = arg.transform { _, value ->
+//    operator fun T.minus(arg: NDElement<T, F>): NDElement<T, F> = arg.transformIndexed { _, value ->
 //        with(arg.context.field) {
 //            this@minus - value
 //        }
@@ -97,38 +97,41 @@ abstract class NDField<T, F : Field<T>>(val shape: IntArray, val field: F) : Fie
 //    /**
 //     * Reverse division operation
 //     */
-//    operator fun T.div(arg: NDElement<T, F>): NDElement<T, F> = arg.transform { _, value ->
+//    operator fun T.div(arg: NDElement<T, F>): NDElement<T, F> = arg.transformIndexed { _, value ->
 //        with(arg.context.field) {
 //            this@div / value
 //        }
 //    }
 }
 
+
+interface NDElement<T, F : Field<T>>: FieldElement<NDElement<T, F>, NDField<T, F>>, NDStructure<T>
+
+inline fun <T, F : Field<T>> NDElement<T, F>.transformIndexed(crossinline action: F.(IntArray, T) -> T): NDElement<T, F> = context.produce { action(it, get(*it)) }
+inline fun <T, F : Field<T>> NDElement<T, F>.transform(crossinline action: F.(T) -> T): NDElement<T, F> = context.produce { action(get(*it)) }
+
+
 /**
- *  Immutable [NDStructure] coupled to the context. Emulates Python ndarray
+ *  Read-only [NDStructure] coupled to the context.
  */
-class NDElement<T, F : Field<T>>(override val context: NDField<T, F>, private val structure: NDStructure<T>) : FieldElement<NDElement<T, F>, NDField<T, F>>, NDStructure<T> by structure {
+class NDStructureElement<T, F : Field<T>>(override val context: NDField<T, F>, private val structure: NDStructure<T>) : NDElement<T,F>, NDStructure<T> by structure {
 
     //TODO ensure structure is immutable
 
-    override val self: NDElement<T, F>
-        get() = this
-
-    inline fun transform(crossinline action: (IntArray, T) -> T): NDElement<T, F> = context.produce { action(it, get(*it)) }
-    inline fun transform(crossinline action: (T) -> T): NDElement<T, F> = context.produce { action(get(*it)) }
+    override val self: NDElement<T, F>  get() = this
 }
 
 /**
  * Element by element application of any operation on elements to the whole array. Just like in numpy
  */
-operator fun <T, F : Field<T>> Function1<T, T>.invoke(ndElement: NDElement<T, F>): NDElement<T, F> = ndElement.transform { _, value -> this(value) }
+operator fun <T, F : Field<T>> Function1<T, T>.invoke(ndElement: NDElement<T, F>): NDElement<T, F> = ndElement.transform {value -> this@invoke(value) }
 
 /* plus and minus */
 
 /**
  * Summation operation for [NDElement] and single element
  */
-operator fun <T, F : Field<T>> NDElement<T, F>.plus(arg: T): NDElement<T, F> = transform { _, value ->
+operator fun <T, F : Field<T>> NDElement<T, F>.plus(arg: T): NDElement<T, F> = transform {value ->
     with(context.field) {
         arg + value
     }
@@ -137,7 +140,7 @@ operator fun <T, F : Field<T>> NDElement<T, F>.plus(arg: T): NDElement<T, F> = t
 /**
  * Subtraction operation between [NDElement] and single element
  */
-operator fun <T, F : Field<T>> NDElement<T, F>.minus(arg: T): NDElement<T, F> = transform { _, value ->
+operator fun <T, F : Field<T>> NDElement<T, F>.minus(arg: T): NDElement<T, F> = transform {value ->
     with(context.field) {
         arg - value
     }
@@ -148,7 +151,7 @@ operator fun <T, F : Field<T>> NDElement<T, F>.minus(arg: T): NDElement<T, F> =
 /**
  * Product operation for [NDElement] and single element
  */
-operator fun <T, F : Field<T>> NDElement<T, F>.times(arg: T): NDElement<T, F> = transform { _, value ->
+operator fun <T, F : Field<T>> NDElement<T, F>.times(arg: T): NDElement<T, F> = transform { value ->
     with(context.field) {
         arg * value
     }
@@ -157,7 +160,7 @@ operator fun <T, F : Field<T>> NDElement<T, F>.times(arg: T): NDElement<T, F> =
 /**
  * Division operation between [NDElement] and single element
  */
-operator fun <T, F : Field<T>> NDElement<T, F>.div(arg: T): NDElement<T, F> = transform { _, value ->
+operator fun <T, F : Field<T>> NDElement<T, F>.div(arg: T): NDElement<T, F> = transform { value ->
     with(context.field) {
         arg / value
     }
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
index ec985eb8b..f9521cf93 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt
@@ -1,7 +1,7 @@
 package scientifik.kmath.structures
 
 
-interface NDStructure<T> : Iterable<Pair<IntArray, T>> {
+interface NDStructure<T> {
 
     val shape: IntArray
 
@@ -9,6 +9,8 @@ interface NDStructure<T> : Iterable<Pair<IntArray, T>> {
         get() = shape.size
 
     operator fun get(index: IntArray): T
+
+    fun elements(): Sequence<Pair<IntArray, T>>
 }
 
 operator fun <T> NDStructure<T>.get(vararg index: Int): T = get(index)
@@ -18,7 +20,7 @@ interface MutableNDStructure<T> : NDStructure<T> {
 }
 
 fun <T> MutableNDStructure<T>.transformInPlace(action: (IntArray, T) -> T) {
-    for ((index, oldValue) in this) {
+    elements().forEach { (index, oldValue) ->
         this[index] = action(index, oldValue)
     }
 }
@@ -87,7 +89,7 @@ class DefaultStrides(override val shape: IntArray) : Strides {
         var current = offset
         var strideIndex = strides.size - 2
         while (strideIndex >= 0) {
-            res[ strideIndex] = (current / strides[strideIndex])
+            res[strideIndex] = (current / strides[strideIndex])
             current %= strides[strideIndex]
             strideIndex--
         }
@@ -107,8 +109,8 @@ abstract class GenericNDStructure<T, B : Buffer<T>> : NDStructure<T> {
     override val shape: IntArray
         get() = strides.shape
 
-    override fun iterator(): Iterator<Pair<IntArray, T>> =
-            strides.indices().map { it to this[it] }.iterator()
+    override fun elements()=
+            strides.indices().map { it to this[it] }
 }
 
 /**
@@ -126,10 +128,10 @@ class BufferNDStructure<T>(
     }
 }
 
-inline fun <reified T: Any> ndStructure(strides: Strides, noinline initializer: (IntArray) -> T) =
-        BufferNDStructure<T>(strides, buffer(strides.linearSize){ i-> initializer(strides.index(i))})
+inline fun <reified T : Any> ndStructure(strides: Strides, noinline initializer: (IntArray) -> T) =
+        BufferNDStructure<T>(strides, buffer(strides.linearSize) { i -> initializer(strides.index(i)) })
 
-inline fun <reified T: Any> ndStructure(shape: IntArray, noinline initializer: (IntArray) -> T) =
+inline fun <reified T : Any> ndStructure(shape: IntArray, noinline initializer: (IntArray) -> T) =
         ndStructure(DefaultStrides(shape), initializer)
 
 
@@ -153,22 +155,22 @@ class MutableBufferNDStructure<T>(
 /**
  * Create optimized mutable structure for given type
  */
-inline fun <reified T: Any> mutableNdStructure(strides: Strides, noinline initializer: (IntArray) -> T) =
-            MutableBufferNDStructure(strides, mutableBuffer(strides.linearSize) { i -> initializer(strides.index(i)) })
+inline fun <reified T : Any> mutableNdStructure(strides: Strides, noinline initializer: (IntArray) -> T) =
+        MutableBufferNDStructure(strides, mutableBuffer(strides.linearSize) { i -> initializer(strides.index(i)) })
 
-inline fun <reified T: Any> mutableNdStructure(shape: IntArray, noinline initializer: (IntArray) -> T) =
+inline fun <reified T : Any> mutableNdStructure(shape: IntArray, noinline initializer: (IntArray) -> T) =
         mutableNdStructure(DefaultStrides(shape), initializer)
 
 /**
  * Create universal mutable structure
  */
-fun <T> genericNdStructure(shape: IntArray, initializer: (IntArray) -> T): MutableBufferNDStructure<T>{
+fun <T> genericNdStructure(shape: IntArray, initializer: (IntArray) -> T): MutableBufferNDStructure<T> {
     val strides = DefaultStrides(shape)
-    val sequence = sequence{
-        strides.indices().forEach{
+    val sequence = sequence {
+        strides.indices().forEach {
             yield(initializer(it))
         }
     }
-    val buffer = MutableListBuffer<T>(sequence.toMutableList())
+    val buffer = MutableListBuffer(sequence.toMutableList())
     return MutableBufferNDStructure(strides, buffer)
 }
diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/BufferSpec.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/BufferSpec.kt
new file mode 100644
index 000000000..a5cc8c5de
--- /dev/null
+++ b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/BufferSpec.kt
@@ -0,0 +1,55 @@
+package scientifik.kmath.structures
+
+import java.nio.ByteBuffer
+
+
+/**
+ * A specification for serialization and deserialization objects to buffer
+ */
+interface BufferSpec<T : Any> {
+    fun fromBuffer(buffer: ByteBuffer): T
+    fun toBuffer(value: T): ByteBuffer
+}
+
+/**
+ * A [BufferSpec] with fixed unit size. Allows storage of any object without boxing.
+ */
+interface FixedSizeBufferSpec<T : Any> : BufferSpec<T> {
+    val unitSize: Int
+
+    /**
+     * Read an object from buffer in current position
+     */
+    fun ByteBuffer.readObject(): T {
+        val buffer = ByteArray(unitSize)
+        get(buffer)
+        return fromBuffer(ByteBuffer.wrap(buffer))
+    }
+
+    /**
+     * Read an object from buffer in given index (not buffer position
+     */
+    fun ByteBuffer.readObject(index: Int): T {
+        val dup = duplicate()
+        dup.position(index*unitSize)
+        return dup.readObject()
+    }
+
+    /**
+     * Write object to [ByteBuffer] in current buffer position
+     */
+    fun ByteBuffer.writeObject(obj: T) {
+        val buffer = toBuffer(obj).apply { rewind() }
+        assert(buffer.limit() == unitSize)
+        put(buffer)
+    }
+
+    /**
+     * Put an object in given index
+     */
+    fun ByteBuffer.writeObject(index: Int, obj: T) {
+        val dup = duplicate()
+        dup.position(index*unitSize)
+        dup.writeObject(obj)
+    }
+}
\ No newline at end of file
diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt
new file mode 100644
index 000000000..68e989e81
--- /dev/null
+++ b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ComplexBufferSpec.kt
@@ -0,0 +1,25 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.Complex
+import java.nio.ByteBuffer
+
+object ComplexBufferSpec : FixedSizeBufferSpec<Complex> {
+    override val unitSize: Int = 16
+
+    override fun fromBuffer(buffer: ByteBuffer): Complex {
+        val re = buffer.getDouble(0)
+        val im = buffer.getDouble(8)
+        return Complex(re, im)
+    }
+
+    override fun toBuffer(value: Complex): ByteBuffer = ByteBuffer.allocate(16).apply {
+        putDouble(value.re)
+        putDouble(value.im)
+    }
+}
+
+/**
+ * Create a mutable buffer which ignores boxing
+ */
+fun Complex.Companion.createBuffer(size: Int) = ObjectBuffer.create(ComplexBufferSpec, size)
+
diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ObjectBuffer.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ObjectBuffer.kt
new file mode 100644
index 000000000..b50ed9674
--- /dev/null
+++ b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/ObjectBuffer.kt
@@ -0,0 +1,28 @@
+package scientifik.kmath.structures
+
+import java.nio.ByteBuffer
+
+class ObjectBuffer<T : Any>(private val buffer: ByteBuffer, private val spec: FixedSizeBufferSpec<T>) : MutableBuffer<T> {
+    override val size: Int
+        get() = buffer.limit() / spec.unitSize
+
+    override fun get(index: Int): T = with(spec) { buffer.readObject(index) }
+
+    override fun iterator(): Iterator<T> = (0 until size).asSequence().map { get(it) }.iterator()
+
+    override fun set(index: Int, value: T) = with(spec) { buffer.writeObject(index, value) }
+
+    override fun copy(): MutableBuffer<T> {
+        val dup = buffer.duplicate()
+        val copy = ByteBuffer.allocate(dup.capacity())
+        dup.rewind()
+        copy.put(dup)
+        copy.flip()
+        return ObjectBuffer(copy, spec)
+    }
+
+    companion object {
+        fun <T : Any> create(spec: FixedSizeBufferSpec<T>, size: Int) =
+                ObjectBuffer<T>(ByteBuffer.allocate(size * spec.unitSize), spec)
+    }
+}
\ No newline at end of file
diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/RealBufferSpec.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/RealBufferSpec.kt
new file mode 100644
index 000000000..761cfc2db
--- /dev/null
+++ b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/RealBufferSpec.kt
@@ -0,0 +1,23 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.Real
+import java.nio.ByteBuffer
+
+object RealBufferSpec : FixedSizeBufferSpec<Real> {
+    override val unitSize: Int = 8
+
+    override fun fromBuffer(buffer: ByteBuffer): Real = Real(buffer.double)
+
+    override fun toBuffer(value: Real): ByteBuffer = ByteBuffer.allocate(8).apply { putDouble(value.value) }
+}
+
+object DoubleBufferSpec : FixedSizeBufferSpec<Double> {
+    override val unitSize: Int = 8
+
+    override fun fromBuffer(buffer: ByteBuffer): Double = buffer.double
+
+    override fun toBuffer(value: Double): ByteBuffer = ByteBuffer.allocate(8).apply { putDouble(value) }
+}
+
+fun Double.Companion.createBuffer(size: Int) = ObjectBuffer.create(DoubleBufferSpec, size)
+fun Real.Companion.createBuffer(size: Int) = ObjectBuffer.create(RealBufferSpec, size)
\ No newline at end of file
diff --git a/kmath-core/src/jvmTest/kotlin/scientifik/kmath/structures/ComplexBufferSpecTest.kt b/kmath-core/src/jvmTest/kotlin/scientifik/kmath/structures/ComplexBufferSpecTest.kt
new file mode 100644
index 000000000..350f06848
--- /dev/null
+++ b/kmath-core/src/jvmTest/kotlin/scientifik/kmath/structures/ComplexBufferSpecTest.kt
@@ -0,0 +1,17 @@
+package scientifik.kmath.structures
+
+import org.junit.Test
+import scientifik.kmath.operations.Complex
+import kotlin.test.assertEquals
+
+class ComplexBufferSpecTest {
+    @Test
+    fun testComplexBuffer() {
+        val buffer = Complex.createBuffer(20)
+        (0 until 20).forEach {
+            buffer[it] = Complex(it.toDouble(), -it.toDouble())
+        }
+
+        assertEquals(Complex(5.0, -5.0), buffer[5])
+    }
+}
\ No newline at end of file
diff --git a/kmath-coroutines/build.gradle b/kmath-coroutines/build.gradle
new file mode 100644
index 000000000..7149e8bb9
--- /dev/null
+++ b/kmath-coroutines/build.gradle
@@ -0,0 +1,42 @@
+plugins {
+    id "org.jetbrains.kotlin.multiplatform"
+}
+
+kotlin {
+    targets {
+        fromPreset(presets.jvm, 'jvm')
+        // For ARM, preset should be changed to presets.iosArm32 or presets.iosArm64
+        // For Linux, preset should be changed to e.g. presets.linuxX64
+        // For MacOS, preset should be changed to e.g. presets.macosX64
+        //fromPreset(presets.mingwX64, 'mingw')
+    }
+    sourceSets {
+        commonMain {
+            dependencies {
+                api project(":kmath-core")
+                api "org.jetbrains.kotlinx:kotlinx-coroutines-core-common:$coroutinesVersion"
+            }
+        }
+        commonTest {
+            dependencies {
+                api 'org.jetbrains.kotlin:kotlin-test-common'
+                api 'org.jetbrains.kotlin:kotlin-test-annotations-common'
+            }
+        }
+        jvmMain {
+            dependencies {
+                api "org.jetbrains.kotlinx:kotlinx-coroutines-core:$coroutinesVersion"
+            }
+        }
+        jvmTest {
+            dependencies {
+                implementation 'org.jetbrains.kotlin:kotlin-test'
+                implementation 'org.jetbrains.kotlin:kotlin-test-junit'
+            }
+        }
+//        mingwMain {
+//        }
+//        mingwTest {
+//        }
+    }
+}
diff --git a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/CoroutinesExtra.kt b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/CoroutinesExtra.kt
new file mode 100644
index 000000000..a837cde01
--- /dev/null
+++ b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/CoroutinesExtra.kt
@@ -0,0 +1,11 @@
+package scientifik.kmath.structures
+
+import kotlinx.coroutines.CoroutineDispatcher
+import kotlinx.coroutines.CoroutineScope
+import kotlinx.coroutines.Dispatchers
+import kotlin.coroutines.CoroutineContext
+import kotlin.coroutines.EmptyCoroutineContext
+
+expect fun <R> runBlocking(context: CoroutineContext = EmptyCoroutineContext, function: suspend CoroutineScope.()->R): R
+
+val Dispatchers.Math: CoroutineDispatcher get() = Dispatchers.Default
\ No newline at end of file
diff --git a/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/LazyNDField.kt b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/LazyNDField.kt
new file mode 100644
index 000000000..823245685
--- /dev/null
+++ b/kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/structures/LazyNDField.kt
@@ -0,0 +1,79 @@
+package scientifik.kmath.structures
+
+import kotlinx.coroutines.*
+import scientifik.kmath.operations.Field
+
+class LazyNDField<T, F : Field<T>>(shape: IntArray, field: F, val scope: CoroutineScope = GlobalScope) : NDField<T, F>(shape, field) {
+
+    override fun produceStructure(initializer: F.(IntArray) -> T): NDStructure<T> = LazyNDStructure(this) { initializer(field, it) }
+
+
+    override fun add(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
+        return LazyNDStructure(this) { index ->
+            val aDeferred = a.deferred(index)
+            val bDeferred = b.deferred(index)
+            aDeferred.await() + bDeferred.await()
+        }
+    }
+
+    override fun multiply(a: NDElement<T, F>, k: Double): NDElement<T, F> {
+        return LazyNDStructure(this) { index -> a.await(index) * k }
+    }
+
+    override fun multiply(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
+        return LazyNDStructure(this) { index ->
+            val aDeferred = a.deferred(index)
+            val bDeferred = b.deferred(index)
+            aDeferred.await() * bDeferred.await()
+        }
+    }
+
+    override fun divide(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
+        return LazyNDStructure(this) { index ->
+            val aDeferred = a.deferred(index)
+            val bDeferred = b.deferred(index)
+            aDeferred.await() / bDeferred.await()
+        }
+    }
+}
+
+class LazyNDStructure<T, F : Field<T>>(override val context: LazyNDField<T, F>, val function: suspend F.(IntArray) -> T) : NDElement<T, F>, NDStructure<T> {
+    override val self: NDElement<T, F> get() = this
+    override val shape: IntArray get() = context.shape
+
+    private val cache = HashMap<IntArray, Deferred<T>>()
+
+    fun deferred(index: IntArray) = cache.getOrPut(index) { context.scope.async(context = Dispatchers.Math) { function.invoke(context.field, index) } }
+
+    suspend fun await(index: IntArray): T = deferred(index).await()
+
+    override fun get(index: IntArray): T = runBlocking {
+        deferred(index).await()
+    }
+
+    override fun elements(): Sequence<Pair<IntArray, T>> {
+        val strides = DefaultStrides(shape)
+        return strides.indices().map { index -> index to runBlocking { await(index) } }
+    }
+}
+
+fun <T> NDElement<T, *>.deferred(index: IntArray) = if (this is LazyNDStructure<T, *>) this.deferred(index) else CompletableDeferred(get(index))
+
+suspend fun <T> NDElement<T, *>.await(index: IntArray) = if (this is LazyNDStructure<T, *>) this.await(index) else get(index)
+
+fun <T, F : Field<T>> NDElement<T, F>.lazy(scope: CoroutineScope = GlobalScope): LazyNDStructure<T, F> {
+    return if (this is LazyNDStructure<T, F>) {
+        this
+    } else {
+        val context = LazyNDField(context.shape, context.field)
+        LazyNDStructure(context) { get(it) }
+    }
+}
+
+inline fun <T, F : Field<T>> LazyNDStructure<T, F>.transformIndexed(crossinline action: suspend F.(IntArray, T) -> T) = LazyNDStructure(context) { index ->
+    action.invoke(this, index, await(index))
+}
+
+inline fun <T, F : Field<T>> LazyNDStructure<T, F>.transform(crossinline action: suspend F.(T) -> T) = LazyNDStructure(context) { index ->
+    action.invoke(this, await(index))
+}
\ No newline at end of file
diff --git a/kmath-coroutines/src/commonTest/kotlin/scientifik/kmath/structures/LazyNDFieldTest.kt b/kmath-coroutines/src/commonTest/kotlin/scientifik/kmath/structures/LazyNDFieldTest.kt
new file mode 100644
index 000000000..403fe2d31
--- /dev/null
+++ b/kmath-coroutines/src/commonTest/kotlin/scientifik/kmath/structures/LazyNDFieldTest.kt
@@ -0,0 +1,20 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.IntField
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+
+class LazyNDFieldTest {
+    @Test
+    fun testLazyStructure() {
+        var counter = 0
+        val regularStructure = NDArrays.create(IntField, intArrayOf(2, 2, 2)) { it[0] + it[1] - it[2] }
+        val result = (regularStructure.lazy() + 2).transform {
+            counter++
+            it * it
+        }
+        assertEquals(4, result[0,0,0])
+        assertEquals(1, counter)
+    }
+}
\ No newline at end of file
diff --git a/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/_CoroutinesExtra.kt b/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/_CoroutinesExtra.kt
new file mode 100644
index 000000000..72a764729
--- /dev/null
+++ b/kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/_CoroutinesExtra.kt
@@ -0,0 +1,6 @@
+package scientifik.kmath.structures
+
+import kotlinx.coroutines.CoroutineScope
+import kotlin.coroutines.CoroutineContext
+
+actual fun <R> runBlocking(context: CoroutineContext, function: suspend CoroutineScope.() -> R): R  = kotlinx.coroutines.runBlocking(context, function)
\ No newline at end of file
diff --git a/kmath-io/build.gradle b/kmath-io/build.gradle
index d0d93262e..28fb7eee5 100644
--- a/kmath-io/build.gradle
+++ b/kmath-io/build.gradle
@@ -1,5 +1,5 @@
 plugins {
-    id 'kotlin-multiplatform'
+    id "org.jetbrains.kotlin.multiplatform"
 }
 
 kotlin {
@@ -16,6 +16,7 @@ kotlin {
             dependencies {
                 api project(":kmath-core")
                 implementation 'org.jetbrains.kotlin:kotlin-stdlib-common'
+                api "org.jetbrains.kotlinx:kotlinx-io:$ioVersion"
             }
         }
         commonTest {
@@ -27,6 +28,7 @@ kotlin {
         jvmMain {
             dependencies {
                 implementation 'org.jetbrains.kotlin:kotlin-stdlib-jdk8'
+                api "org.jetbrains.kotlinx:kotlinx-io-jvm:$ioVersion"
             }
         }
         jvmTest {
diff --git a/kmath-jmh/build.gradle b/kmath-jmh/build.gradle
index 97337cb40..1876d86d7 100644
--- a/kmath-jmh/build.gradle
+++ b/kmath-jmh/build.gradle
@@ -5,6 +5,6 @@ plugins {
 }
 
 dependencies {
-    implementation project(':kmath-core')
-    jmh 'org.jetbrains.kotlin:kotlin-stdlib-jdk8'
+    compile project(':kmath-core')
+    //jmh project(':kmath-core')
 }
\ No newline at end of file
diff --git a/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt b/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt
index eb40195e0..fe10fbd75 100644
--- a/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt
+++ b/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt
@@ -4,8 +4,7 @@ import org.openjdk.jmh.annotations.*
 import java.nio.IntBuffer
 
 
-@Fork(1)
-@Warmup(iterations = 2)
+@Warmup(iterations = 1)
 @Measurement(iterations = 5)
 @State(Scope.Benchmark)
 open class ArrayBenchmark {
@@ -30,7 +29,6 @@ open class ArrayBenchmark {
         for (i in 1..10000) {
             res += array[10000 - i]
         }
-        print(res)
     }
 
     @Benchmark
@@ -39,7 +37,6 @@ open class ArrayBenchmark {
         for (i in 1..10000) {
             res += arrayBuffer.get(10000 - i)
         }
-        print(res)
     }
 
     @Benchmark
@@ -48,6 +45,5 @@ open class ArrayBenchmark {
         for (i in 1..10000) {
             res += nativeBuffer.get(10000 - i)
         }
-        print(res)
     }
 }
\ No newline at end of file
diff --git a/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/BufferBenchmark.kt b/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/BufferBenchmark.kt
new file mode 100644
index 000000000..fd27ae3e0
--- /dev/null
+++ b/kmath-jmh/src/jmh/kotlin/scientifik/kmath/structures/BufferBenchmark.kt
@@ -0,0 +1,38 @@
+package scientifik.kmath.structures
+
+import org.openjdk.jmh.annotations.*
+import scientifik.kmath.operations.Complex
+
+@Warmup(iterations = 1)
+@Measurement(iterations = 5)
+@State(Scope.Benchmark)
+open class BufferBenchmark {
+
+    @Benchmark
+    fun genericDoubleBufferReadWrite() {
+        val buffer = Double.createBuffer(size)
+        (0 until size).forEach {
+            buffer[it] = it.toDouble()
+        }
+
+        (0 until size).forEach {
+            buffer[it]
+        }
+    }
+
+    @Benchmark
+    fun complexBufferReadWrite() {
+        val buffer = Complex.createBuffer(size/2)
+        (0 until size/2).forEach {
+            buffer[it] = Complex(it.toDouble(), -it.toDouble())
+        }
+
+        (0 until size/2).forEach {
+            buffer[it]
+        }
+    }
+
+    companion object {
+        const val size = 1000
+    }
+}
\ No newline at end of file
diff --git a/settings.gradle.kts b/settings.gradle.kts
index 0f90d4c0b..0c91b8446 100644
--- a/settings.gradle.kts
+++ b/settings.gradle.kts
@@ -5,10 +5,12 @@ pluginManagement {
     }
 }
 
-enableFeaturePreview("GRADLE_METADATA")
+//enableFeaturePreview("GRADLE_METADATA")
 
 rootProject.name = "kmath"
 include(
         ":kmath-core",
+        ":kmath-io",
+        ":kmath-coroutines",
         ":kmath-jmh"
 )