Matrix and linear algebra redone

build.gradle.kts

@@ -1,13 +1,14 @@
 buildscript {
-    extra["kotlinVersion"] = "1.3.11"
+    extra["kotlinVersion"] = "1.3.20-eap-52"
     extra["ioVersion"] = "0.1.2-dev-2"
-    extra["coroutinesVersion"] = "1.0.1"
+    extra["coroutinesVersion"] = "1.1.0"
 
     val kotlinVersion: String by extra
     val ioVersion: String by extra
     val coroutinesVersion: String by extra
 
     repositories {
+        maven ("https://dl.bintray.com/kotlin/kotlin-eap")
         jcenter()
     }
 
@@ -28,6 +29,11 @@ allprojects {
 
     group = "scientifik"
     version = "0.0.2-dev-1"
+
+    repositories{
+        maven ("https://dl.bintray.com/kotlin/kotlin-eap")
+        jcenter()
+    }
 }
 
 if(file("artifactory.gradle").exists()){

kmath-core/src/commonMain/kotlin/scientifik/kmath/histogram/FastHistogram.kt

@@ -20,7 +20,7 @@ class FastHistogram(
 
     private val strides = DefaultStrides(IntArray(binNums.size) { binNums[it] + 2 })
 
-    private val values: NDStructure<LongCounter> = ndStructure(strides) { LongCounter() }
+    private val values: NDStructure<LongCounter> = inlineNDStructure(strides) { LongCounter() }
 
     //private val weight: NDStructure<DoubleCounter?> = ndStructure(strides){null}
 
@@ -93,7 +93,7 @@ class FastHistogram(
      * Convert this histogram into NDStructure containing bin values but not bin descriptions
      */
     fun asNDStructure(): NDStructure<Number> {
-        return ndStructure(this.values.shape) { values[it].sum() }
+        return inlineNdStructure(this.values.shape) { values[it].sum() }
     }
 
     /**

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt

@@ -2,10 +2,7 @@ package scientifik.kmath.linear
 
 import scientifik.kmath.operations.DoubleField
 import scientifik.kmath.operations.Field
-import scientifik.kmath.structures.MutableNDStructure
+import scientifik.kmath.structures.*
-import scientifik.kmath.structures.NDStructure
-import scientifik.kmath.structures.get
-import scientifik.kmath.structures.mutableNdStructure
 import kotlin.math.absoluteValue
 
 /**
@@ -13,7 +10,7 @@ import kotlin.math.absoluteValue
  */
 abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matrix<T, F>) {
 
-    private val field get() = matrix.context.field
+    private val field get() = matrix.context.ring
     /** Entries of LU decomposition.  */
     internal val lu: NDStructure<T>
     /** Pivot permutation associated with LU decomposition.  */
@@ -34,11 +31,11 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
      * @return the L matrix (or null if decomposed matrix is singular)
      */
     val l: Matrix<out T, F> by lazy {
-        matrix.context.produce { i, j ->
+        matrix.context.produce(matrix.numRows, matrix.numCols) { i, j ->
             when {
                 j < i -> lu[i, j]
-                j == i -> matrix.context.field.one
+                j == i -> matrix.context.ring.one
-                else -> matrix.context.field.zero
+                else -> matrix.context.ring.zero
             }
         }
     }
@@ -51,7 +48,7 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
      * @return the U matrix (or null if decomposed matrix is singular)
      */
     val u: Matrix<out T, F> by lazy {
-        matrix.context.produce { i, j ->
+        matrix.context.produce(matrix.numRows, matrix.numCols) { i, j ->
             if (j >= i) lu[i, j] else field.zero
         }
     }
@@ -67,7 +64,7 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
      * @see .getPivot
      */
     val p: Matrix<out T, F> by lazy {
-        matrix.context.produce { i, j ->
+        matrix.context.produce(matrix.numRows, matrix.numCols) { i, j ->
             //TODO ineffective. Need sparse matrix for that
             if (j == pivot[i]) field.one else field.zero
         }
@@ -79,9 +76,9 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
      */
     val determinant: T
         get() {
-            with(matrix.context.field) {
+            with(matrix.context.ring) {
                 var determinant = if (even) one else -one
-                for (i in 0 until matrix.rows) {
+                for (i in 0 until matrix.numRows) {
                     determinant *= lu[i, i]
                 }
                 return determinant
@@ -97,20 +94,20 @@ abstract class LUDecomposition<T : Comparable<T>, F : Field<T>>(val matrix: Matr
 
     abstract fun isSingular(value: T): Boolean
 
-    private fun abs(value: T) = if (value > matrix.context.field.zero) value else with(matrix.context.field) { -value }
+    private fun abs(value: T) = if (value > matrix.context.ring.zero) value else with(matrix.context.ring) { -value }
 
     private fun calculateLU(): Pair<NDStructure<T>, IntArray> {
-        if (matrix.rows != matrix.columns) {
+        if (matrix.numRows != matrix.numCols) {
             error("LU decomposition supports only square matrices")
         }
 
-        val m = matrix.columns
+        val m = matrix.numCols
-        val pivot = IntArray(matrix.rows)
+        val pivot = IntArray(matrix.numRows)
         //TODO fix performance
-        val lu: MutableNDStructure<T> = mutableNdStructure(intArrayOf(matrix.rows, matrix.columns)) { index -> matrix[index[0], index[1]] }
+        val lu: MutableNDStructure<T> = MutableNdStructure(intArrayOf(matrix.numRows, matrix.numCols), ::boxingMutableBuffer) { index: IntArray -> matrix[index[0], index[1]] }
 
 
-        with(matrix.context.field) {
+        with(matrix.context.ring) {
             // Initialize permutation array and parity
             for (row in 0 until m) {
                 pivot[row] = row
@@ -201,50 +198,50 @@ object RealLUSolver : LinearSolver<Double, DoubleField> {
     fun decompose(mat: Matrix<Double, DoubleField>, threshold: Double = 1e-11): RealLUDecomposition = RealLUDecomposition(mat, threshold)
 
     override fun solve(a: RealMatrix, b: RealMatrix): RealMatrix {
-        val decomposition = decompose(a, a.context.field.zero)
+        val decomposition = decompose(a, a.context.ring.zero)
 
-        if (b.rows != a.rows) {
+        if (b.numRows != a.numCols) {
-            error("Matrix dimension mismatch expected ${a.rows}, but got ${b.rows}")
+            error("Matrix dimension mismatch expected ${a.numRows}, but got ${b.numCols}")
         }
 
         // Apply permutations to b
-        val bp = Array(a.rows) { DoubleArray(b.columns) }
+        val bp = Array(a.numRows) { DoubleArray(b.numCols) }
-        for (row in 0 until a.rows) {
+        for (row in 0 until a.numRows) {
             val bpRow = bp[row]
             val pRow = decomposition.pivot[row]
-            for (col in 0 until b.columns) {
+            for (col in 0 until b.numCols) {
                 bpRow[col] = b[pRow, col]
             }
         }
 
         // Solve LY = b
-        for (col in 0 until a.rows) {
+        for (col in 0 until a.numRows) {
             val bpCol = bp[col]
-            for (i in col + 1 until a.rows) {
+            for (i in col + 1 until a.numRows) {
                 val bpI = bp[i]
                 val luICol = decomposition.lu[i, col]
-                for (j in 0 until b.columns) {
+                for (j in 0 until b.numCols) {
                     bpI[j] -= bpCol[j] * luICol
                 }
             }
         }
 
         // Solve UX = Y
-        for (col in a.rows - 1 downTo 0) {
+        for (col in a.numRows - 1 downTo 0) {
             val bpCol = bp[col]
             val luDiag = decomposition.lu[col, col]
-            for (j in 0 until b.columns) {
+            for (j in 0 until b.numCols) {
                 bpCol[j] /= luDiag
             }
             for (i in 0 until col) {
                 val bpI = bp[i]
                 val luICol = decomposition.lu[i, col]
-                for (j in 0 until b.columns) {
+                for (j in 0 until b.numCols) {
                     bpI[j] -= bpCol[j] * luICol
                 }
             }
         }
 
-        return a.context.produce { i, j -> bp[i][j] }
+        return a.context.produce(a.numRows, a.numCols) { i, j -> bp[i][j] }
     }
 }

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt

@@ -3,6 +3,9 @@ package scientifik.kmath.linear
 import scientifik.kmath.operations.DoubleField
 import scientifik.kmath.operations.Field
 import scientifik.kmath.operations.Norm
+import scientifik.kmath.operations.Ring
+import scientifik.kmath.structures.asSequence
+import scientifik.kmath.structures.boxingBuffer
 
 
 /**
@@ -11,23 +14,22 @@ import scientifik.kmath.operations.Norm
 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))
+    fun inverse(a: Matrix<T, F>): Matrix<T, F> = solve(a, a.context.one)
 }
 
 /**
  * Convert vector to array (copying content of array)
  */
-fun <T : Any> Array<T>.toVector(field: Field<T>) = Vector.of(size, field) { this[it] }
+fun <T : Any> Array<T>.toVector(field: Field<T>) = Vector.generic(size, field) { this[it] }
 
-fun DoubleArray.toVector() = Vector.ofReal(this.size) { this[it] }
+fun DoubleArray.toVector() = Vector.real(this.size) { this[it] }
-fun List<Double>.toVector() = Vector.ofReal(this.size) { this[it] }
+fun List<Double>.toVector() = Vector.real(this.size) { this[it] }
 
 /**
  * Convert matrix to vector if it is possible
  */
-fun <T : Any, F : Field<T>> Matrix<T, F>.toVector(): Vector<T, F> {
+fun <T : Any, F : Ring<T>> Matrix<T, F>.toVector(): Vector<T, F> {
-    return when {
+    return if (this.numCols == 1) {
-        this.columns == 1 -> {
 //            if (this is ArrayMatrix) {
 //                //Reuse existing underlying array
 //                ArrayVector(ArrayVectorSpace(rows, context.field, context.ndFactory), array)
@@ -35,26 +37,27 @@ fun <T : Any, F : Field<T>> Matrix<T, F>.toVector(): Vector<T, F> {
 //                //Generic vector
 //                vector(rows, context.field) { get(it, 0) }
 //            }
-            Vector.of(rows, context.field) { get(it, 0) }
+        Vector.generic(numRows, context.ring) { get(it, 0) }
-        }
+    } else error("Can't convert matrix with more than one column to vector")
-        else -> error("Can't convert matrix with more than one column to vector")
-    }
 }
 
-fun <T : Any, F : Field<T>> Vector<T, F>.toMatrix(): Matrix<T, F> {
+fun <T : Any, R : Ring<T>> Vector<T, R>.toMatrix(): Matrix<T, R> {
+//    val context = StructureMatrixContext(size, 1, context.space)
+//
 //    return if (this is ArrayVector) {
 //        //Reuse existing underlying array
-//        ArrayMatrix(ArrayMatrixSpace(size, 1, context.field, context.ndFactory), array)
+//        StructureMatrix(context,this.buffer)
 //    } else {
 //        //Generic vector
 //        matrix(size, 1, context.field) { i, j -> get(i) }
 //    }
-    return Matrix.of(size, 1, context.space) { i, _ -> get(i) }
+    //return Matrix.of(size, 1, context.space) { i, _ -> get(i) }
+    return StructureMatrixSpace(size, 1, context.space, ::boxingBuffer).produce { 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() })
+        return kotlin.math.sqrt(arg.asSequence().sumByDouble { it.toDouble() })
     }
 }
 

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Matrix.kt

@@ -1,198 +1,156 @@
 package scientifik.kmath.linear
 
 import scientifik.kmath.histogram.Point
-import scientifik.kmath.operations.*
+import scientifik.kmath.operations.DoubleField
+import scientifik.kmath.operations.Ring
+import scientifik.kmath.operations.Space
+import scientifik.kmath.operations.SpaceElement
 import scientifik.kmath.structures.*
 
-/**
+
- * The space for linear elements. Supports scalar product alongside with standard linear operations.
+interface MatrixSpace<T : Any, R : Ring<T>> : Space<Matrix<T, R>> {
- * @param T type of individual element of the vector or matrix
+    /**
- * @param V the type of vector space element
+     * The ring context for matrix elements
- */
+     */
-abstract class MatrixSpace<T : Any, F : Ring<T>>(val rows: Int, val columns: Int, val field: F) : Space<Matrix<T, F>> {
+    val ring: R
+
+    val rowNum: Int
+    val colNum: Int
+
+    val shape get() = intArrayOf(rowNum, colNum)
 
     /**
-     * Produce the element of this space
+     * Produce a matrix with this context and given dimensions
      */
-    abstract fun produce(initializer: (Int, Int) -> T): Matrix<T, F>
+    fun produce(rows: Int = rowNum, columns: Int = colNum, initializer: (i: Int, j: Int) -> T): Matrix<T, R>
 
     /**
-     * Produce new matrix space with given dimensions. The space produced could be raised from cache since [MatrixSpace] does not have mutable elements
+     * Produce a point compatible with matrix space
      */
-    abstract fun produceSpace(rows: Int, columns: Int): MatrixSpace<T, F>
+    fun point(size: Int, initializer: (Int) -> T): Point<T>
 
-    override val zero: Matrix<T, F> by lazy {
+    override val zero: Matrix<T, R> get() = produce { _, _ -> ring.zero }
-        produce { _, _ -> field.zero }
-    }
 
-//    val one: Matrix<T> by lazy {
+    val one get() = produce { i, j -> if (i == j) ring.one else ring.zero }
-//        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> {
+    override fun add(a: Matrix<T, R>, b: Matrix<T, R>): Matrix<T, R> = produce(rowNum, colNum) { i, j -> ring.run { a[i, j] + b[i, j] } }
-        return produce { i, j -> with(field) { a[i, j] + b[i, j] } }
-    }
 
-    override fun multiply(a: Matrix<T, F>, k: Double): Matrix<T, F> {
+    override fun multiply(a: Matrix<T, R>, k: Double): Matrix<T, R> = produce(rowNum, colNum) { i, j -> ring.run { a[i, j] * k } }
-        //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 } }
-    }
-
-    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
-    }
-}
-
-/**
- * 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 {
+        /**
+         * Non-boxing double matrix
+         */
+        fun real(rows: Int, columns: Int): MatrixSpace<Double, DoubleField> = StructureMatrixSpace(rows, columns, DoubleField, DoubleBufferFactory)
 
         /**
-         * Create [ArrayMatrix] with custom field
+         * A structured matrix with custom buffer
          */
-        fun <T : Any, F: Field<T>> of(rows: Int, columns: Int, field: F, initializer: (Int, Int) -> T) =
+        fun <T : Any, R : Ring<T>> buffered(rows: Int, columns: Int, ring: R, bufferFactory: BufferFactory<T> = ::boxingBuffer): MatrixSpace<T, R> = StructureMatrixSpace(rows, columns, ring, bufferFactory)
-                ArrayMatrix(ArrayMatrixSpace(rows, columns, field), initializer)
 
         /**
-         * Create [ArrayMatrix] of doubles. The implementation in general should be faster than generic one due to boxing.
+         * Automatic buffered matrix, unboxed if it is possible
          */
-        fun ofReal(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
+        inline fun <reified T : Any, R : Ring<T>> smart(rows: Int, columns: Int, ring: R): MatrixSpace<T, R> = buffered(rows, columns, ring, ::inlineBuffer)
-                ArrayMatrix(ArrayMatrixSpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
+    }
+}
 
-        /**
+
-         * Create a diagonal value matrix. By default value equals [Field.one].
+/**
-         */
+ * Specialized 2-d structure
-        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 }
+interface Matrix<T : Any, R : Ring<T>> : NDStructure<T>, SpaceElement<Matrix<T, R>, MatrixSpace<T, R>> {
+    operator fun get(i: Int, j: Int): T
+
+    override fun get(index: IntArray): T = get(index[0], index[1])
+
+    override val shape: IntArray get() = context.shape
+
+    val numRows get() = context.rowNum
+    val numCols get() = context.colNum
+
+    //TODO replace by lazy buffers
+    val rows: List<Point<T>>
+        get() = (0 until numRows).map { i ->
+            context.point(numCols) { j -> get(i, j) }
         }
 
-        /**
+    val columns: List<Point<T>>
-         * Equality check on two generic matrices
+        get() = (0 until numCols).map { j ->
-         */
+            context.point(numRows) { i -> get(i, j) }
-        fun equals(mat1: Matrix<*, *>, mat2: Matrix<*, *>): Boolean {
+        }
-            if (mat1 === mat2) return true
+
-            if (mat1.context != mat2.context) return false
+    companion object {
-            for (i in 0 until mat1.rows) {
+        fun real(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
-                for (j in 0 until mat2.columns) {
+                MatrixSpace.real(rows, columns).produce(rows, columns, initializer)
-                    if (mat1[i, j] != mat2[i, j]) return false
+    }
-                }
+}
-            }
+
-            return true
+infix fun <T : Any, R : Ring<T>> Matrix<T, R>.dot(other: Matrix<T, R>): Matrix<T, R> {
+    //TODO add typed error
+    if (this.numCols != other.numRows) error("Matrix dot operation dimension mismatch: ($numRows, $numCols) x (${other.numRows}, ${other.numCols})")
+    return context.produce(numRows, other.numCols) { i, j ->
+        val row = rows[i]
+        val column = other.columns[j]
+        with(context.ring) {
+            row.asSequence().zip(column.asSequence(), ::multiply).sum()
         }
     }
 }
 
-/**
+infix fun <T : Any, R : Ring<T>> Matrix<T, R>.dot(vector: Point<T>): Point<T> {
- * Dot product. Throws exception on dimension mismatch
+    //TODO add typed error
- */
+    if (this.numCols != vector.size) error("Matrix dot vector operation dimension mismatch: ($numRows, $numCols) x (${vector.size})")
-infix fun <T : Any, F : Field<T>> Matrix<T, F>.dot(b: Matrix<T, F>): Matrix<T, F> {
+    return context.point(numRows) { i ->
-    if (columns != b.rows) {
+        val row = rows[i]
-        //TODO replace by specific exception
+        with(context.ring) {
-        error("Dimension mismatch in linear structure dot product: [$rows,$columns]*[${b.rows},${b.columns}]")
+            row.asSequence().zip(vector.asSequence(), ::multiply).sum()
-    }
+        }
-    return context.produceSpace(rows, b.columns).produce { i, j ->
-        (0 until columns).asSequence().map { k -> context.field.multiply(this[i, k], b[k, j]) }.reduce { first, second -> context.field.add(first, second) }
     }
 }
 
-/**
+data class StructureMatrixSpace<T : Any, R : Ring<T>>(
- * Matrix x Vector dot product.
+        override val rowNum: Int,
- */
+        override val colNum: Int,
-infix fun <T : Any, F : Field<T>> Matrix<T, F>.dot(b: Point<T>): Matrix<T, F> {
+        override val ring: R,
-    if (columns != b.size) {
+        private val bufferFactory: BufferFactory<T>
-        //TODO replace by specific exception
+) : MatrixSpace<T, R> {
-        error("Dimension mismatch in linear structure dot product: [$rows,$columns]*[${b.size},1]")
+
-    }
+    override val shape: IntArray = intArrayOf(rowNum, colNum)
-    return context.produceSpace(rows, 1).produce { i, j ->
+
-        (0 until columns).asSequence().map { k -> context.field.multiply(this[i, k], b[k]) }.reduce { first, second -> context.field.add(first, second) }
+    private val strides = DefaultStrides(shape)
+
+    override fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> T): Matrix<T, R> {
+        return if (rows == rowNum && columns == colNum) {
+            val structure = NdStructure(strides, bufferFactory) { initializer(it[0], it[1]) }
+            StructureMatrix(this, structure)
+        } else {
+            val context = StructureMatrixSpace(rows, columns, ring, bufferFactory)
+            val structure = NdStructure(context.strides, bufferFactory) { initializer(it[0], it[1]) }
+            StructureMatrix(context, structure)
+        }
     }
+
+    override fun point(size: Int, initializer: (Int) -> T): Point<T> = bufferFactory(size, initializer)
 }
 
-
+data class StructureMatrix<T : Any, R : Ring<T>>(override val context: StructureMatrixSpace<T, R>, val structure: NDStructure<T>) : Matrix<T, R> {
-
+    init {
-typealias NDFieldFactory<T, F> = (IntArray) -> NDField<T, F>
+        if (structure.shape.size != 2 || structure.shape[0] != context.rowNum || structure.shape[1] != context.colNum) {
-
+            error("Dimension mismatch for structure, (${context.rowNum}, ${context.colNum}) expected, but ${structure.shape} found")
-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 val shape: IntArray get() = structure.shape
+    override val self: Matrix<T, R> get() = this
 
-    override fun produceSpace(rows: Int, columns: Int): ArrayMatrixSpace<T, F> {
+    override fun get(index: IntArray): T = structure[index]
-        return ArrayMatrixSpace(rows, columns, field, ndFactory)
+
-    }
+    override fun get(i: Int, j: Int): T = structure[i, j]
+
+    override fun elements(): Sequence<Pair<IntArray, T>> = structure.elements()
 }
 
-/**
+//TODO produce transposed matrix via reference without creating new space and structure
- * Member of [ArrayMatrixSpace] which wraps 2-D array
+fun <T : Any, R : Ring<T>> Matrix<T, R>.transpose(): Matrix<T, R> =
- */
+        context.produce(numCols, numRows) { i, j -> get(j, i) }
-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
-}

kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/Vector.kt

@@ -2,32 +2,58 @@ 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.*
-import scientifik.kmath.structures.get
 
 /**
  * A linear space for vectors.
  * Could be used on any point-like structure
  */
-abstract class VectorSpace<T : Any, S : Space<T>>(val size: Int, val space: S) : Space<Point<T>> {
+interface VectorSpace<T : Any, S : Space<T>> : Space<Point<T>> {
 
-    abstract fun produce(initializer: (Int) -> T): Vector<T, S>
+    val size: Int
 
-    override val zero: Vector<T, S> by lazy { produce { space.zero } }
+    val space: S
+
+    fun produce(initializer: (Int) -> T): Vector<T, S>
+
+    override val zero: Vector<T, S> get() = produce { space.zero }
 
     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, S> = produce { with(space) { a[it] * k } }
+
+    //TODO add basis
+
+    companion object {
+
+        private val realSpaceCache = HashMap<Int, BufferVectorSpace<Double, DoubleField>>()
+
+        /**
+         * Non-boxing double vector space
+         */
+        fun real(size: Int): BufferVectorSpace<Double, DoubleField> {
+            return realSpaceCache.getOrPut(size) { BufferVectorSpace(size, DoubleField, DoubleBufferFactory) }
+        }
+
+        /**
+         * A structured vector space with custom buffer
+         */
+        fun <T : Any, S : Space<T>> buffered(size: Int, space: S, bufferFactory: BufferFactory<T> = ::boxingBuffer): VectorSpace<T, S> = BufferVectorSpace(size, space, bufferFactory)
+
+        /**
+         * Automatic buffered vector, unboxed if it is possible
+         */
+        inline fun <reified T : Any, S : Space<T>> smart(size: Int, space: S): VectorSpace<T, S> = buffered(size, space, ::inlineBuffer)
+    }
 }
 
 
 /**
  * A point coupled to the linear space
  */
-interface Vector<T : Any, S : Space<T>> : SpaceElement<Point<T>, VectorSpace<T, S>>, Point<T>, Iterable<T> {
+interface Vector<T : Any, S : Space<T>> : SpaceElement<Point<T>, VectorSpace<T, S>>, Point<T> {
     override val size: Int get() = context.size
 
     override operator fun plus(b: Point<T>): Vector<T, S> = context.add(self, b)
@@ -39,65 +65,40 @@ 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 : Space<T>> generic(size: Int, field: F, initializer: (Int) -> T) =
-                ArrayVector(ArrayVectorSpace(size, field), initializer)
+                VectorSpace.buffered(size, field).produce(initializer)
 
-        private val realSpaceCache = HashMap<Int, ArrayVectorSpace<Double, DoubleField>>()
+        fun real(size: Int, initializer: (Int) -> Double) = VectorSpace.real(size).produce(initializer)
+        fun ofReal(vararg elements: Double) = VectorSpace.real(elements.size).produce{elements[it]}
 
-        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(getRealSpace(size), 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, F : Field<T>>(
+data class BufferVectorSpace<T : Any, S : Space<T>>(
-        size: Int,
+        override val size: Int,
-        field: F,
+        override val space: S,
-        val ndFactory: NDFieldFactory<T, F> = genericNDFieldFactory(field)
+        val bufferFactory: BufferFactory<T>
-) : VectorSpace<T, F>(size, field) {
+) : VectorSpace<T, S> {
-    val ndField by lazy {
+    override fun produce(initializer: (Int) -> T): Vector<T, S> = BufferVector(this, bufferFactory(size, initializer))
-        ndFactory(intArrayOf(size))
-    }
-
-    override fun produce(initializer: (Int) -> T): Vector<T, F> = ArrayVector(this, initializer)
 }
 
 
-class ArrayVector<T : Any, F : Field<T>> internal constructor(override val context: VectorSpace<T, F>, val element: NDElement<T, F>) : Vector<T, F> {
+data class BufferVector<T : Any, S : Space<T>>(override val context: VectorSpace<T, S>, val buffer: Buffer<T>) : Vector<T, S> {
-
-    constructor(context: ArrayVectorSpace<T, F>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
 
     init {
-        if (context.size != element.shape[0]) {
+        if (context.size != buffer.size) {
             error("Array dimension mismatch")
         }
     }
 
     override fun get(index: Int): T {
-        return element[index]
+        return buffer[index]
     }
 
-    override val self: ArrayVector<T, F> get() = this
+    override val self: BufferVector<T, S> get() = this
 
-    override fun iterator(): Iterator<T> = (0 until size).map { element[it] }.iterator()
+    override fun iterator(): Iterator<T> = (0 until size).map { buffer[it] }.iterator()
 
-    override fun toString(): String = this.joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() }
+    override fun toString(): String = this.asSequence().joinToString(prefix = "[", postfix = "]", separator = ", ") { it.toString() }
 }
 

kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt

@@ -51,6 +51,9 @@ interface Space<T> {
     operator fun T.div(k: Number) = multiply(this, 1.0 / k.toDouble())
     operator fun Number.times(b: T) = b * this
 
+    //TODO move to external extensions when they are available
+    fun Iterable<T>.sum(): T = fold(zero) { left, right -> left + right }
+    fun Sequence<T>.sum(): T = fold(zero) { left, right -> left + right }
 }
 
 /**
@@ -111,8 +114,8 @@ interface Field<T> : Ring<T> {
 /**
  * Field element
  */
-interface FieldElement<T, S : Field<T>> : RingElement<T, S> {
+interface FieldElement<T, F : Field<T>> : RingElement<T, F> {
-    override val context: S
+    override val context: F
 
     operator fun div(b: T): T = context.divide(self, b)
 }

kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt

@@ -5,11 +5,11 @@ import kotlin.math.pow
 /**
  * Advanced Number-like field that implements basic operations
  */
-interface ExtendedField<N : Any> :
+interface ExtendedField<T : Any> :
-        Field<N>,
+        Field<T>,
-        TrigonometricOperations<N>,
+        TrigonometricOperations<T>,
-        PowerOperations<N>,
+        PowerOperations<T>,
-        ExponentialOperations<N>
+        ExponentialOperations<T>
 
 
 /**

kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/Buffers.kt

@@ -11,10 +11,14 @@ interface Buffer<T> {
     operator fun get(index: Int): T
 
     operator fun iterator(): Iterator<T>
+
+    fun contentEquals(other: Buffer<*>): Boolean = asSequence().mapIndexed { index, value -> value == other[index] }.all { it }
 }
 
 fun <T> Buffer<T>.asSequence(): Sequence<T> = iterator().asSequence()
 
+fun <T> Buffer<T>.asIterable(): Iterable<T> = iterator().asSequence().asIterable()
+
 interface MutableBuffer<T> : Buffer<T> {
     operator fun set(index: Int, value: T)
 
@@ -95,6 +99,20 @@ inline class IntBuffer(private val array: IntArray) : MutableBuffer<Int> {
     override fun copy(): MutableBuffer<Int> = IntBuffer(array.copyOf())
 }
 
+inline class LongBuffer(private val array: LongArray) : MutableBuffer<Long> {
+    override val size: Int get() = array.size
+
+    override fun get(index: Int): Long = array[index]
+
+    override fun set(index: Int, value: Long) {
+        array[index] = value
+    }
+
+    override fun iterator(): Iterator<Long> = array.iterator()
+
+    override fun copy(): MutableBuffer<Long> = LongBuffer(array.copyOf())
+}
+
 inline class ReadOnlyBuffer<T>(private val buffer: MutableBuffer<T>) : Buffer<T> {
     override val size: Int get() = buffer.size
 
@@ -112,26 +130,46 @@ fun <T> Buffer<T>.asReadOnly(): Buffer<T> = if (this is MutableBuffer) {
     this
 }
 
+/**
+ * Create a boxing buffer of given type
+ */
+fun <T : Any> boxingBuffer(size: Int, initializer: (Int) -> T): Buffer<T> = ListBuffer(List(size, initializer))
+
 /**
  * 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> {
+inline fun <reified T : Any> inlineBuffer(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))
+        Long::class -> LongBuffer(LongArray(size) { initializer(it) as Long }) as Buffer<T>
+        else -> boxingBuffer(size, initializer)
     }
 }
 
+/**
+ * Create a boxing mutable buffer of given type
+ */
+fun <T : Any> boxingMutableBuffer(size: Int, initializer: (Int) -> T): MutableBuffer<T> = MutableListBuffer(MutableList(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> {
+inline fun <reified T : Any> inlineMutableBuffer(size: Int, noinline initializer: (Int) -> T): MutableBuffer<T> {
     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))
+        Long::class -> LongBuffer(LongArray(size) { initializer(it) as Long }) as MutableBuffer<T>
+        else -> boxingMutableBuffer(size, initializer)
     }
-}
+}
+
+typealias BufferFactory<T> = (Int, (Int) -> T) -> Buffer<T>
+typealias MutableBufferFactory<T> = (Int, (Int) -> T) -> MutableBuffer<T>
+
+val DoubleBufferFactory: BufferFactory<Double> = { size, initializer -> DoubleBuffer(DoubleArray(size, initializer)) }
+val IntBufferFactory: BufferFactory<Int> = { size, initializer -> IntBuffer(IntArray(size, initializer)) }
+val LongBufferFactory: BufferFactory<Long> = { size, initializer -> LongBuffer(LongArray(size, initializer)) }
+

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, F : ExtendedField<N>>(shape: IntArray, field: F) : NDField<N, F>(shape, field),
+class ExtendedNDField<T : Any, F : ExtendedField<T>>(shape: IntArray, field: F) : NDField<T, F>(shape, field),
-        TrigonometricOperations<NDElement<N, F>>,
+        TrigonometricOperations<NDStructure<T>>,
-        PowerOperations<NDElement<N, F>>,
+        PowerOperations<NDStructure<T>>,
-        ExponentialOperations<NDElement<N, F>> {
+        ExponentialOperations<NDStructure<T>> {
 
-    override fun produceStructure(initializer: F.(IntArray) -> N): NDStructure<N> {
+    override fun produceStructure(initializer: F.(IntArray) -> T): NDStructure<T> {
-        return ndStructure(shape) { field.initializer(it) }
+        return NdStructure(shape, ::boxingBuffer) { field.initializer(it) }
     }
 
-    override fun power(arg: NDElement<N, F>, pow: Double): NDElement<N, F> {
+    override fun power(arg: NDStructure<T>, pow: Double): NDElement<T, F> {
-        return arg.transform { d -> with(field) { power(d, pow) } }
+        return produce { with(field) { power(arg[it], pow) } }
     }
 
-    override fun exp(arg: NDElement<N, F>): NDElement<N, F> {
+    override fun exp(arg: NDStructure<T>): NDElement<T, F> {
-        return arg.transform { d -> with(field) { exp(d) } }
+        return produce { with(field) { exp(arg[it]) } }
     }
 
-    override fun ln(arg: NDElement<N, F>): NDElement<N, F> {
+    override fun ln(arg: NDStructure<T>): NDElement<T, F> {
-        return arg.transform { d -> with(field) { ln(d) } }
+        return produce { with(field) { ln(arg[it]) } }
     }
 
-    override fun sin(arg: NDElement<N, F>): NDElement<N, F> {
+    override fun sin(arg: NDStructure<T>): NDElement<T, F> {
-        return arg.transform { d -> with(field) { sin(d) } }
+        return produce { with(field) { sin(arg[it]) } }
     }
 
-    override fun cos(arg: NDElement<N, F>): NDElement<N, F> {
+    override fun cos(arg: NDStructure<T>): NDElement<T, F> {
-        return arg.transform { d -> with(field) { cos(d) } }
+        return produce { with(field) { cos(arg[it]) } }
     }
 }
 

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, F : Field<T>>(val shape: IntArray, val field: F) : Field<NDElement<T, F>> {
+abstract class NDField<T, F : Field<T>>(val shape: IntArray, val field: F) : Field<NDStructure<T>> {
 
     abstract fun produceStructure(initializer: F.(IntArray) -> T): NDStructure<T>
 
@@ -29,10 +29,14 @@ abstract class NDField<T, F : Field<T>>(val shape: IntArray, val field: F) : Fie
         produce { zero }
     }
 
+    override val one: NDElement<T, F> by lazy {
+        produce { one }
+    }
+
     /**
      * 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, F>) {
+    private fun checkShape(vararg elements: NDStructure<T>) {
         elements.forEach {
             if (!shape.contentEquals(it.shape)) {
                 throw ShapeMismatchException(shape, it.shape)
@@ -43,7 +47,7 @@ abstract class NDField<T, F : Field<T>>(val shape: IntArray, val field: F) : Fie
     /**
      * Element-by-element addition
      */
-    override fun add(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
+    override fun add(a: NDStructure<T>, b: NDStructure<T>): NDElement<T, F> {
         checkShape(a, b)
         return produce { with(field) { a[it] + b[it] } }
     }
@@ -51,18 +55,15 @@ abstract class NDField<T, F : Field<T>>(val shape: IntArray, val field: F) : Fie
     /**
      * Multiply all elements by cinstant
      */
-    override fun multiply(a: NDElement<T, F>, k: Double): NDElement<T, F> {
+    override fun multiply(a: NDStructure<T>, k: Double): NDElement<T, F> {
         checkShape(a)
         return produce { with(field) { a[it] * k } }
     }
 
-    override val one: NDElement<T, F>
-        get() = produce { one }
-
     /**
      * Element-by-element multiplication
      */
-    override fun multiply(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
+    override fun multiply(a: NDStructure<T>, b: NDStructure<T>): NDElement<T, F> {
         checkShape(a)
         return produce { with(field) { a[it] * b[it] } }
     }
@@ -70,7 +71,7 @@ abstract class NDField<T, F : Field<T>>(val shape: IntArray, val field: F) : Fie
     /**
      * Element-by-element division
      */
-    override fun divide(a: NDElement<T, F>, b: NDElement<T, F>): NDElement<T, F> {
+    override fun divide(a: NDStructure<T>, b: NDStructure<T>): NDElement<T, F> {
         checkShape(a)
         return produce { with(field) { a[it] / b[it] } }
     }
@@ -105,7 +106,7 @@ abstract class NDField<T, F : Field<T>>(val shape: IntArray, val field: F) : Fie
 }
 
 
-interface NDElement<T, F : Field<T>>: FieldElement<NDElement<T, F>, NDField<T, F>>, NDStructure<T>
+interface NDElement<T, F : Field<T>> : FieldElement<NDStructure<T>, 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)) }
@@ -114,24 +115,24 @@ inline fun <T, F : Field<T>> NDElement<T, F>.transform(crossinline action: F.(T)
 /**
  *  Read-only [NDStructure] coupled to the context.
  */
-class NDStructureElement<T, F : Field<T>>(override val context: NDField<T, F>, private val structure: NDStructure<T>) : NDElement<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
+    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@invoke(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
     }
@@ -140,7 +141,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
     }
@@ -167,7 +168,7 @@ operator fun <T, F : Field<T>> NDElement<T, F>.div(arg: T): NDElement<T, F> = tr
 }
 
 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> = ndStructure(shape) { field.initializer(it) }
+    override fun produceStructure(initializer: F.(IntArray) -> T): NDStructure<T> = NdStructure(shape, ::boxingBuffer) { field.initializer(it) }
 }
 
 //typealias NDFieldFactory<T> = (IntArray)->NDField<T>
@@ -192,8 +193,10 @@ object NDArrays {
         return realNDArray(intArrayOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
     }
 
-    inline fun produceReal(shape: IntArray, block: ExtendedNDField<Double, DoubleField>.() -> NDElement<Double, DoubleField>) =
+    inline fun produceReal(shape: IntArray, block: ExtendedNDField<Double, DoubleField>.() -> NDStructure<Double>): NDElement<Double, DoubleField> {
-            ExtendedNDField(shape, DoubleField).run(block)
+        val field = ExtendedNDField(shape, DoubleField)
+        return NDStructureElement(field, field.run(block))
+    }
 
     /**
      * Simple boxing NDArray

kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDStructure.kt

@@ -80,7 +80,7 @@ class DefaultStrides private constructor(override val shape: IntArray) : Strides
 
     override fun offset(index: IntArray): Int {
         return index.mapIndexed { i, value ->
-            if (value < 0 || value >= shape[i]) {
+            if (value < 0 || value >= this.shape[i]) {
                 throw RuntimeException("Index $value out of shape bounds: (0,${this.shape[i]})")
             }
             value * strides[i]
@@ -138,24 +138,48 @@ class BufferNDStructure<T>(
             error("Expected buffer side of ${strides.linearSize}, but found ${buffer.size}")
         }
     }
+
+    override fun equals(other: Any?): Boolean {
+        return when {
+            this === other -> true
+            other is BufferNDStructure<*> -> this.strides == other.strides && this.buffer.contentEquals(other.buffer)
+            other is NDStructure<*> -> elements().all { (index, value) -> value == other[index] }
+            else -> false
+        }
+    }
+
+    override fun hashCode(): Int {
+        var result = strides.hashCode()
+        result = 31 * result + buffer.hashCode()
+        return result
+    }
+
 }
 
 /**
- * Create a most suitable nd-structure avoiding boxing if possible using given strides.
+ * Create a NDStructure with explicit buffer factory
  *
  * Strides should be reused if possible
  */
-inline fun <reified T : Any> ndStructure(strides: Strides, noinline initializer: (IntArray) -> T) =
+@Suppress("FunctionName")
-        BufferNDStructure<T>(strides, buffer(strides.linearSize) { i -> initializer(strides.index(i)) })
+fun <T : Any> NdStructure(strides: Strides, bufferFactory: BufferFactory<T>, initializer: (IntArray) -> T) =
+        BufferNDStructure(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
 
 /**
- * Create a most suitable nd-structure avoiding boxing if possible using default strides with given shape
+ * Inline create NDStructure with non-boxing buffer implementation if it is possible
  */
-inline fun <reified T : Any> ndStructure(shape: IntArray, noinline initializer: (IntArray) -> T) =
+inline fun <reified T : Any> inlineNDStructure(strides: Strides, crossinline initializer: (IntArray) -> T) =
-        ndStructure(DefaultStrides(shape), initializer)
+        BufferNDStructure(strides, inlineBuffer(strides.linearSize) { i -> initializer(strides.index(i)) })
+
+@Suppress("FunctionName")
+fun <T : Any> NdStructure(shape: IntArray, bufferFactory: BufferFactory<T>, initializer: (IntArray) -> T) =
+        NdStructure(DefaultStrides(shape), bufferFactory, initializer)
+
+inline fun <reified T : Any> inlineNdStructure(shape: IntArray, crossinline initializer: (IntArray) -> T) =
+        inlineNDStructure(DefaultStrides(shape), initializer)
 
 /**
- * Mutable ND buffer based on linear [Buffer]
+ * Mutable ND buffer based on linear [inlineBuffer]
  */
 class MutableBufferNDStructure<T>(
         override val strides: Strides,
@@ -172,13 +196,27 @@ class MutableBufferNDStructure<T>(
 }
 
 /**
- * The same as [ndStructure], but mutable
+ * The same as [inlineNDStructure], but mutable
  */
-inline fun <reified T : Any> mutableNdStructure(strides: Strides, noinline initializer: (IntArray) -> T) =
+@Suppress("FunctionName")
-        MutableBufferNDStructure(strides, mutableBuffer(strides.linearSize) { i -> initializer(strides.index(i)) })
+fun <T : Any> MutableNdStructure(strides: Strides, bufferFactory: MutableBufferFactory<T>, initializer: (IntArray) -> T) =
+        MutableBufferNDStructure(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
+
+inline fun <reified T : Any> inlineMutableNdStructure(strides: Strides, crossinline initializer: (IntArray) -> T) =
+        MutableBufferNDStructure(strides, inlineMutableBuffer(strides.linearSize) { i -> initializer(strides.index(i)) })
+
+@Suppress("FunctionName")
+fun <T : Any> MutableNdStructure(shape: IntArray, bufferFactory: MutableBufferFactory<T>, initializer: (IntArray) -> T) =
+        MutableNdStructure(DefaultStrides(shape), bufferFactory, initializer)
+
+inline fun <reified T : Any> inlineMutableNdStructure(shape: IntArray, crossinline initializer: (IntArray) -> T) =
+        inlineMutableNdStructure(DefaultStrides(shape), initializer)
+
+inline fun <reified T : Any> NDStructure<T>.combine(struct: NDStructure<T>, crossinline block: (T, T) -> T): NDStructure<T> {
+    if (!this.shape.contentEquals(struct.shape)) error("Shape mismatch in structure combination")
+    return inlineNdStructure(shape) { block(this[it], struct[it]) }
+}
 
-inline fun <reified T : Any> mutableNdStructure(shape: IntArray, noinline initializer: (IntArray) -> T) =
-        mutableNdStructure(DefaultStrides(shape), initializer)
 
 ///**
 // * Create universal mutable structure

kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/ArrayMatrixTest.kt

@@ -1,34 +0,0 @@
-package scientifik.kmath.linear
-
-import kotlin.test.Test
-import kotlin.test.assertEquals
-
-class ArrayMatrixTest {
-
-    @Test
-    fun testSum() {
-        val vector1 = Vector.ofReal(5) { it.toDouble() }
-        val vector2 = Vector.ofReal(5) { 5 - it.toDouble() }
-        val sum = vector1 + vector2
-        assertEquals(5.0, sum[2])
-    }
-
-    @Test
-    fun testVectorToMatrix() {
-        val vector = Vector.ofReal(5) { it.toDouble() }
-        val matrix = vector.toMatrix()
-        assertEquals(4.0, matrix[4, 0])
-    }
-
-
-    @Test
-    fun testDot() {
-        val vector1 = Vector.ofReal(5) { it.toDouble() }
-        val vector2 = Vector.ofReal(5) { 5 - it.toDouble() }
-        val product = vector1.toMatrix() dot (vector2.toMatrix().transpose())
-
-
-        assertEquals(5.0, product[1, 0])
-        assertEquals(6.0, product[2, 2])
-    }
-}

kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/MatrixTest.kt

@@ -0,0 +1,46 @@
+package scientifik.kmath.linear
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class MatrixTest {
+
+    @Test
+    fun testSum() {
+        val vector1 = Vector.real(5) { it.toDouble() }
+        val vector2 = Vector.real(5) { 5 - it.toDouble() }
+        val sum = vector1 + vector2
+        assertEquals(5.0, sum[2])
+    }
+
+    @Test
+    fun testVectorToMatrix() {
+        val vector = Vector.real(5) { it.toDouble() }
+        val matrix = vector.toMatrix()
+        assertEquals(4.0, matrix[4, 0])
+    }
+
+    @Test
+    fun testTranspose(){
+        val matrix = MatrixSpace.real(3,3).one
+        val transposed = matrix.transpose()
+        assertEquals(matrix.context, transposed.context)
+        assertEquals((matrix as StructureMatrix).structure, (transposed as StructureMatrix).structure)
+        assertEquals(matrix, transposed)
+    }
+
+
+    @Test
+    fun testDot() {
+        val vector1 = Vector.real(5) { it.toDouble() }
+        val vector2 = Vector.real(5) { 5 - it.toDouble() }
+
+        val matrix1 = vector1.toMatrix()
+        val matrix2 = vector2.toMatrix().transpose()
+        val product = matrix1 dot matrix2
+
+
+        assertEquals(5.0, product[1, 0])
+        assertEquals(6.0, product[2, 2])
+    }
+}

kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt

@@ -1,15 +1,14 @@
 package scientifik.kmath.linear
 
-import scientifik.kmath.operations.DoubleField
 import kotlin.test.Test
-import kotlin.test.assertTrue
+import kotlin.test.assertEquals
 
 class RealLUSolverTest {
     @Test
     fun testInvertOne() {
-        val matrix = Matrix.diagonal(2, 2, DoubleField)
+        val matrix = MatrixSpace.real(2,2).one
         val inverted = RealLUSolver.inverse(matrix)
-        assertTrue { Matrix.equals(matrix,inverted) }
+        assertEquals(matrix,inverted)
     }
 
 //    @Test

kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt

@@ -6,7 +6,6 @@ import kotlin.test.assertEquals
 class RealFieldTest {
     @Test
     fun testSqrt() {
-
         //fails because KT-27586
         val sqrt = with(RealField) {
             sqrt(  25 * one)

kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/NumberNDFieldTest.kt

@@ -51,9 +51,14 @@ class NumberNDFieldTest {
         assertEquals(2.0, result[0, 2])
     }
 
-    object L2Norm : Norm<NDElement<out Number, *>, Double> {
+    @Test
-        override fun norm(arg: NDElement<out Number, *>): Double {
+    fun combineTest(){
-            return kotlin.math.sqrt(arg.sumByDouble { it.second.toDouble() })
+        val division = array1.combine(array2, Double::div)
+    }
+
+    object L2Norm : Norm<NDStructure<out Number>, Double> {
+        override fun norm(arg: NDStructure<out Number>): Double {
+            return kotlin.math.sqrt(arg.elements().sumByDouble { it.second.toDouble() })
         }
     }