.gitignore

@@ -8,3 +8,4 @@
 # Cache of project
 .gradletasknamecache
 
+gradle.properties

README.md

@@ -9,13 +9,27 @@ and `scipy` it is modular and has a lightweight core.
     * Basic linear algebra operations (summs products, etc) backed by `Space` API.
     * Complex numbers backed by `Field` API (meaning that they will be useable in any structures like vectors and NDArrays).
     * [In progress] advanced linear algebra operations like matrix inversions.
-* **Array-like structures** Full support of numpy-like ndarray including mixed ariphmetic operations and function operations
+* **Array-like structures** Full support of numpy-like ndarray including mixed arithmetic operations and function operations
 on arrays and numbers just like it works in python (with benefit of static type checking).
 
+* **Expressions** Expressions are one of the ultimate goals of kmath. It is planned to be able to write some mathematical
+expression once an then apply it to different types of objects by providing different context. Exception could be used
+for a wide variety of purposes from high performance calculations to code generation.
+
+## Planned features
+
+* **Common mathematics** It is planned to gradually wrap most parts of [Apache commons-math](http://commons.apache.org/proper/commons-math/) 
+library in kotlin code and maybe rewrite some parts to better suite kotlin programming paradigm. There is no fixed priority list for that. Feel free
+to submit a future request if you want something to be done first.
+
+* **Messaging** A mathematical notation to support multi-language and multi-node communication for mathematical tasks.
+
 ## Multi-platform support
 KMath is developed as a multi-platform library, which means that most of interfaces are declared in common module.
 Implementation is also done in common module wherever it is possible. In some cases features are delegated to 
 platform even if they could be done in common module because of platform performance optimization. 
+Currently the main focus of development is the JVM platform, contribution of implementations for Kotlin - Native and
+Kotlin - JS is welcome. 
 
 ## Performance
 The calculation performance is one of major goals of KMath in the future, but in some cases it is not possible to achieve 
@@ -25,8 +39,33 @@ but worse than optimized native/scipy (mostly due to boxing operations on primit
 of optimized parts should be better than scipy.
 
 ## Releases
-The project is currently in pre-release stage. Work builds could be obtained with 
-[![](https://jitpack.io/v/altavir/kmath.svg)](https://jitpack.io/#altavir/kmath). 
+
+The project is currently in pre-release stage. Nightly builds could be used by adding additional repository to (groovy) gradle config:
+```groovy
+repositories {
+    maven { url = "http://npm.mipt.ru:8081/artifactory/gradle-dev" }
+    mavenCentral()
+} 
+```
+or for kotlin gradle dsl:
+
+```kotlin
+repositories {
+    maven { setUrl("http://npm.mipt.ru:8081/artifactory/gradle-dev") }
+    mavenCentral()
+} 
+```
+
+Then use regular dependency like
+```groovy
+compile(group: 'scientifik', name: 'kmath-core-jvm', version: '0.0.1-SNAPSHOT')
+```
+or in kotlin
+```kotlin
+compile(group = "scientifik", name = "kmath-core-jvm", version = "0.0.1-SNAPSHOT")
+```
+
+Work builds could be obtained with [![](https://jitpack.io/v/altavir/kmath.svg)](https://jitpack.io/#altavir/kmath). 
 
 ## Contributing
 The project requires a lot of additional work. Please fill free to contribute in any way and propose new features.

build.gradle

@@ -1,14 +1,49 @@
 buildscript {
-    ext.kotlin_version = '1.2.60'
+    ext.kotlin_version = '1.3.0-rc-146'
 
     repositories {
-        mavenCentral()
+        jcenter()
+        maven {
+            url = "http://dl.bintray.com/kotlin/kotlin-eap"
+        }
     }
+
     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.1 - SNAPSHOT'
+    group = 'scientifik'
+    version = '0.0.1-SNAPSHOT'
+}
 
+
+artifactory {
+    contextUrl = "${artifactory_contextUrl}"   //The base Artifactory URL if not overridden by the publisher/resolver
+    publish {
+        repository {
+            repoKey = 'gradle-dev-local'
+            username = "${artifactory_user}"
+            password = "${artifactory_password}"
+        }
+
+        defaults {
+            publications('jvm', 'js', 'kotlinMultiplatform', 'metadata')
+            publishBuildInfo = false
+            publishArtifacts = true
+            publishPom = true
+            publishIvy = false
+        }
+    }
+    resolve {
+        repository {
+            repoKey = 'gradle-dev'
+            username = "${artifactory_user}"
+            password = "${artifactory_password}"
+        }
+    }
+}

kmath-common/build.gradle

@@ -1,14 +0,0 @@
-description = "Platform-independent interfaces for kotlin maths"
-
-apply plugin: 'kotlin-platform-common'
-
-repositories {
-    mavenCentral()
-}
-
-dependencies {
-    compile "org.jetbrains.kotlin:kotlin-stdlib-common:$kotlin_version"
-    testCompile "org.jetbrains.kotlin:kotlin-test-annotations-common:$kotlin_version"
-    testCompile "org.jetbrains.kotlin:kotlin-test-common:$kotlin_version"
-}
-

kmath-common/src/main/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt

@@ -1,144 +0,0 @@
-package scientifik.kmath.structures
-
-import scientifik.kmath.operations.DoubleField
-import scientifik.kmath.operations.Field
-import scientifik.kmath.operations.Space
-import scientifik.kmath.operations.SpaceElement
-
-/**
- * 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 LinearSpace<T : Any, V : LinearStructure<out T>>(val rows: Int, val columns: Int, val field: Field<T>) : Space<V> {
-
-    /**
-     * Produce the element of this space
-     */
-    abstract fun produce(initializer: (Int, Int) -> T): V
-
-    /**
-     * Produce new linear space with given dimensions
-     */
-    abstract fun produceSpace(rows: Int, columns: Int): LinearSpace<T, V>
-
-    override val zero: V by lazy {
-        produce { _, _ -> field.zero }
-    }
-
-    override fun add(a: V, b: V): V {
-        return produce { i, j -> with(field) { a[i, j] + b[i, j] } }
-    }
-
-    override fun multiply(a: V, k: Double): V {
-        //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
-     */
-    fun multiply(a: V, b: V): V {
-        if (a.rows != b.columns) {
-            //TODO replace by specific exception
-            error("Dimension mismatch in vector dot product")
-        }
-        return produceSpace(a.rows, b.columns).produce { i, j ->
-            (0..a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) }
-        }
-    }
-
-    infix fun V.dot(b: V): V = multiply(this, b)
-}
-
-/**
- * A matrix-like structure that is not dependent on specific space implementation
- */
-interface LinearStructure<T : Any> {
-    val rows: Int
-    val columns: Int
-
-    operator fun get(i: Int, j: Int): T
-
-    fun transpose(): LinearStructure<T> {
-        return object : LinearStructure<T> {
-            override val rows: Int = this@LinearStructure.columns
-            override val columns: Int = this@LinearStructure.rows
-            override fun get(i: Int, j: Int): T = this@LinearStructure.get(j, i)
-        }
-    }
-}
-
-interface Vector<T : Any> : LinearStructure<T> {
-    override val columns: Int
-        get() = 1
-
-    operator fun get(i: Int) = get(i, 0)
-}
-
-
-/**
- * DoubleArray-based implementation of vector space
- */
-class ArraySpace<T : Any>(rows: Int, columns: Int, field: Field<T>) : LinearSpace<T, LinearStructure<out T>>(rows, columns, field) {
-
-    override fun produce(initializer: (Int, Int) -> T): LinearStructure<T> = ArrayMatrix<T>(this, initializer)
-
-
-    override fun produceSpace(rows: Int, columns: Int): LinearSpace<T, LinearStructure<out T>> {
-        return ArraySpace(rows, columns, field)
-    }
-}
-
-/**
- * Member of [ArraySpace] which wraps 2-D array
- */
-class ArrayMatrix<T : Any>(override val context: ArraySpace<T>, initializer: (Int, Int) -> T) : LinearStructure<T>, SpaceElement<LinearStructure<out T>, ArraySpace<T>> {
-
-    val list: List<List<T>> = (0 until rows).map { i -> (0 until columns).map { j -> initializer(i, j) } }
-
-    override val rows: Int get() = context.rows
-
-    override val columns: Int get() = context.columns
-
-    override fun get(i: Int, j: Int): T {
-        return list[i][j]
-    }
-
-    override val self: ArrayMatrix<T> get() = this
-}
-
-
-class ArrayVector<T : Any>(override val context: ArraySpace<T>, initializer: (Int) -> T) : Vector<T>, SpaceElement<LinearStructure<out T>, ArraySpace<T>> {
-
-    init {
-        if (context.columns != 1) {
-            error("Vector must have single column")
-        }
-    }
-
-    val list: List<T> = (0 until context.rows).map(initializer)
-
-
-    override val rows: Int get() = context.rows
-
-    override val columns: Int = 1
-
-    override fun get(i: Int, j: Int): T {
-        return list[i]
-    }
-
-    override val self: ArrayVector<T> get() = this
-
-}
-
-fun <T : Any> vector(size: Int, field: Field<T>, initializer: (Int) -> T) = ArrayVector(ArraySpace(size, 1, field), initializer)
-//TODO replace by primitive array version
-fun realVector(size: Int, initializer: (Int) -> Double) = vector(size, DoubleField, initializer)
-
-fun <T : Any> Array<T>.asVector(field: Field<T>) = vector(size, field) { this[it] }
-//TODO add inferred field from field element
-fun DoubleArray.asVector() = realVector(this.size) { this[it] }
-
-fun <T : Any> matrix(rows: Int, columns: Int, field: Field<T>, initializer: (Int, Int) -> T) = ArrayMatrix<T>(ArraySpace(rows, columns, field), initializer)
-fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) = matrix(rows, columns, DoubleField, initializer)

kmath-core/build.gradle

@@ -0,0 +1,58 @@
+plugins {
+    id 'kotlin-multiplatform'// version '1.3.0-rc-116'
+    id "me.champeau.gradle.jmh" version "0.4.5"
+}
+
+repositories {
+    maven { url = 'http://dl.bintray.com/kotlin/kotlin-eap' }
+    mavenCentral()
+}
+
+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 {
+                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 {
+//        }
+    }
+}

kmath-core/src/commonMain/kotlin/scientifik/kmath/expressions/Expression.kt

@@ -0,0 +1,62 @@
+package scientifik.kmath.expressions
+
+import scientifik.kmath.operations.Field
+import scientifik.kmath.operations.Space
+
+
+interface Expression<T> {
+    operator fun invoke(arguments: Map<String, T>): T
+}
+
+operator fun <T> Expression<T>.invoke(vararg pairs: Pair<String, T>): T = invoke(mapOf(*pairs))
+
+interface ExpressionContext<T> {
+    fun variable(name: String, default: T? = null): Expression<T>
+
+    fun const(value: T): Expression<T>
+}
+
+internal class VariableExpression<T>(val name: String, val default: T? = null) : Expression<T> {
+    override fun invoke(arguments: Map<String, T>): T {
+        return arguments[name] ?: default ?: error("The parameter not found: $name")
+    }
+}
+
+internal class ConstantExpression<T>(val value: T) : Expression<T> {
+    override fun invoke(arguments: Map<String, T>): T = value
+}
+
+internal class SumExpression<T>(val context: Space<T>, val first: Expression<T>, val second: Expression<T>) : Expression<T> {
+    override fun invoke(arguments: Map<String, T>): T = context.add(first.invoke(arguments), second.invoke(arguments))
+}
+
+internal class ProductExpression<T>(val context: Field<T>, val first: Expression<T>, val second: Expression<T>) : Expression<T> {
+    override fun invoke(arguments: Map<String, T>): T = context.multiply(first.invoke(arguments), second.invoke(arguments))
+}
+
+internal class ConstProductExpession<T>(val context: Field<T>, val expr: Expression<T>, val const: Double) : Expression<T> {
+    override fun invoke(arguments: Map<String, T>): T = context.multiply(expr.invoke(arguments), const)
+}
+
+internal class DivExpession<T>(val context: Field<T>, val expr: Expression<T>, val second: Expression<T>) : Expression<T> {
+    override fun invoke(arguments: Map<String, T>): T = context.divide(expr.invoke(arguments), second.invoke(arguments))
+}
+
+class FieldExpressionContext<T>(val field: Field<T>) : Field<Expression<T>>, ExpressionContext<T> {
+
+    override val zero: Expression<T> = ConstantExpression(field.zero)
+
+    override val one: Expression<T> = ConstantExpression(field.one)
+
+    override fun const(value: T): Expression<T> = ConstantExpression(value)
+
+    override fun variable(name: String, default: T?): Expression<T> = VariableExpression(name, default)
+
+    override fun add(a: Expression<T>, b: Expression<T>): Expression<T> = SumExpression(field, a, b)
+
+    override fun multiply(a: Expression<T>, k: Double): Expression<T> = ConstProductExpession(field, a, k)
+
+    override fun multiply(a: Expression<T>, b: Expression<T>): Expression<T> = ProductExpression(field, a, b)
+
+    override fun divide(a: Expression<T>, b: Expression<T>): Expression<T> = DivExpession(field, a, b)
+}

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

@@ -0,0 +1,294 @@
+package scientifik.kmath.linear
+
+import scientifik.kmath.structures.MutableNDArray
+import scientifik.kmath.structures.NDArray
+import scientifik.kmath.structures.NDArrays
+import kotlin.math.absoluteValue
+
+/**
+ * Implementation copier from Apache common-maths
+ */
+abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
+
+    private val field get() = matrix.context.field
+    /** Entries of LU decomposition.  */
+    internal val lu: NDArray<T>
+    /** Pivot permutation associated with LU decomposition.  */
+    internal val pivot: IntArray
+    /** Parity of the permutation associated with the LU decomposition.  */
+    private var even: Boolean = false
+
+    init {
+        val pair = calculateLU()
+        lu = pair.first
+        pivot = pair.second
+    }
+
+    /**
+     * Returns the matrix L of the decomposition.
+     *
+     * L is a lower-triangular matrix
+     * @return the L matrix (or null if decomposed matrix is singular)
+     */
+    val l: Matrix<out T> by lazy {
+        matrix.context.produce { i, j ->
+            when {
+                j < i -> lu[i, j]
+                j == i -> matrix.context.field.one
+                else -> matrix.context.field.zero
+            }
+        }
+    }
+
+
+    /**
+     * Returns the matrix U of the decomposition.
+     *
+     * U is an upper-triangular matrix
+     * @return the U matrix (or null if decomposed matrix is singular)
+     */
+    val u: Matrix<out T> by lazy {
+        matrix.context.produce { i, j ->
+            if (j >= i) lu[i, j] else field.zero
+        }
+    }
+
+    /**
+     * Returns the P rows permutation matrix.
+     *
+     * P is a sparse matrix with exactly one element set to 1.0 in
+     * each row and each column, all other elements being set to 0.0.
+     *
+     * The positions of the 1 elements are given by the [ pivot permutation vector][.getPivot].
+     * @return the P rows permutation matrix (or null if decomposed matrix is singular)
+     * @see .getPivot
+     */
+    val p: Matrix<out T> by lazy {
+        matrix.context.produce { i, j ->
+            //TODO ineffective. Need sparse matrix for that
+            if (j == pivot[i]) field.one else field.zero
+        }
+    }
+
+    /**
+     * Return the determinant of the matrix
+     * @return determinant of the matrix
+     */
+    val determinant: T
+        get() {
+            with(matrix.context.field) {
+                var determinant = if (even) one else -one
+                for (i in 0 until matrix.rows) {
+                    determinant *= lu[i, i]
+                }
+                return determinant
+            }
+        }
+
+//    /**
+//     * Get a solver for finding the A  X = B solution in exact linear
+//     * sense.
+//     * @return a solver
+//     */
+//    val solver: DecompositionSolver
+//        get() = Solver(lu, pivot, singular)
+
+    /**
+     * In-place transformation for [MutableNDArray], using given transformation for each element
+     */
+    operator fun <T> MutableNDArray<T>.set(i: Int, j: Int, value: T) {
+        this[listOf(i, j)] = value
+    }
+
+    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 calculateLU(): Pair<NDArray<T>, IntArray> {
+        if (matrix.rows != matrix.columns) {
+            error("LU decomposition supports only square matrices")
+        }
+
+        val m = matrix.columns
+        val pivot = IntArray(matrix.rows)
+        //TODO fix performance
+        val lu: MutableNDArray<T> = NDArrays.createMutable(matrix.context.field, listOf(matrix.rows, matrix.columns)) { index -> matrix[index[0], index[1]] }
+
+
+        with(matrix.context.field) {
+            // Initialize permutation array and parity
+            for (row in 0 until m) {
+                pivot[row] = row
+            }
+            even = true
+
+            // Loop over columns
+            for (col in 0 until m) {
+
+                // upper
+                for (row in 0 until col) {
+                    var sum = lu[row, col]
+                    for (i in 0 until row) {
+                        sum -= lu[row, i] * lu[i, col]
+                    }
+                    lu[row, col] = sum
+                }
+
+                // lower
+                val max = (col until m).maxBy { row ->
+                    var sum = lu[row, col]
+                    for (i in 0 until col) {
+                        sum -= lu[row, i] * lu[i, col]
+                    }
+                    //luRow[col] = sum
+                    lu[row, col] = sum
+
+                    abs(sum)
+                } ?: col
+
+                // Singularity check
+                if (isSingular(lu[max, col])) {
+                    error("Singular matrix")
+                }
+
+                // Pivot if necessary
+                if (max != col) {
+                    //var tmp = zero
+                    //val luMax = lu[max]
+                    //val luCol = lu[col]
+                    for (i in 0 until m) {
+                        lu[max, i] = lu[col, i]
+                        lu[col, i] = lu[max, i]
+                    }
+                    val temp = pivot[max]
+                    pivot[max] = pivot[col]
+                    pivot[col] = temp
+                    even = !even
+                }
+
+                // Divide the lower elements by the "winning" diagonal elt.
+                val luDiag = lu[col, col]
+                for (row in col + 1 until m) {
+                    lu[row, col] = lu[row, col] / luDiag
+//                    lu[row, col] /= luDiag
+                }
+            }
+        }
+        return Pair(lu, pivot)
+    }
+
+    /**
+     * Returns the pivot permutation vector.
+     * @return the pivot permutation vector
+     * @see .getP
+     */
+    fun getPivot(): IntArray {
+        return pivot.copyOf()
+    }
+
+}
+
+class RealLUDecomposition(matrix: Matrix<Double>, private val singularityThreshold: Double = DEFAULT_TOO_SMALL) : LUDecomposition<Double>(matrix) {
+    override fun isSingular(value: Double): Boolean {
+        return value.absoluteValue < singularityThreshold
+    }
+
+    companion object {
+        /** Default bound to determine effective singularity in LU decomposition.  */
+        private const val DEFAULT_TOO_SMALL = 1e-11
+    }
+}
+
+
+/** Specialized solver.  */
+object RealLUSolver : LinearSolver<Double> {
+
+//
+//    /** {@inheritDoc}  */
+//    override fun solve(b: RealVector): RealVector {
+//        val m = pivot.size
+//        if (b.getDimension() != m) {
+//            throw DimensionMismatchException(b.getDimension(), m)
+//        }
+//        if (singular) {
+//            throw SingularMatrixException()
+//        }
+//
+//        val bp = DoubleArray(m)
+//
+//        // Apply permutations to b
+//        for (row in 0 until m) {
+//            bp[row] = b.getEntry(pivot[row])
+//        }
+//
+//        // Solve LY = b
+//        for (col in 0 until m) {
+//            val bpCol = bp[col]
+//            for (i in col + 1 until m) {
+//                bp[i] -= bpCol * lu[i][col]
+//            }
+//        }
+//
+//        // Solve UX = Y
+//        for (col in m - 1 downTo 0) {
+//            bp[col] /= lu[col][col]
+//            val bpCol = bp[col]
+//            for (i in 0 until col) {
+//                bp[i] -= bpCol * lu[i][col]
+//            }
+//        }
+//
+//        return ArrayRealVector(bp, false)
+//    }
+
+
+    fun decompose(mat: Matrix<Double>, threshold: Double = 1e-11): RealLUDecomposition = RealLUDecomposition(mat, threshold)
+
+    override fun solve(a: Matrix<Double>, b: Matrix<Double>): Matrix<Double> {
+        val decomposition = decompose(a, a.context.field.zero)
+
+        if (b.rows != a.rows) {
+            error("Matrix dimension mismatch expected ${a.rows}, but got ${b.rows}")
+        }
+
+        // Apply permutations to b
+        val bp = Array(a.rows) { DoubleArray(b.columns) }
+        for (row in 0 until a.rows) {
+            val bpRow = bp[row]
+            val pRow = decomposition.pivot[row]
+            for (col in 0 until b.columns) {
+                bpRow[col] = b[pRow, col]
+            }
+        }
+
+        // Solve LY = b
+        for (col in 0 until a.rows) {
+            val bpCol = bp[col]
+            for (i in col + 1 until a.rows) {
+                val bpI = bp[i]
+                val luICol = decomposition.lu[i, col]
+                for (j in 0 until b.columns) {
+                    bpI[j] -= bpCol[j] * luICol
+                }
+            }
+        }
+
+        // Solve UX = Y
+        for (col in a.rows - 1 downTo 0) {
+            val bpCol = bp[col]
+            val luDiag = decomposition.lu[col, col]
+            for (j in 0 until b.columns) {
+                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) {
+                    bpI[j] -= bpCol[j] * luICol
+                }
+            }
+        }
+
+        return a.context.produce { i, j -> bp[i][j] }
+    }
+}

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

@@ -0,0 +1,311 @@
+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.structures.NDArray
+import scientifik.kmath.structures.NDArrays.createFactory
+import scientifik.kmath.structures.NDFieldFactory
+import scientifik.kmath.structures.realNDFieldFactory
+
+/**
+ * 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
+        }
+    }
+}
+
+
+/**
+ * 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>> {
+    val size: Int
+        get() = context.size
+
+    operator fun get(i: Int): T
+
+    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 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
+        }
+    }
+}
+
+
+/**
+ * NDArray-based implementation of vector space. By default uses slow [SimpleNDField], but could be overridden with custom [NDField] factory.
+ */
+class ArrayMatrixSpace<T : Any>(
+        rows: Int,
+        columns: Int,
+        field: Field<T>,
+        val ndFactory: NDFieldFactory<T> = createFactory(field)
+) : MatrixSpace<T>(rows, columns, field) {
+
+    val ndField by lazy {
+        ndFactory(listOf(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)
+    }
+}
+
+class ArrayVectorSpace<T : Any>(
+        size: Int,
+        field: Field<T>,
+        val ndFactory: NDFieldFactory<T> = createFactory(field)
+) : VectorSpace<T>(size, field) {
+    val ndField by lazy {
+        ndFactory(listOf(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 array: NDArray<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 array[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> {
+
+    constructor(context: ArrayVectorSpace<T>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
+
+    init {
+        if (context.size != array.shape[0]) {
+            error("Array dimension mismatch")
+        }
+    }
+
+    override fun get(i: Int): T {
+        return array[i]
+    }
+
+    override val self: ArrayVector<T> get() = this
+}
+
+/**
+ * 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))
+}
+
+/**
+ * Convert vector to array (copying content of array)
+ */
+fun <T : Any> Array<T>.toVector(field: Field<T>) = Vector.of(size, field) { this[it] }
+
+fun DoubleArray.toVector() = Vector.ofReal(this.size) { this[it] }
+
+/**
+ * Convert matrix to vector if it is possible
+ */
+fun <T : Any> Matrix<T>.toVector(): Vector<T> {
+    return when {
+        this.columns == 1 -> {
+//            if (this is ArrayMatrix) {
+//                //Reuse existing underlying array
+//                ArrayVector(ArrayVectorSpace(rows, context.field, context.ndFactory), array)
+//            } else {
+//                //Generic vector
+//                vector(rows, context.field) { get(it, 0) }
+//            }
+            Vector.of(rows, context.field) { get(it, 0) }
+        }
+        else -> error("Can't convert matrix with more than one column to vector")
+    }
+}
+
+fun <T : Any> Vector<T>.toMatrix(): Matrix<T> {
+//    return if (this is ArrayVector) {
+//        //Reuse existing underlying array
+//        ArrayMatrix(ArrayMatrixSpace(size, 1, context.field, context.ndFactory), array)
+//    } else {
+//        //Generic vector
+//        matrix(size, 1, context.field) { i, j -> get(i) }
+//    }
+    return Matrix.of(size, 1, context.field) { i, j -> get(i) }
+}
+

kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Cumulative.kt

@@ -0,0 +1,59 @@
+package scientifik.kmath.misc
+
+import kotlin.jvm.JvmName
+
+
+/**
+ * Generic cumulative operation on iterator
+ * @param T type of initial iterable
+ * @param R type of resulting iterable
+ * @param initial lazy evaluated
+ */
+fun <T, R> Iterator<T>.cumulative(initial: R, operation: (T, R) -> R): Iterator<R> = object : Iterator<R> {
+    var state: R = initial
+    override fun hasNext(): Boolean = this@cumulative.hasNext()
+
+    override fun next(): R {
+        state = operation.invoke(this@cumulative.next(), state)
+        return state
+    }
+}
+
+fun <T, R> Iterable<T>.cumulative(initial: R, operation: (T, R) -> R): Iterable<R> = object : Iterable<R> {
+    override fun iterator(): Iterator<R> = this@cumulative.iterator().cumulative(initial, operation)
+}
+
+fun <T, R> Sequence<T>.cumulative(initial: R, operation: (T, R) -> R): Sequence<R> = object : Sequence<R> {
+    override fun iterator(): Iterator<R> = this@cumulative.iterator().cumulative(initial, operation)
+}
+
+fun <T, R> List<T>.cumulative(initial: R, operation: (T, R) -> R): List<R> = this.iterator().cumulative(initial, operation).asSequence().toList()
+
+//Cumulative sum
+
+@JvmName("cumulativeSumOfDouble")
+fun Iterable<Double>.cumulativeSum() = this.cumulative(0.0){ element, sum -> sum + element}
+
+@JvmName("cumulativeSumOfInt")
+fun Iterable<Int>.cumulativeSum() = this.cumulative(0){ element, sum -> sum + element}
+
+@JvmName("cumulativeSumOfLong")
+fun Iterable<Long>.cumulativeSum() = this.cumulative(0L){ element, sum -> sum + element}
+
+@JvmName("cumulativeSumOfDouble")
+fun Sequence<Double>.cumulativeSum() = this.cumulative(0.0){ element, sum -> sum + element}
+
+@JvmName("cumulativeSumOfInt")
+fun Sequence<Int>.cumulativeSum() = this.cumulative(0){ element, sum -> sum + element}
+
+@JvmName("cumulativeSumOfLong")
+fun Sequence<Long>.cumulativeSum() = this.cumulative(0L){ element, sum -> sum + element}
+
+@JvmName("cumulativeSumOfDouble")
+fun List<Double>.cumulativeSum() = this.cumulative(0.0){ element, sum -> sum + element}
+
+@JvmName("cumulativeSumOfInt")
+fun List<Int>.cumulativeSum() = this.cumulative(0){ element, sum -> sum + element}
+
+@JvmName("cumulativeSumOfLong")
+fun List<Long>.cumulativeSum() = this.cumulative(0L){ element, sum -> sum + element}

kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Grids.kt

@@ -0,0 +1,37 @@
+package scientifik.kmath.misc
+
+/**
+ * Convert double range to sequence.
+ *
+ * If the step is positive, than the sequence starts with the lower boundary and increments by [step] until current value is lower than upper boundary.
+ * The boundary itself is not necessary included.
+ *
+ * If step is negative, the same goes from upper boundary downwards
+ */
+fun ClosedFloatingPointRange<Double>.toSequence(step: Double): Sequence<Double> {
+    return when {
+        step == 0.0 -> error("Zero step in double progression")
+        step > 0 -> sequence {
+            var current = start
+            while (current <= endInclusive) {
+                yield(current)
+                current += step
+            }
+        }
+        else -> sequence {
+            var current = endInclusive
+            while (current >= start) {
+                yield(current)
+                current += step
+            }
+        }
+    }
+}
+
+/**
+ * Convert double range to array of evenly spaced doubles, where the size of array equals [numPoints]
+ */
+fun ClosedFloatingPointRange<Double>.toGrid(numPoints: Int): DoubleArray {
+    if (numPoints < 2) error("Can't create grid with less than two points")
+    return DoubleArray(numPoints) { i -> start + (endInclusive - start) / (numPoints - 1) * i }
+}

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

@@ -3,8 +3,15 @@ package scientifik.kmath.operations
 
 /**
  * The generic mathematics elements which is able to store its context
+ * @param T the self type of the element
+ * @param S the type of mathematical context for this element
  */
-interface MathElement<T, S>{
+interface MathElement<T, S> {
+    /**
+     * Self value. Needed for static type checking.
+     */
+    val self: T
+
     /**
      * The context this element belongs to
      */
@@ -51,13 +58,7 @@ interface Space<T> {
  * @param T  self type of the element. Needed for static type checking
  * @param S the type of space
  */
-interface SpaceElement<T, S : Space<T>>: MathElement<T,S> {
-
-    /**
-     * Self value. Needed for static type checking. Needed to avoid type erasure on JVM.
-     */
-    val self: T
-
+interface SpaceElement<T, S : Space<T>> : MathElement<T, S> {
     operator fun plus(b: T): T = context.add(self, b)
     operator fun minus(b: T): T = context.add(self, context.multiply(b, -1.0))
     operator fun times(k: Number): T = context.multiply(self, k.toDouble())
@@ -99,6 +100,12 @@ interface Field<T> : Ring<T> {
 
     operator fun T.div(b: T): T = divide(this, b)
     operator fun Number.div(b: T) = this * divide(one, b)
+
+    operator fun T.plus(b: Number) = this.plus(b * one)
+    operator fun Number.plus(b: T) = b + this
+
+    operator fun T.minus(b: Number) = this.minus(b * one)
+    operator fun Number.minus(b: T) = -b + this
 }
 
 /**

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

@@ -0,0 +1,41 @@
+package scientifik.kmath.operations
+
+/**
+ * A field for complex numbers
+ */
+object ComplexField : Field<Complex> {
+    override val zero: Complex = Complex(0.0, 0.0)
+
+    override fun add(a: Complex, b: Complex): Complex = Complex(a.re + b.re, a.im + b.im)
+
+    override fun multiply(a: Complex, k: Double): Complex = Complex(a.re * k, a.im * k)
+
+    override val one: Complex = Complex(1.0, 0.0)
+
+    override fun multiply(a: Complex, b: Complex): Complex = Complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re)
+
+    override fun divide(a: Complex, b: Complex): Complex = Complex(a.re * b.re + a.im * b.im, a.re * b.im - a.im * b.re) / b.square
+
+}
+
+/**
+ * Complex number class
+ */
+data class Complex(val re: Double, val im: Double) : FieldElement<Complex, ComplexField> {
+    override val self: Complex get() = this
+    override val context: ComplexField
+        get() = ComplexField
+
+    /**
+     * A complex conjugate
+     */
+    val conjugate: Complex
+        get() = Complex(re, -im)
+
+    val square: Double
+        get() = re * re + im * im
+
+    val abs: Double
+        get() = kotlin.math.sqrt(square)
+
+}

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

@@ -1,7 +1,6 @@
 package scientifik.kmath.operations
 
 import kotlin.math.pow
-import kotlin.math.sqrt
 
 /**
  * Field for real values
@@ -29,7 +28,7 @@ object RealField : Field<Real>, TrigonometricOperations<Real>, PowerOperations<R
  *
  * TODO could be replaced by inline class in kotlin 1.3 if it would allow to avoid boxing
  */
-class Real(val value: Double) : Number(), FieldElement<Real, RealField> {
+data class Real(val value: Double) : Number(), FieldElement<Real, RealField> {
     /*
      * The class uses composition instead of inheritance since Double is final
      */
@@ -51,53 +50,22 @@ class Real(val value: Double) : Number(), FieldElement<Real, RealField> {
 
 }
 
-/**
- * A field for complex numbers
- */
-object ComplexField : Field<Complex> {
-    override val zero: Complex = Complex(0.0, 0.0)
-
-    override fun add(a: Complex, b: Complex): Complex = Complex(a.re + b.re, a.im + b.im)
-
-    override fun multiply(a: Complex, k: Double): Complex = Complex(a.re * k, a.im * k)
-
-    override val one: Complex = Complex(1.0, 0.0)
-
-    override fun multiply(a: Complex, b: Complex): Complex = Complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re)
-
-    override fun divide(a: Complex, b: Complex): Complex = Complex(a.re * b.re + a.im * b.im, a.re * b.im - a.im * b.re) / b.square
-
-}
-
-/**
- * Complex number class
- */
-data class Complex(val re: Double, val im: Double) : FieldElement<Complex, ComplexField> {
-    override val self: Complex get() = this
-    override val context: ComplexField
-        get() = ComplexField
-
-    /**
-     * A complex conjugate
-     */
-    val conjugate: Complex
-        get() = Complex(re, -im)
-
-    val square: Double
-        get() = re * re + im * im
-
-    val module: Double
-        get() = sqrt(square)
-
-}
-
 /**
  * A field for double without boxing. Does not produce appropriate field element
  */
-object DoubleField : Field<Double> {
+object DoubleField : Field<Double>, TrigonometricOperations<Double>, PowerOperations<Double>, ExponentialOperations<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 cos(arg: Double): Double = kotlin.math.cos(arg)
+
+    override fun power(arg: Double, pow: Double): Double  = arg.pow(pow)
+
+    override fun exp(arg: Double): Double  =kotlin.math.exp(arg)
+
+    override fun ln(arg: Double): Double  = kotlin.math.ln(arg)
 }

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

@@ -10,7 +10,7 @@ package scientifik.kmath.operations
  * It also allows to override behavior for optional operations
  *
  */
-interface TrigonometricOperations<T>: Field<T> {
+interface TrigonometricOperations<T> : Field<T> {
     fun sin(arg: T): T
     fun cos(arg: T): T
 
@@ -39,10 +39,10 @@ fun <T : MathElement<T, out PowerOperations<T>>> sqr(arg: T): T = arg pow 2.0
 
 /* Exponential */
 
-interface ExponentialOperations<T>{
+interface ExponentialOperations<T> {
     fun exp(arg: T): T
     fun ln(arg: T): 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)
+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)

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

@@ -0,0 +1,112 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.Field
+
+
+/**
+ * A generic buffer for both primitives and objects
+ */
+interface Buffer<T> {
+    operator fun get(index: Int): T
+    operator fun set(index: Int, value: T)
+
+    /**
+     * A shallow copy of the buffer
+     */
+    fun copy(): Buffer<T>
+}
+
+/**
+ * Generic implementation of NDField based on continuous buffer
+ */
+abstract class BufferNDField<T>(shape: List<Int>, field: Field<T>) : NDField<T>(shape, field) {
+
+    /**
+     * Strides for memory access
+     */
+    private val strides: List<Int> by lazy {
+        ArrayList<Int>(shape.size).apply {
+            var current = 1
+            add(1)
+            shape.forEach {
+                current *= it
+                add(current)
+            }
+        }
+    }
+
+    protected fun offset(index: List<Int>): Int {
+        return index.mapIndexed { i, value ->
+            if (value < 0 || value >= shape[i]) {
+                throw RuntimeException("Index out of shape bounds: ($i,$value)")
+            }
+            value * strides[i]
+        }.sum()
+    }
+
+    //TODO introduce a fast way to calculate index of the next element?
+    protected fun index(offset: Int): List<Int> {
+        return sequence {
+            var current = offset
+            var strideIndex = strides.size - 2
+            while (strideIndex >= 0) {
+                yield(current / strides[strideIndex])
+                current %= strides[strideIndex]
+                strideIndex--
+            }
+        }.toList().reversed()
+    }
+
+    private val capacity: Int
+        get() = strides[shape.size]
+
+
+    protected abstract fun createBuffer(capacity: Int, initializer: (Int) -> T): Buffer<T>
+
+    override fun produce(initializer: (List<Int>) -> T): NDArray<T> {
+        val buffer = createBuffer(capacity) { initializer(index(it)) }
+        return BufferNDArray(this, buffer)
+    }
+
+    /**
+     * Produce mutable NDArray instance
+     */
+    fun produceMutable(initializer: (List<Int>) -> T): MutableNDArray<T> {
+        val buffer = createBuffer(capacity) { initializer(index(it)) }
+        return MutableBufferedNDArray(this, buffer)
+    }
+
+
+    private open class BufferNDArray<T>(override val context: BufferNDField<T>, val data: Buffer<T>) : NDArray<T> {
+
+        override fun get(vararg index: Int): T {
+            return data[context.offset(index.asList())]
+        }
+
+        override fun equals(other: Any?): Boolean {
+            if (this === other) return true
+            if (other !is BufferNDArray<*>) return false
+
+            if (context != other.context) return false
+            if (data != other.data) return false
+
+            return true
+        }
+
+        override fun hashCode(): Int {
+            var result = context.hashCode()
+            result = 31 * result + data.hashCode()
+            return result
+        }
+
+        override val self: NDArray<T> get() = this
+    }
+
+    private class MutableBufferedNDArray<T>(context: BufferNDField<T>, data: Buffer<T>): BufferNDArray<T>(context,data), MutableNDArray<T>{
+        override operator fun set(index: List<Int>, value: T){
+            data[context.offset(index)] = value
+        }
+    }
+}
+
+

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

@@ -73,7 +73,9 @@ abstract class NDField<T>(val shape: List<Int>, val field: Field<T>) : Field<NDA
     }
 }
 
-
+/**
+ * Many-dimensional array
+ */
 interface NDArray<T> : FieldElement<NDArray<T>, NDField<T>> {
 
     /**
@@ -125,6 +127,22 @@ interface NDArray<T> : FieldElement<NDArray<T>, NDField<T>> {
     }
 }
 
+/**
+ * In-place mutable [NDArray]
+ */
+interface MutableNDArray<T> : NDArray<T> {
+    operator fun set(index: List<Int>, value: T)
+}
+
+/**
+ * In-place transformation for [MutableNDArray], using given transformation for each element
+ */
+fun <T> MutableNDArray<T>.transformInPlace(action: (List<Int>, T) -> T) {
+    for ((index, oldValue) in this) {
+        this[index] = action(index, oldValue)
+    }
+}
+
 /**
  * Element by element application of any operation on elements to the whole array. Just like in numpy
  */
@@ -136,7 +154,7 @@ operator fun <T> Function1<T, T>.invoke(ndArray: NDArray<T>): NDArray<T> = ndArr
  * Summation operation for [NDArray] and single element
  */
 operator fun <T> NDArray<T>.plus(arg: T): NDArray<T> = transform { _, value ->
-    with(context.field){
+    with(context.field) {
         arg + value
     }
 }
@@ -149,8 +167,8 @@ operator fun <T> T.plus(arg: NDArray<T>): NDArray<T> = arg + this
 /**
  * Subtraction operation between [NDArray] and single element
  */
-operator fun <T> NDArray<T>.minus(arg: T): NDArray<T> = transform  { _, value ->
-    with(context.field){
+operator fun <T> NDArray<T>.minus(arg: T): NDArray<T> = transform { _, value ->
+    with(context.field) {
         arg - value
     }
 }
@@ -159,7 +177,7 @@ operator fun <T> NDArray<T>.minus(arg: T): NDArray<T> = transform  { _, value ->
  * Reverse minus operation
  */
 operator fun <T> T.minus(arg: NDArray<T>): NDArray<T> = arg.transform { _, value ->
-    with(arg.context.field){
+    with(arg.context.field) {
         this@minus - value
     }
 }
@@ -170,7 +188,7 @@ operator fun <T> T.minus(arg: NDArray<T>): NDArray<T> = arg.transform { _, value
  * Product operation for [NDArray] and single element
  */
 operator fun <T> NDArray<T>.times(arg: T): NDArray<T> = transform { _, value ->
-    with(context.field){
+    with(context.field) {
         arg * value
     }
 }
@@ -183,8 +201,8 @@ operator fun <T> T.times(arg: NDArray<T>): NDArray<T> = arg * this
 /**
  * Division operation between [NDArray] and single element
  */
-operator fun <T> NDArray<T>.div(arg: T): NDArray<T> = transform  { _, value ->
-    with(context.field){
+operator fun <T> NDArray<T>.div(arg: T): NDArray<T> = transform { _, value ->
+    with(context.field) {
         arg / value
     }
 }
@@ -193,20 +211,8 @@ operator fun <T> NDArray<T>.div(arg: T): NDArray<T> = transform  { _, value ->
  * Reverse division operation
  */
 operator fun <T> T.div(arg: NDArray<T>): NDArray<T> = arg.transform { _, value ->
-    with(arg.context.field){
-        this@div/ value
+    with(arg.context.field) {
+        this@div / value
     }
 }
 
-/**
- * Create a platform-specific NDArray of doubles
- */
-expect fun realNDArray(shape: List<Int>, initializer: (List<Int>) -> Double = { 0.0 }): NDArray<Double>
-
-fun real2DArray(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }): NDArray<Double> {
-    return realNDArray(listOf(dim1, dim2)) { initializer(it[0], it[1]) }
-}
-
-fun real3DArray(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }): NDArray<Double> {
-    return realNDArray(listOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
-}

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

@@ -0,0 +1,72 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.Field
+
+typealias NDFieldFactory<T> = (shape: List<Int>) -> NDField<T>
+
+/**
+ * The factory class for fast platform-dependent implementation of NDField of doubles
+ */
+expect val realNDFieldFactory: NDFieldFactory<Double>
+
+
+class SimpleNDField<T : Any>(field: Field<T>, shape: List<Int>) : BufferNDField<T>(shape, field) {
+    override fun createBuffer(capacity: Int, initializer: (Int) -> T): Buffer<T> {
+        val array = ArrayList<T>(capacity)
+        (0 until capacity).forEach {
+            array.add(initializer(it))
+        }
+
+        return BufferOfObjects(array)
+    }
+
+    private class BufferOfObjects<T>(val array: ArrayList<T>) : Buffer<T> {
+        override fun get(index: Int): T = array[index]
+
+        override fun set(index: Int, value: T) {
+            array[index] = value
+        }
+
+        override fun copy(): Buffer<T> = BufferOfObjects(ArrayList(array))
+    }
+}
+
+object NDArrays {
+    /**
+     * Create a platform-optimized NDArray of doubles
+     */
+    fun realNDArray(shape: List<Int>, initializer: (List<Int>) -> Double = { 0.0 }): NDArray<Double> {
+        return realNDFieldFactory(shape).produce(initializer)
+    }
+
+    fun real1DArray(dim: Int, initializer: (Int) -> Double = { _ -> 0.0 }): NDArray<Double> {
+        return realNDArray(listOf(dim)) { initializer(it[0]) }
+    }
+
+    fun real2DArray(dim1: Int, dim2: Int, initializer: (Int, Int) -> Double = { _, _ -> 0.0 }): NDArray<Double> {
+        return realNDArray(listOf(dim1, dim2)) { initializer(it[0], it[1]) }
+    }
+
+    fun real3DArray(dim1: Int, dim2: Int, dim3: Int, initializer: (Int, Int, Int) -> Double = { _, _, _ -> 0.0 }): NDArray<Double> {
+        return realNDArray(listOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
+    }
+
+    /**
+     * Simple boxing NDField
+     */
+    fun <T : Any> createFactory(field: Field<T>): NDFieldFactory<T> = { shape -> SimpleNDField(field, shape) }
+
+    /**
+     * Simple boxing NDArray
+     */
+    fun <T : Any> create(field: Field<T>, shape: List<Int>, initializer: (List<Int>) -> T): NDArray<T> {
+        return SimpleNDField(field, shape).produce { initializer(it) }
+    }
+
+    /**
+     * Mutable boxing NDArray
+     */
+    fun <T : Any> createMutable(field: Field<T>, shape: List<Int>, initializer: (List<Int>) -> T): MutableNDArray<T> {
+        return SimpleNDField(field, shape).produceMutable { initializer(it) }
+    }
+}

kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt

@@ -0,0 +1,42 @@
+package scientifik.kmath.expressions
+
+import scientifik.kmath.operations.Complex
+import scientifik.kmath.operations.ComplexField
+import scientifik.kmath.operations.DoubleField
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class FieldExpressionContextTest {
+    @Test
+    fun testExpression() {
+        val context = FieldExpressionContext(DoubleField)
+        val expression = with(context) {
+            val x = variable("x", 2.0)
+            x * x + 2 * x + 1.0
+        }
+        assertEquals(expression("x" to 1.0), 4.0)
+        assertEquals(expression(), 9.0)
+    }
+
+    @Test
+    fun testComplex() {
+        val context = FieldExpressionContext(ComplexField)
+        val expression = with(context) {
+            val x = variable("x", Complex(2.0, 0.0))
+            x * x + 2 * x + 1.0
+        }
+        assertEquals(expression("x" to Complex(1.0, 0.0)), Complex(4.0, 0.0))
+        assertEquals(expression(), Complex(9.0, 0.0))
+    }
+
+    @Test
+    fun separateContext() {
+        fun <T> FieldExpressionContext<T>.expression(): Expression<T>{
+            val x = variable("x")
+            return x * x + 2 * x + 1.0
+        }
+
+        val expression = FieldExpressionContext(DoubleField).expression()
+        assertEquals(expression("x" to 1.0), 4.0)
+    }
+}

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

@@ -0,0 +1,34 @@
+package scientifik.kmath.linear
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class ArrayMatrixTest {
+
+    @Test
+    fun testSum() {
+        val vector1 = realVector(5) { it.toDouble() }
+        val vector2 = realVector(5) { 5 - it.toDouble() }
+        val sum = vector1 + vector2
+        assertEquals(5.0, sum[2])
+    }
+
+    @Test
+    fun testVectorToMatrix() {
+        val vector = realVector(5) { it.toDouble() }
+        val matrix = vector.toMatrix()
+        assertEquals(4.0, matrix[4, 0])
+    }
+
+
+    @Test
+    fun testDot() {
+        val vector1 = realVector(5) { it.toDouble() }
+        val vector2 = realVector(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/RealLUSolverTest.kt

@@ -0,0 +1,21 @@
+package scientifik.kmath.linear
+
+import scientifik.kmath.operations.DoubleField
+import kotlin.test.Test
+import kotlin.test.assertTrue
+
+class RealLUSolverTest {
+    @Test
+    fun testInvertOne() {
+        val matrix = Matrix.diagonal(2, 2, DoubleField)
+        val inverted = RealLUSolver.inverse(matrix)
+        assertTrue { Matrix.equals(matrix,inverted) }
+    }
+
+//    @Test
+//    fun testInvert() {
+//        val matrix = realMatrix(2,2){}
+//        val inverted = RealLUSolver.inverse(matrix)
+//        assertTrue { Matrix.equals(matrix,inverted) }
+//    }
+}

kmath-core/src/commonTest/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt

@@ -0,0 +1,13 @@
+package scientifik.kmath.misc
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class CumulativeKtTest {
+    @Test
+    fun testCumulativeSum() {
+        val initial = listOf(-1.0, 2.0, 1.0, 1.0)
+        val cumulative = initial.cumulativeSum()
+        assertEquals(listOf(-1.0, 1.0, 2.0, 3.0), cumulative)
+    }
+}

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

@@ -1,8 +1,9 @@
 package scientifik.kmath.structures
 
-import org.junit.Assert.assertEquals
+import scientifik.kmath.structures.NDArrays.real2DArray
 import kotlin.math.pow
 import kotlin.test.Test
+import kotlin.test.assertEquals
 
 class RealNDFieldTest {
     val array1 = real2DArray(3, 3) { i, j -> (i + j).toDouble() }
@@ -11,13 +12,13 @@ class RealNDFieldTest {
     @Test
     fun testSum() {
         val sum = array1 + array2
-        assertEquals(4.0, sum[2, 2], 0.1)
+        assertEquals(4.0, sum[2, 2])
     }
 
     @Test
     fun testProduct() {
         val product = array1 * array2
-        assertEquals(0.0, product[2, 2], 0.1)
+        assertEquals(0.0, product[2, 2])
     }
 
     @Test
@@ -28,7 +29,7 @@ class RealNDFieldTest {
         for (i in 0..2) {
             for (j in 0..2) {
                 val expected = (i * 10 + j).toDouble()
-                assertEquals("Error at index [$i, $j]", expected, array[i, j], 0.1)
+                assertEquals(expected, array[i, j],"Error at index [$i, $j]")
             }
         }
     }
@@ -37,6 +38,6 @@ class RealNDFieldTest {
     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],0.01)
+        assertEquals(10.0, result[1,1])
     }
 }

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

@@ -0,0 +1,16 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.DoubleField
+import scientifik.kmath.structures.NDArrays.create
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+
+class SimpleNDFieldTest{
+    @Test
+    fun testStrides(){
+        val ndArray = create(DoubleField, listOf(10,10)){(it[0]+it[1]).toDouble()}
+        assertEquals(ndArray[5,5], 10.0)
+    }
+
+}

kmath-core/src/jmh/kotlin/scietifik/kmath/structures/ArrayBenchmark.kt

@@ -0,0 +1,53 @@
+package scietifik.kmath.structures
+
+import org.openjdk.jmh.annotations.*
+import java.nio.IntBuffer
+
+
+@Fork(1)
+@Warmup(iterations = 2)
+@Measurement(iterations = 50)
+@State(Scope.Benchmark)
+open class ArrayBenchmark {
+
+    lateinit var array: IntArray
+    lateinit var arrayBuffer: IntBuffer
+    lateinit var nativeBuffer: IntBuffer
+
+    @Setup
+    fun setup() {
+        array = IntArray(10000) { it }
+        arrayBuffer = IntBuffer.wrap(array)
+        nativeBuffer = IntBuffer.allocate(10000)
+        for (i in 0 until 10000) {
+            nativeBuffer.put(i,i)
+        }
+    }
+
+    @Benchmark
+    fun benchmarkArrayRead() {
+        var res = 0
+        for (i in 1..10000) {
+            res += array[10000 - i]
+        }
+        print(res)
+    }
+
+    @Benchmark
+    fun benchmarkBufferRead() {
+        var res = 0
+        for (i in 1..10000) {
+            res += arrayBuffer.get(10000 - i)
+        }
+        print(res)
+    }
+
+    @Benchmark
+    fun nativeBufferRead() {
+        var res = 0
+        for (i in 1..10000) {
+            res += nativeBuffer.get(10000 - i)
+        }
+        print(res)
+    }
+}

kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt

@@ -0,0 +1,8 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.DoubleField
+
+/**
+ * Using boxing implementation for js
+ */
+actual val realNDFieldFactory: NDFieldFactory<Double> = NDArrays.createFactory(DoubleField)

kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt

@@ -0,0 +1,26 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.DoubleField
+import java.nio.DoubleBuffer
+
+private class RealNDField(shape: List<Int>) : BufferNDField<Double>(shape, DoubleField) {
+    override fun createBuffer(capacity: Int, initializer: (Int) -> Double): Buffer<Double> {
+        val array = DoubleArray(capacity, initializer)
+        val buffer = DoubleBuffer.wrap(array)
+        return BufferOfDoubles(buffer)
+    }
+
+    private class BufferOfDoubles(val buffer: DoubleBuffer): Buffer<Double>{
+        override fun get(index: Int): Double = buffer.get(index)
+
+        override fun set(index: Int, value: Double) {
+            buffer.put(index, value)
+        }
+
+        override fun copy(): Buffer<Double> {
+            return BufferOfDoubles(buffer)
+        }
+    }
+}
+
+actual val realNDFieldFactory: NDFieldFactory<Double> = { shape -> RealNDField(shape) }

kmath-jvm/build.gradle

@@ -1,21 +0,0 @@
-apply plugin: 'kotlin-platform-jvm'
-
-repositories {
-    mavenCentral()
-}
-
-dependencies {
-    expectedBy project(":kmath-common")
-    compile "org.jetbrains.kotlin:kotlin-stdlib-jdk8:$kotlin_version"
-    testCompile "junit:junit:4.12"
-    testCompile "org.jetbrains.kotlin:kotlin-test-junit:$kotlin_version"
-    testCompile "org.jetbrains.kotlin:kotlin-test:$kotlin_version"
-}
-
-compileKotlin {
-    kotlinOptions.jvmTarget = "1.8"
-}
-compileTestKotlin {
-    kotlinOptions.jvmTarget = "1.8"
-}
-sourceCompatibility = "1.8"

kmath-jvm/src/main/kotlin/scientifik/kmath/structures/RealNDArray.kt

@@ -1,80 +0,0 @@
-package scientifik.kmath.structures
-
-import scientifik.kmath.operations.DoubleField
-import java.nio.DoubleBuffer
-
-private class RealNDField(shape: List<Int>) : NDField<Double>(shape, DoubleField) {
-
-    /**
-     * Strides for memory access
-     */
-    private val strides: List<Int> by lazy {
-        ArrayList<Int>(shape.size).apply {
-            var current = 1
-            add(1)
-            shape.forEach {
-                current *= it
-                add(current)
-            }
-        }
-    }
-
-    fun offset(index: List<Int>): Int {
-        return index.mapIndexed { i, value ->
-            if (value < 0 || value >= shape[i]) {
-                throw RuntimeException("Index out of shape bounds: ($i,$value)")
-            }
-            value * strides[i]
-        }.sum()
-    }
-
-    val capacity: Int
-        get() = strides[shape.size]
-
-
-    override fun produce(initializer: (List<Int>) -> Double): NDArray<Double> {
-        //TODO use sparse arrays for large capacities
-        val buffer = DoubleBuffer.allocate(capacity)
-        //FIXME there could be performance degradation due to iteration procedure. Replace by straight iteration
-        NDArray.iterateIndexes(shape).forEach {
-            buffer.put(offset(it), initializer(it))
-        }
-        return RealNDArray(this, buffer)
-    }
-
-    class RealNDArray(override val context: RealNDField, val data: DoubleBuffer) : NDArray<Double> {
-
-        override fun get(vararg index: Int): Double {
-            return data.get(context.offset(index.asList()))
-        }
-
-        override fun equals(other: Any?): Boolean {
-            if (this === other) return true
-            if (javaClass != other?.javaClass) return false
-
-            other as RealNDArray
-
-            if (context.shape != other.context.shape) return false
-            if (data != other.data) return false
-
-            return true
-        }
-
-        override fun hashCode(): Int {
-            var result = context.shape.hashCode()
-            result = 31 * result + data.hashCode()
-            return result
-        }
-
-        //TODO generate fixed hash code for quick comparison?
-
-
-        override val self: NDArray<Double> get() = this
-    }
-}
-
-
-actual fun realNDArray(shape: List<Int>, initializer: (List<Int>) -> Double): NDArray<Double> {
-    //TODO cache fields?
-    return RealNDField(shape).produce { initializer(it) }
-}

kmath-jvm/src/test/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt

@@ -1,26 +0,0 @@
-package scientifik.kmath.structures
-
-import org.junit.Assert.assertEquals
-import org.junit.Test
-
-class ArrayMatrixTest {
-
-    @Test
-    fun testSum() {
-        val vector1 = realVector(5) { it.toDouble() }
-        val vector2 = realVector(5) { 5 - it.toDouble() }
-        val sum = vector1 + vector2
-        assertEquals(5.0, sum[2, 0], 0.1)
-    }
-
-    @Test
-    fun testDot() {
-        val vector1 = realVector(5) { it.toDouble() }
-        val vector2 = realVector(5) { 5 - it.toDouble() }
-        val product = with(vector1.context) {
-            vector1 dot (vector2.transpose())
-        }
-
-        assertEquals(10.0, product[1, 0], 0.1)
-    }
-}

settings.gradle

@@ -1,4 +1,13 @@
-rootProject.name = 'kmath'
-include 'kmath-common'
-include 'kmath-jvm'
+pluginManagement {
+    repositories {
+        maven { url = 'http://dl.bintray.com/kotlin/kotlin-eap' }
+        mavenCentral()
+        maven { url = 'https://plugins.gradle.org/m2/' }
+    }
+}
+
+enableFeaturePreview('GRADLE_METADATA')
+
+rootProject.name = 'kmath'
+include ':kmath-core'