From 982bee02ba89700eee0d9a8da0cac8ddbb752b4e Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Tue, 21 Aug 2018 21:27:46 +0300
Subject: [PATCH 01/17] Basic cumulative operation and cumulative sum #10

---
 build.gradle                                  |  2 +-
 .../scientifik/kmath/misc/Cumulative.kt       | 59 +++++++++++++++++++
 .../scientifik.kmath.misc/CumulativeKtTest.kt | 13 ++++
 3 files changed, 73 insertions(+), 1 deletion(-)
 create mode 100644 kmath-common/src/main/kotlin/scientifik/kmath/misc/Cumulative.kt
 create mode 100644 kmath-common/src/test/kotlin/scientifik.kmath.misc/CumulativeKtTest.kt

diff --git a/build.gradle b/build.gradle
index bce9c1492..1d08c5527 100644
--- a/build.gradle
+++ b/build.gradle
@@ -1,5 +1,5 @@
 buildscript {
-    ext.kotlin_version = '1.2.60'
+    ext.kotlin_version = '1.2.61'
 
     repositories {
         mavenCentral()
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/misc/Cumulative.kt b/kmath-common/src/main/kotlin/scientifik/kmath/misc/Cumulative.kt
new file mode 100644
index 000000000..4d4f8ced6
--- /dev/null
+++ b/kmath-common/src/main/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}
\ No newline at end of file
diff --git a/kmath-common/src/test/kotlin/scientifik.kmath.misc/CumulativeKtTest.kt b/kmath-common/src/test/kotlin/scientifik.kmath.misc/CumulativeKtTest.kt
new file mode 100644
index 000000000..e7c99e7d0
--- /dev/null
+++ b/kmath-common/src/test/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)
+    }
+}
\ No newline at end of file
-- 
2.34.1


From 544f65d9a202f2637346ee3082b8ddba0b6eda53 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Wed, 22 Aug 2018 12:45:05 +0300
Subject: [PATCH 02/17] added maven plugin

---
 build.gradle | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/build.gradle b/build.gradle
index 1d08c5527..a57ba6cc8 100644
--- a/build.gradle
+++ b/build.gradle
@@ -9,6 +9,10 @@ buildscript {
     }
 }
 
+subprojects{
+    apply plugin: 'maven'
+}
+
 group = 'scientifik'
 version = '0.1 - SNAPSHOT'
 
-- 
2.34.1


From ce5c281427330793b61b52ec6bed9ff4f5d78d0c Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Fri, 24 Aug 2018 20:55:44 +0300
Subject: [PATCH 03/17] Fixed version syntax

---
 build.gradle                                              | 2 +-
 .../scientifik/kmath/operations/OptionalOperations.kt     | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/build.gradle b/build.gradle
index a57ba6cc8..dbe2ef5de 100644
--- a/build.gradle
+++ b/build.gradle
@@ -14,5 +14,5 @@ subprojects{
 }
 
 group = 'scientifik'
-version = '0.1 - SNAPSHOT'
+version = '0.0.1-SNAPSHOT'
 
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/operations/OptionalOperations.kt b/kmath-common/src/main/kotlin/scientifik/kmath/operations/OptionalOperations.kt
index 18ef89407..7628cb4fc 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/operations/OptionalOperations.kt
+++ b/kmath-common/src/main/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)
\ No newline at end of file
+fun <T : MathElement<T, out ExponentialOperations<T>>> exp(arg: T): T = arg.context.exp(arg)
+fun <T : MathElement<T, out ExponentialOperations<T>>> ln(arg: T): T = arg.context.ln(arg)
\ No newline at end of file
-- 
2.34.1


From 8da54786a0495ece41fa89c3ea9d3a36ff3525a0 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Fri, 24 Aug 2018 21:11:59 +0300
Subject: [PATCH 04/17] Minor documentation update

---
 README.md                                     | 20 ++++++++++++++++---
 .../scientifik/kmath/operations/Algebra.kt    |  2 ++
 2 files changed, 19 insertions(+), 3 deletions(-)

diff --git a/README.md b/README.md
index 13db36602..ca7d1d09a 100644
--- a/README.md
+++ b/README.md
@@ -6,15 +6,29 @@ and `scipy` it is modular and has a lightweight core.
 
 * **Algebra**
     * Mathematical operation entities like rings, spaces and fields with (**TODO** add example to wiki) 
-    * Basic linear algebra operations (summs products, etc) backed by `Space` API.
+    * Basic linear algebra operations (sums products, etc) backed by `Space` API.
     * [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).
 
+## 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.
+
+* **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.
+
+* **Messaging** A mathematical notation to support multilanguage and multinod 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.
+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 
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt b/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt
index 6f06f5980..985cf3b14 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt
+++ b/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt
@@ -3,6 +3,8 @@ 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>{
     /**
-- 
2.34.1


From c344380120123c91aead15c4014512dfc84bfaf5 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sat, 25 Aug 2018 20:58:26 +0300
Subject: [PATCH 05/17] Expressions!

---
 README.md                                     |  8 +--
 .../kmath/expressions/Expression.kt           | 62 +++++++++++++++++++
 .../scientifik/kmath/operations/Fields.kt     | 13 +++-
 .../expressions/FieldExpressionContextTest.kt | 18 ++++++
 .../kmath/misc}/CumulativeKtTest.kt           |  0
 5 files changed, 95 insertions(+), 6 deletions(-)
 create mode 100644 kmath-common/src/main/kotlin/scientifik/kmath/expressions/Expression.kt
 create mode 100644 kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
 rename kmath-common/src/test/kotlin/{scientifik.kmath.misc => scientifik/kmath/misc}/CumulativeKtTest.kt (100%)

diff --git a/README.md b/README.md
index ca7d1d09a..7defe0482 100644
--- a/README.md
+++ b/README.md
@@ -11,16 +11,16 @@ and `scipy` it is modular and has a lightweight core.
 * **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.
 
-* **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.
-
 * **Messaging** A mathematical notation to support multilanguage and multinod communication for mathematical tasks.
 
 ## Multi-platform support
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/expressions/Expression.kt b/kmath-common/src/main/kotlin/scientifik/kmath/expressions/Expression.kt
new file mode 100644
index 000000000..0a34b536c
--- /dev/null
+++ b/kmath-common/src/main/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)
+}
\ No newline at end of file
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/operations/Fields.kt b/kmath-common/src/main/kotlin/scientifik/kmath/operations/Fields.kt
index 2f3eb99ff..0d2600d82 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/operations/Fields.kt
+++ b/kmath-common/src/main/kotlin/scientifik/kmath/operations/Fields.kt
@@ -29,7 +29,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
      */
@@ -94,10 +94,19 @@ data class Complex(val re: Double, val im: Double) : FieldElement<Complex, Compl
 /**
  * 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)
 }
\ No newline at end of file
diff --git a/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt b/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
new file mode 100644
index 000000000..5e9832604
--- /dev/null
+++ b/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
@@ -0,0 +1,18 @@
+package scientifik.kmath.expressions
+
+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 * one
+        }
+        assertEquals(expression("x" to 1.0), 4.0)
+        assertEquals(expression(), 9.0)
+    }
+}
\ No newline at end of file
diff --git a/kmath-common/src/test/kotlin/scientifik.kmath.misc/CumulativeKtTest.kt b/kmath-common/src/test/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt
similarity index 100%
rename from kmath-common/src/test/kotlin/scientifik.kmath.misc/CumulativeKtTest.kt
rename to kmath-common/src/test/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt
-- 
2.34.1


From 3576ddaf317418be793f8516afd2f3d67a408615 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sat, 25 Aug 2018 21:15:55 +0300
Subject: [PATCH 06/17] Added extension to field to use numbers as first class
 citizens

---
 .../src/main/kotlin/scientifik/kmath/operations/Algebra.kt | 7 +++++--
 .../kmath/expressions/FieldExpressionContextTest.kt        | 2 +-
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt b/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt
index 985cf3b14..4613dda7c 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt
+++ b/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt
@@ -6,7 +6,7 @@ package scientifik.kmath.operations
  * @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> {
     /**
      * The context this element belongs to
      */
@@ -53,7 +53,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> {
+interface SpaceElement<T, S : Space<T>> : MathElement<T, S> {
 
     /**
      * Self value. Needed for static type checking. Needed to avoid type erasure on JVM.
@@ -101,6 +101,9 @@ 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
 }
 
 /**
diff --git a/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt b/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
index 5e9832604..3d68a4674 100644
--- a/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
+++ b/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
@@ -10,7 +10,7 @@ class FieldExpressionContextTest {
         val context = FieldExpressionContext(DoubleField)
         val expression = with(context) {
             val x = variable("x", 2.0)
-            x * x + 2 * x + 1.0 * one
+            x * x + 2 * x + 1.0
         }
         assertEquals(expression("x" to 1.0), 4.0)
         assertEquals(expression(), 9.0)
-- 
2.34.1


From 520765abb9d296401d47570a96cb67fa1ae4ac8c Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sat, 25 Aug 2018 21:18:17 +0300
Subject: [PATCH 07/17] The same for minus

---
 .../src/main/kotlin/scientifik/kmath/operations/Algebra.kt     | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt b/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt
index 4613dda7c..9938f4856 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt
+++ b/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt
@@ -104,6 +104,9 @@ interface Field<T> : Ring<T> {
 
     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
 }
 
 /**
-- 
2.34.1


From 6df416f064828156b7323e907e1490786ed3c362 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sat, 25 Aug 2018 21:22:07 +0300
Subject: [PATCH 08/17] Added ComplexField test for expressions.

---
 .../kmath/expressions/FieldExpressionContextTest.kt | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt b/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
index 3d68a4674..dda267ec8 100644
--- a/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
+++ b/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
@@ -1,5 +1,7 @@
 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
@@ -15,4 +17,15 @@ class FieldExpressionContextTest {
         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))
+    }
 }
\ No newline at end of file
-- 
2.34.1


From 6c25042f0f4feaf2bf6fbb358be1e0e08a020fcf Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sat, 25 Aug 2018 21:36:50 +0300
Subject: [PATCH 09/17] Context separation test for expressions

---
 .../kmath/expressions/FieldExpressionContextTest.kt   | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt b/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
index dda267ec8..543e13d79 100644
--- a/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
+++ b/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
@@ -28,4 +28,15 @@ class FieldExpressionContextTest {
         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)
+    }
 }
\ No newline at end of file
-- 
2.34.1


From ba63b2e373d8b8199f97a523d03123e7563d8f6f Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Wed, 29 Aug 2018 13:45:06 +0300
Subject: [PATCH 10/17] Generic definition for NDArray

---
 kmath-common/build.gradle                     | 12 ++-
 .../kmath/structures/BufferNDField.kt         | 93 +++++++++++++++++++
 .../scientifik/kmath/structures/NDArray.kt    | 12 ---
 .../kmath/structures/NDArrayFactories.kt      | 38 ++++++++
 .../kmath/structures/SimpleNDFieldTest.kt     | 15 +++
 kmath-jvm/build.gradle                        |  5 +-
 .../kmath/structures/ArrayBenchmark.kt        | 53 +++++++++++
 .../kmath/structures/RealNDArray.kt           | 76 ++-------------
 8 files changed, 223 insertions(+), 81 deletions(-)
 create mode 100644 kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt
 create mode 100644 kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArrayFactories.kt
 create mode 100644 kmath-common/src/test/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
 create mode 100644 kmath-jvm/src/jmh/kotlin/scietifik/kmath/structures/ArrayBenchmark.kt

diff --git a/kmath-common/build.gradle b/kmath-common/build.gradle
index 2712fb0eb..eecc32818 100644
--- a/kmath-common/build.gradle
+++ b/kmath-common/build.gradle
@@ -1,6 +1,8 @@
-description = "Platform-independent interfaces for kotlin maths"
+plugins{
+    id "kotlin-platform-common"
+}
 
-apply plugin: 'kotlin-platform-common'
+description = "Platform-independent interfaces for kotlin maths"
 
 repositories {
     mavenCentral()
@@ -12,3 +14,9 @@ dependencies {
     testCompile "org.jetbrains.kotlin:kotlin-test-common:$kotlin_version"
 }
 
+kotlin {
+    experimental {
+        coroutines "enable"
+    }
+}
+
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt b/kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt
new file mode 100644
index 000000000..6a2df9acc
--- /dev/null
+++ b/kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt
@@ -0,0 +1,93 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.Field
+import kotlin.coroutines.experimental.buildSequence
+
+
+/**
+ * A generic buffer for both primitives and objects
+ */
+interface Buffer<T> {
+    operator fun get(index: Int): T
+    operator fun set(index: Int, value: 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()
+    }
+
+    protected fun index(offset: Int): List<Int>{
+        return buildSequence {
+            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)
+    }
+
+
+    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
+    }
+}
+
+
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArray.kt b/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArray.kt
index 5acac4ff1..462cd6adb 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArray.kt
+++ b/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArray.kt
@@ -198,15 +198,3 @@ operator fun <T> T.div(arg: NDArray<T>): NDArray<T> = arg.transform { _, 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]) }
-}
\ No newline at end of file
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArrayFactories.kt b/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArrayFactories.kt
new file mode 100644
index 000000000..b3e3850e9
--- /dev/null
+++ b/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArrayFactories.kt
@@ -0,0 +1,38 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.Field
+
+/**
+ * Create a platform-optimized 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]) }
+}
+
+
+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 object : Buffer<T> {
+            override fun get(index: Int): T = array[index]
+
+            override fun set(index: Int, value: T) {
+                array[index] = initializer(index)
+            }
+        }
+    }
+}
+
+fun <T: Any> simpleNDArray(field: Field<T>, shape: List<Int>, initializer: (List<Int>) -> T): NDArray<T> {
+    return SimpleNDField(field, shape).produce { initializer(it) }
+}
\ No newline at end of file
diff --git a/kmath-common/src/test/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt b/kmath-common/src/test/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
new file mode 100644
index 000000000..5a4317251
--- /dev/null
+++ b/kmath-common/src/test/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
@@ -0,0 +1,15 @@
+package scientifik.kmath.structures
+
+import scientifik.kmath.operations.DoubleField
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+
+class SimpleNDFieldTest{
+    @Test
+    fun testStrides(){
+        val ndArray = simpleNDArray(DoubleField, listOf(10,10)){(it[0]+it[1]).toDouble()}
+        assertEquals(ndArray[5,5], 10.0)
+    }
+
+}
\ No newline at end of file
diff --git a/kmath-jvm/build.gradle b/kmath-jvm/build.gradle
index 82bda0270..0f1240a3a 100644
--- a/kmath-jvm/build.gradle
+++ b/kmath-jvm/build.gradle
@@ -1,4 +1,7 @@
-apply plugin: 'kotlin-platform-jvm'
+plugins{
+    id "kotlin-platform-jvm"
+    id "me.champeau.gradle.jmh" version "0.4.5"
+}
 
 repositories {
     mavenCentral()
diff --git a/kmath-jvm/src/jmh/kotlin/scietifik/kmath/structures/ArrayBenchmark.kt b/kmath-jvm/src/jmh/kotlin/scietifik/kmath/structures/ArrayBenchmark.kt
new file mode 100644
index 000000000..3430b14d4
--- /dev/null
+++ b/kmath-jvm/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)
+    }
+}
\ No newline at end of file
diff --git a/kmath-jvm/src/main/kotlin/scientifik/kmath/structures/RealNDArray.kt b/kmath-jvm/src/main/kotlin/scientifik/kmath/structures/RealNDArray.kt
index 2c213f91e..8099e8700 100644
--- a/kmath-jvm/src/main/kotlin/scientifik/kmath/structures/RealNDArray.kt
+++ b/kmath-jvm/src/main/kotlin/scientifik/kmath/structures/RealNDArray.kt
@@ -3,78 +3,22 @@ package scientifik.kmath.structures
 import scientifik.kmath.operations.DoubleField
 import java.nio.DoubleBuffer
 
-private class RealNDField(shape: List<Int>) : NDField<Double>(shape, DoubleField) {
+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 object : Buffer<Double> {
+            override fun get(index: Int): Double = buffer.get(index)
 
-    /**
-     * 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)
+            override fun set(index: Int, value: Double) {
+                buffer.put(index, value)
             }
         }
     }
-
-    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?
+    //TODO create a cache for fields to save time generating strides?
+
     return RealNDField(shape).produce { initializer(it) }
 }
\ No newline at end of file
-- 
2.34.1


From 9c064b00b14092a0e3f45a0ac53b255827c6e450 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Tue, 11 Sep 2018 22:00:00 +0300
Subject: [PATCH 11/17] Grid generators for 1D  double array

---
 .../kotlin/scientifik/kmath/misc/Grids.kt     | 40 +++++++++++++++++++
 .../kmath/structures/BufferNDField.kt         |  9 +++--
 kmath-jvm/build.gradle                        |  8 +++-
 3 files changed, 52 insertions(+), 5 deletions(-)
 create mode 100644 kmath-common/src/main/kotlin/scientifik/kmath/misc/Grids.kt

diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/misc/Grids.kt b/kmath-common/src/main/kotlin/scientifik/kmath/misc/Grids.kt
new file mode 100644
index 000000000..5ea4bd949
--- /dev/null
+++ b/kmath-common/src/main/kotlin/scientifik/kmath/misc/Grids.kt
@@ -0,0 +1,40 @@
+package scientifik.kmath.misc
+
+import kotlin.coroutines.experimental.buildSequence
+
+
+/**
+ * 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 -> buildSequence {
+            var current = start
+            while (current <= endInclusive) {
+                yield(current)
+                current += step
+            }
+        }
+        else -> buildSequence {
+            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 }
+}
\ No newline at end of file
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt b/kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt
index 6a2df9acc..5b3ef4c6a 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt
+++ b/kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt
@@ -40,11 +40,12 @@ abstract class BufferNDField<T>(shape: List<Int>, field: Field<T>) : NDField<T>(
         }.sum()
     }
 
-    protected fun index(offset: Int): List<Int>{
+    //TODO introduce a fast way to calculate index of the next element?
+    protected fun index(offset: Int): List<Int> {
         return buildSequence {
             var current = offset
-            var strideIndex = strides.size-2
-            while (strideIndex>=0){
+            var strideIndex = strides.size - 2
+            while (strideIndex >= 0) {
                 yield(current / strides[strideIndex])
                 current %= strides[strideIndex]
                 strideIndex--
@@ -59,7 +60,7 @@ abstract class BufferNDField<T>(shape: List<Int>, field: Field<T>) : NDField<T>(
     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))}
+        val buffer = createBuffer(capacity) { initializer(index(it)) }
         return BufferNDArray(this, buffer)
     }
 
diff --git a/kmath-jvm/build.gradle b/kmath-jvm/build.gradle
index 0f1240a3a..110dc0287 100644
--- a/kmath-jvm/build.gradle
+++ b/kmath-jvm/build.gradle
@@ -21,4 +21,10 @@ compileKotlin {
 compileTestKotlin {
     kotlinOptions.jvmTarget = "1.8"
 }
-sourceCompatibility = "1.8"
\ No newline at end of file
+sourceCompatibility = "1.8"
+
+kotlin {
+    experimental {
+        coroutines "enable"
+    }
+}
\ No newline at end of file
-- 
2.34.1


From 38c7f4382aebf7cd728da1266b5ecc4e6ffa6e33 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Fri, 28 Sep 2018 21:07:42 +0300
Subject: [PATCH 12/17] Artifactory deploy

---
 .gitignore             |  1 +
 README.md              |  2 ++
 build.gradle           | 42 +++++++++++++++++++++++++++++++++++++-----
 kmath-jvm/build.gradle |  8 ++++++++
 4 files changed, 48 insertions(+), 5 deletions(-)

diff --git a/.gitignore b/.gitignore
index 0fa59f100..5b55fc854 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,3 +8,4 @@
 # Cache of project
 .gradletasknamecache
 
+gradle.properties
\ No newline at end of file
diff --git a/README.md b/README.md
index 7defe0482..8f914d791 100644
--- a/README.md
+++ b/README.md
@@ -30,6 +30,8 @@ platform even if they could be done in common module because of platform perform
 Currently the main focus of development is the JVM platform, contribution of implementations for Kotlin - Native and
 Kotlin - JS is welcome. 
 
+It is planned to move project structure to [kotlin gradle mpp structure](https://github.com/h0tk3y/k-new-mpp-samples) after release of kotlin 1.3. 
+
 ## Performance
 The calculation performance is one of major goals of KMath in the future, but in some cases it is not possible to achieve 
 both performance and flexibility. We expect to firstly focus on creating convenient universal API and then work on 
diff --git a/build.gradle b/build.gradle
index dbe2ef5de..a3af48ee0 100644
--- a/build.gradle
+++ b/build.gradle
@@ -1,18 +1,50 @@
 buildscript {
-    ext.kotlin_version = '1.2.61'
+    ext.kotlin_version = '1.2.71'
 
     repositories {
-        mavenCentral()
+        jcenter()
     }
     dependencies {
         classpath "org.jetbrains.kotlin:kotlin-gradle-plugin:$kotlin_version"
+        classpath "org.jfrog.buildinfo:build-info-extractor-gradle:4+"
     }
 }
 
-subprojects{
+allprojects{
     apply plugin: 'maven'
+    apply plugin: 'maven-publish'
+    apply plugin: "com.jfrog.artifactory"
+
+    group = 'scientifik'
+    version = '0.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}"
+            maven = true
+
+        }
+
+        defaults {
+            publications 'defaultPublication'
+            publishBuildInfo = true
+            publishArtifacts = true
+            publishPom = true
+            publishIvy = false
+        }
+    }
+    resolve {
+        repository {
+            repoKey = 'gradle-dev'
+            username = "${artifactory_user}"
+            password = "${artifactory_password}"
+            maven = true
+        }
+    }
+}
\ No newline at end of file
diff --git a/kmath-jvm/build.gradle b/kmath-jvm/build.gradle
index 110dc0287..b8952a379 100644
--- a/kmath-jvm/build.gradle
+++ b/kmath-jvm/build.gradle
@@ -27,4 +27,12 @@ kotlin {
     experimental {
         coroutines "enable"
     }
+}
+
+publishing {
+    publications {
+        defaultPublication(MavenPublication) {
+            from components.java
+        }
+    }
 }
\ No newline at end of file
-- 
2.34.1


From c4b334976acfe4ddb0d0ba713bd653d2ec860d26 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sun, 30 Sep 2018 17:18:04 +0300
Subject: [PATCH 13/17] Switch to kotlin 1.3 eap and mpp build

---
 README.md                                     | 29 ++++++++-
 build.gradle                                  | 17 +++---
 kmath-common/build.gradle                     | 22 -------
 .../kmath/structures/NDArrayFactories.kt      | 38 ------------
 kmath-core/build.gradle                       | 58 ++++++++++++++++++
 .../kmath/expressions/Expression.kt           |  0
 .../scientifik/kmath/misc/Cumulative.kt       |  0
 .../kotlin/scientifik/kmath/misc/Grids.kt     |  7 +--
 .../scientifik/kmath/operations/Algebra.kt    |  0
 .../scientifik/kmath/operations/Fields.kt     |  0
 .../kmath/operations/OptionalOperations.kt    |  0
 .../kmath/structures/BufferNDField.kt         |  3 +-
 .../kmath/structures/LinearAlgrebra.kt        | 61 +++++++++++++------
 .../scientifik/kmath/structures/NDArray.kt    | 18 +++---
 .../scientifik/kmath/structures/NDArrays.kt   | 54 ++++++++++++++++
 .../expressions/FieldExpressionContextTest.kt |  0
 .../scientifik/kmath/misc/CumulativeKtTest.kt |  0
 .../kmath/operations/RealFieldTest.kt         |  0
 .../kmath/structures/ArrayMatrixTest.kt       |  8 +--
 .../kmath/structures/RealNDFieldTest.kt       | 11 ++--
 .../kmath/structures/SimpleNDFieldTest.kt     |  1 +
 .../kmath/structures/ArrayBenchmark.kt        |  0
 .../scientifik/kmath/structures/_NDArrays.kt  |  8 +++
 .../scientifik/kmath/structures/_NDArrays.kt  |  6 +-
 kmath-jvm/build.gradle                        | 38 ------------
 settings.gradle                               | 15 ++++-
 26 files changed, 234 insertions(+), 160 deletions(-)
 delete mode 100644 kmath-common/build.gradle
 delete mode 100644 kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArrayFactories.kt
 create mode 100644 kmath-core/build.gradle
 rename {kmath-common/src/main => kmath-core/src/commonMain}/kotlin/scientifik/kmath/expressions/Expression.kt (100%)
 rename {kmath-common/src/main => kmath-core/src/commonMain}/kotlin/scientifik/kmath/misc/Cumulative.kt (100%)
 rename {kmath-common/src/main => kmath-core/src/commonMain}/kotlin/scientifik/kmath/misc/Grids.kt (90%)
 rename {kmath-common/src/main => kmath-core/src/commonMain}/kotlin/scientifik/kmath/operations/Algebra.kt (100%)
 rename {kmath-common/src/main => kmath-core/src/commonMain}/kotlin/scientifik/kmath/operations/Fields.kt (100%)
 rename {kmath-common/src/main => kmath-core/src/commonMain}/kotlin/scientifik/kmath/operations/OptionalOperations.kt (100%)
 rename {kmath-common/src/main => kmath-core/src/commonMain}/kotlin/scientifik/kmath/structures/BufferNDField.kt (96%)
 rename {kmath-common/src/main => kmath-core/src/commonMain}/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt (66%)
 rename {kmath-common/src/main => kmath-core/src/commonMain}/kotlin/scientifik/kmath/structures/NDArray.kt (96%)
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt
 rename {kmath-common/src/test => kmath-core/src/commonTest}/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt (100%)
 rename {kmath-common/src/test => kmath-core/src/commonTest}/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt (100%)
 rename {kmath-common/src/test => kmath-core/src/commonTest}/kotlin/scientifik/kmath/operations/RealFieldTest.kt (100%)
 rename {kmath-jvm/src/test => kmath-core/src/commonTest}/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt (77%)
 rename {kmath-jvm/src/test => kmath-core/src/commonTest}/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt (74%)
 rename {kmath-common/src/test => kmath-core/src/commonTest}/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt (85%)
 rename {kmath-jvm => kmath-core}/src/jmh/kotlin/scietifik/kmath/structures/ArrayBenchmark.kt (100%)
 create mode 100644 kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
 rename kmath-jvm/src/main/kotlin/scientifik/kmath/structures/RealNDArray.kt => kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt (73%)
 delete mode 100644 kmath-jvm/build.gradle

diff --git a/README.md b/README.md
index 8f914d791..d5831b03a 100644
--- a/README.md
+++ b/README.md
@@ -40,8 +40,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.
\ No newline at end of file
diff --git a/build.gradle b/build.gradle
index a3af48ee0..655acfe2d 100644
--- a/build.gradle
+++ b/build.gradle
@@ -1,17 +1,19 @@
 buildscript {
-    ext.kotlin_version = '1.2.71'
+    ext.kotlin_version = '1.3.0-rc-116'
 
     repositories {
         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'
+allprojects {
     apply plugin: 'maven-publish'
     apply plugin: "com.jfrog.artifactory"
 
@@ -27,13 +29,11 @@ artifactory {
             repoKey = 'gradle-dev-local'
             username = "${artifactory_user}"
             password = "${artifactory_password}"
-            maven = true
-
         }
 
         defaults {
-            publications 'defaultPublication'
-            publishBuildInfo = true
+            publications('jvm', 'kotlinMultiplatform', 'metadata')
+            publishBuildInfo = false
             publishArtifacts = true
             publishPom = true
             publishIvy = false
@@ -44,7 +44,6 @@ artifactory {
             repoKey = 'gradle-dev'
             username = "${artifactory_user}"
             password = "${artifactory_password}"
-            maven = true
         }
     }
 }
\ No newline at end of file
diff --git a/kmath-common/build.gradle b/kmath-common/build.gradle
deleted file mode 100644
index eecc32818..000000000
--- a/kmath-common/build.gradle
+++ /dev/null
@@ -1,22 +0,0 @@
-plugins{
-    id "kotlin-platform-common"
-}
-
-description = "Platform-independent interfaces for kotlin maths"
-
-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"
-}
-
-kotlin {
-    experimental {
-        coroutines "enable"
-    }
-}
-
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArrayFactories.kt b/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArrayFactories.kt
deleted file mode 100644
index b3e3850e9..000000000
--- a/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArrayFactories.kt
+++ /dev/null
@@ -1,38 +0,0 @@
-package scientifik.kmath.structures
-
-import scientifik.kmath.operations.Field
-
-/**
- * Create a platform-optimized 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]) }
-}
-
-
-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 object : Buffer<T> {
-            override fun get(index: Int): T = array[index]
-
-            override fun set(index: Int, value: T) {
-                array[index] = initializer(index)
-            }
-        }
-    }
-}
-
-fun <T: Any> simpleNDArray(field: Field<T>, shape: List<Int>, initializer: (List<Int>) -> T): NDArray<T> {
-    return SimpleNDField(field, shape).produce { initializer(it) }
-}
\ No newline at end of file
diff --git a/kmath-core/build.gradle b/kmath-core/build.gradle
new file mode 100644
index 000000000..8e0e6ca2e
--- /dev/null
+++ b/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 {
+//        }
+    }
+}
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/expressions/Expression.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/expressions/Expression.kt
similarity index 100%
rename from kmath-common/src/main/kotlin/scientifik/kmath/expressions/Expression.kt
rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/expressions/Expression.kt
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/misc/Cumulative.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Cumulative.kt
similarity index 100%
rename from kmath-common/src/main/kotlin/scientifik/kmath/misc/Cumulative.kt
rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Cumulative.kt
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/misc/Grids.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Grids.kt
similarity index 90%
rename from kmath-common/src/main/kotlin/scientifik/kmath/misc/Grids.kt
rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Grids.kt
index 5ea4bd949..b48afbdcc 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/misc/Grids.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/misc/Grids.kt
@@ -1,8 +1,5 @@
 package scientifik.kmath.misc
 
-import kotlin.coroutines.experimental.buildSequence
-
-
 /**
  * Convert double range to sequence.
  *
@@ -14,14 +11,14 @@ import kotlin.coroutines.experimental.buildSequence
 fun ClosedFloatingPointRange<Double>.toSequence(step: Double): Sequence<Double> {
     return when {
         step == 0.0 -> error("Zero step in double progression")
-        step > 0 -> buildSequence {
+        step > 0 -> sequence {
             var current = start
             while (current <= endInclusive) {
                 yield(current)
                 current += step
             }
         }
-        else -> buildSequence {
+        else -> sequence {
             var current = endInclusive
             while (current >= start) {
                 yield(current)
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt
similarity index 100%
rename from kmath-common/src/main/kotlin/scientifik/kmath/operations/Algebra.kt
rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/operations/Fields.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
similarity index 100%
rename from kmath-common/src/main/kotlin/scientifik/kmath/operations/Fields.kt
rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/operations/OptionalOperations.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt
similarity index 100%
rename from kmath-common/src/main/kotlin/scientifik/kmath/operations/OptionalOperations.kt
rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/OptionalOperations.kt
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt
similarity index 96%
rename from kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt
rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt
index 5b3ef4c6a..9fe7efe87 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/structures/BufferNDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt
@@ -1,7 +1,6 @@
 package scientifik.kmath.structures
 
 import scientifik.kmath.operations.Field
-import kotlin.coroutines.experimental.buildSequence
 
 
 /**
@@ -42,7 +41,7 @@ abstract class BufferNDField<T>(shape: List<Int>, field: Field<T>) : NDField<T>(
 
     //TODO introduce a fast way to calculate index of the next element?
     protected fun index(offset: Int): List<Int> {
-        return buildSequence {
+        return sequence {
             var current = offset
             var strideIndex = strides.size - 2
             while (strideIndex >= 0) {
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt
similarity index 66%
rename from kmath-common/src/main/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt
rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt
index 62bad9b75..a5e6f2ce1 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt
@@ -4,6 +4,7 @@ import scientifik.kmath.operations.DoubleField
 import scientifik.kmath.operations.Field
 import scientifik.kmath.operations.Space
 import scientifik.kmath.operations.SpaceElement
+import scientifik.kmath.structures.NDArrays.createSimpleNDFieldFactory
 
 /**
  * The space for linear elements. Supports scalar product alongside with standard linear operations.
@@ -18,7 +19,7 @@ abstract class LinearSpace<T : Any, V : LinearStructure<out T>>(val rows: Int, v
     abstract fun produce(initializer: (Int, Int) -> T): V
 
     /**
-     * Produce new linear space with given dimensions
+     * Produce new linear space with given dimensions. The space produced could be raised from cache since [LinearSpace] does not have mutable elements
      */
     abstract fun produceSpace(rows: Int, columns: Int): LinearSpace<T, V>
 
@@ -36,12 +37,12 @@ abstract class LinearSpace<T : Any, V : LinearStructure<out T>>(val rows: Int, v
     }
 
     /**
-     * Dot product
+     * Dot product. Throws exception on dimension mismatch
      */
     fun multiply(a: V, b: V): V {
         if (a.rows != b.columns) {
             //TODO replace by specific exception
-            error("Dimension mismatch in vector dot product")
+            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..a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) }
@@ -55,16 +56,25 @@ abstract class LinearSpace<T : Any, V : LinearStructure<out T>>(val rows: Int, v
  * A matrix-like structure that is not dependent on specific space implementation
  */
 interface LinearStructure<T : Any> {
+    /**
+     * 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
 
     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)
+            override fun get(i: Int, j: Int): T = this@LinearStructure[j, i]
         }
     }
 }
@@ -78,15 +88,23 @@ interface Vector<T : Any> : LinearStructure<T> {
 
 
 /**
- * DoubleArray-based implementation of vector space
+ * NDArray-based implementation of vector space. By default uses slow [SimpleNDField], but could be overridden with custom [NDField] factory.
  */
-class ArraySpace<T : Any>(rows: Int, columns: Int, field: Field<T>) : LinearSpace<T, LinearStructure<out T>>(rows, columns, field) {
+class ArraySpace<T : Any>(
+        rows: Int,
+        columns: Int,
+        field: Field<T>,
+        val ndFactory: NDFieldFactory<T> = createSimpleNDFieldFactory(field)
+) : LinearSpace<T, LinearStructure<out T>>(rows, columns, field) {
+
+    val ndField by lazy {
+        ndFactory(listOf(rows, columns))
+    }
 
     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)
+        return ArraySpace(rows, columns, field, ndFactory)
     }
 }
 
@@ -95,14 +113,16 @@ class ArraySpace<T : Any>(rows: Int, columns: Int, field: Field<T>) : LinearSpac
  */
 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) } }
+    private val array = context.ndField.produce { list -> initializer(list[0], list[1]) }
+
+    //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]
+        return array[i, j]
     }
 
     override val self: ArrayMatrix<T> get() = this
@@ -117,7 +137,7 @@ class ArrayVector<T : Any>(override val context: ArraySpace<T>, initializer: (In
         }
     }
 
-    val list: List<T> = (0 until context.rows).map(initializer)
+    private val array = context.ndField.produce { list -> initializer(list[0]) }
 
 
     override val rows: Int get() = context.rows
@@ -125,20 +145,25 @@ class ArrayVector<T : Any>(override val context: ArraySpace<T>, initializer: (In
     override val columns: Int = 1
 
     override fun get(i: Int, j: Int): T {
-        return list[i]
+        return array[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> vector(size: Int, field: Field<T>, initializer: (Int) -> T) =
+        ArrayVector(ArraySpace(size, 1, field), initializer)
+
+fun realVector(size: Int, initializer: (Int) -> Double) =
+        ArrayVector(ArraySpace(size, 1, DoubleField, realNDFieldFactory), 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)
\ No newline at end of file
+fun <T : Any> matrix(rows: Int, columns: Int, field: Field<T>, initializer: (Int, Int) -> T) =
+        ArrayMatrix(ArraySpace(rows, columns, field), initializer)
+
+fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
+        ArrayMatrix(ArraySpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
diff --git a/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArray.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArray.kt
similarity index 96%
rename from kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArray.kt
rename to kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArray.kt
index 462cd6adb..b1ceed5b1 100644
--- a/kmath-common/src/main/kotlin/scientifik/kmath/structures/NDArray.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArray.kt
@@ -136,7 +136,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 +149,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 +159,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 +170,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 +183,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,8 +193,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
     }
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt
new file mode 100644
index 000000000..ee5ff0a74
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt
@@ -0,0 +1,54 @@
+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>
+
+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]) }
+    }
+
+    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 object : Buffer<T> {
+                override fun get(index: Int): T = array[index]
+
+                override fun set(index: Int, value: T) {
+                    array[index] = initializer(index)
+                }
+            }
+        }
+    }
+
+    fun <T : Any> createSimpleNDFieldFactory(field: Field<T>): NDFieldFactory<T> = { list -> SimpleNDField(field, list) }
+
+    fun <T : Any> simpleNDArray(field: Field<T>, shape: List<Int>, initializer: (List<Int>) -> T): NDArray<T> {
+        return SimpleNDField(field, shape).produce { initializer(it) }
+    }
+}
\ No newline at end of file
diff --git a/kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
similarity index 100%
rename from kmath-common/src/test/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
rename to kmath-core/src/commonTest/kotlin/scientifik/kmath/expressions/FieldExpressionContextTest.kt
diff --git a/kmath-common/src/test/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt
similarity index 100%
rename from kmath-common/src/test/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt
rename to kmath-core/src/commonTest/kotlin/scientifik/kmath/misc/CumulativeKtTest.kt
diff --git a/kmath-common/src/test/kotlin/scientifik/kmath/operations/RealFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt
similarity index 100%
rename from kmath-common/src/test/kotlin/scientifik/kmath/operations/RealFieldTest.kt
rename to kmath-core/src/commonTest/kotlin/scientifik/kmath/operations/RealFieldTest.kt
diff --git a/kmath-jvm/src/test/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt
similarity index 77%
rename from kmath-jvm/src/test/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt
rename to kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt
index 70fd8a6a9..10d383a22 100644
--- a/kmath-jvm/src/test/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt
@@ -1,7 +1,7 @@
 package scientifik.kmath.structures
 
-import org.junit.Assert.assertEquals
-import org.junit.Test
+import kotlin.test.Test
+import kotlin.test.assertEquals
 
 class ArrayMatrixTest {
 
@@ -10,7 +10,7 @@ class ArrayMatrixTest {
         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)
+        assertEquals(5.0, sum[2, 0])
     }
 
     @Test
@@ -21,6 +21,6 @@ class ArrayMatrixTest {
             vector1 dot (vector2.transpose())
         }
 
-        assertEquals(10.0, product[1, 0], 0.1)
+        assertEquals(10.0, product[1, 0])
     }
 }
\ No newline at end of file
diff --git a/kmath-jvm/src/test/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt
similarity index 74%
rename from kmath-jvm/src/test/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt
rename to kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt
index 9bed502f9..329f13189 100644
--- a/kmath-jvm/src/test/kotlin/scientifik/kmath/structures/RealNDFieldTest.kt
+++ b/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])
     }
 }
diff --git a/kmath-common/src/test/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
similarity index 85%
rename from kmath-common/src/test/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
rename to kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
index 5a4317251..21431dad5 100644
--- a/kmath-common/src/test/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
@@ -1,6 +1,7 @@
 package scientifik.kmath.structures
 
 import scientifik.kmath.operations.DoubleField
+import scientifik.kmath.structures.NDArrays.simpleNDArray
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
diff --git a/kmath-jvm/src/jmh/kotlin/scietifik/kmath/structures/ArrayBenchmark.kt b/kmath-core/src/jmh/kotlin/scietifik/kmath/structures/ArrayBenchmark.kt
similarity index 100%
rename from kmath-jvm/src/jmh/kotlin/scietifik/kmath/structures/ArrayBenchmark.kt
rename to kmath-core/src/jmh/kotlin/scietifik/kmath/structures/ArrayBenchmark.kt
diff --git a/kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt b/kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
new file mode 100644
index 000000000..29e64a575
--- /dev/null
+++ b/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.createSimpleNDFieldFactory(DoubleField)
\ No newline at end of file
diff --git a/kmath-jvm/src/main/kotlin/scientifik/kmath/structures/RealNDArray.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
similarity index 73%
rename from kmath-jvm/src/main/kotlin/scientifik/kmath/structures/RealNDArray.kt
rename to kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
index 8099e8700..c46d4210a 100644
--- a/kmath-jvm/src/main/kotlin/scientifik/kmath/structures/RealNDArray.kt
+++ b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
@@ -17,8 +17,4 @@ private class RealNDField(shape: List<Int>) : BufferNDField<Double>(shape, Doubl
     }
 }
 
-actual fun realNDArray(shape: List<Int>, initializer: (List<Int>) -> Double): NDArray<Double> {
-    //TODO create a cache for fields to save time generating strides?
-
-    return RealNDField(shape).produce { initializer(it) }
-}
\ No newline at end of file
+actual val realNDFieldFactory: NDFieldFactory<Double> = { shape -> RealNDField(shape) }
\ No newline at end of file
diff --git a/kmath-jvm/build.gradle b/kmath-jvm/build.gradle
deleted file mode 100644
index b8952a379..000000000
--- a/kmath-jvm/build.gradle
+++ /dev/null
@@ -1,38 +0,0 @@
-plugins{
-    id "kotlin-platform-jvm"
-    id "me.champeau.gradle.jmh" version "0.4.5"
-}
-
-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"
-
-kotlin {
-    experimental {
-        coroutines "enable"
-    }
-}
-
-publishing {
-    publications {
-        defaultPublication(MavenPublication) {
-            from components.java
-        }
-    }
-}
\ No newline at end of file
diff --git a/settings.gradle b/settings.gradle
index 9dce0ac08..03760221a 100644
--- a/settings.gradle
+++ b/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'
 
-- 
2.34.1


From 76c3025acc5a48b9274055d834b22f26859523c1 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Sun, 30 Sep 2018 17:25:11 +0300
Subject: [PATCH 14/17] Readme fix. Added js deploy target

---
 README.md    | 4 +---
 build.gradle | 2 +-
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/README.md b/README.md
index d5831b03a..7dfd1a038 100644
--- a/README.md
+++ b/README.md
@@ -21,7 +21,7 @@ for a wide variety of purposes from high performance calculations to code genera
 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 multilanguage and multinod communication for mathematical tasks.
+* **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.
@@ -30,8 +30,6 @@ platform even if they could be done in common module because of platform perform
 Currently the main focus of development is the JVM platform, contribution of implementations for Kotlin - Native and
 Kotlin - JS is welcome. 
 
-It is planned to move project structure to [kotlin gradle mpp structure](https://github.com/h0tk3y/k-new-mpp-samples) after release of kotlin 1.3. 
-
 ## Performance
 The calculation performance is one of major goals of KMath in the future, but in some cases it is not possible to achieve 
 both performance and flexibility. We expect to firstly focus on creating convenient universal API and then work on 
diff --git a/build.gradle b/build.gradle
index 655acfe2d..8db73094a 100644
--- a/build.gradle
+++ b/build.gradle
@@ -32,7 +32,7 @@ artifactory {
         }
 
         defaults {
-            publications('jvm', 'kotlinMultiplatform', 'metadata')
+            publications('jvm', 'js', 'kotlinMultiplatform', 'metadata')
             publishBuildInfo = false
             publishArtifacts = true
             publishPom = true
-- 
2.34.1


From a9c084773b3931c37ce5ec3db5d7d58abcaffc6d Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Tue, 9 Oct 2018 12:22:35 +0300
Subject: [PATCH 15/17] Linear algebra in progress

---
 build.gradle                                  |   2 +-
 .../scientifik/kmath/commons/CommonsMatrix.kt |   3 +
 .../kmath/linear/LUDecomposition.kt           | 321 ++++++++++++++++++
 .../scientifik/kmath/linear/LinearAlgrebra.kt | 250 ++++++++++++++
 .../scientifik/kmath/operations/Complex.kt    |  41 +++
 .../scientifik/kmath/operations/Fields.kt     |  41 ---
 .../kmath/structures/BufferNDField.kt         |  21 +-
 .../kmath/structures/LinearAlgrebra.kt        | 169 ---------
 .../scientifik/kmath/structures/NDArray.kt    |  20 +-
 .../scientifik/kmath/structures/NDArrays.kt   |  56 +--
 .../kmath/structures/ArrayMatrixTest.kt       |   1 +
 .../kmath/structures/SimpleNDFieldTest.kt     |   4 +-
 .../scientifik/kmath/structures/_NDArrays.kt  |   2 +-
 .../scientifik/kmath/structures/_NDArrays.kt  |  16 +-
 settings.gradle                               |   1 +
 15 files changed, 708 insertions(+), 240 deletions(-)
 create mode 100644 kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
 create mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt
 delete mode 100644 kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt

diff --git a/build.gradle b/build.gradle
index 8db73094a..04700f377 100644
--- a/build.gradle
+++ b/build.gradle
@@ -1,5 +1,5 @@
 buildscript {
-    ext.kotlin_version = '1.3.0-rc-116'
+    ext.kotlin_version = '1.3.0-rc-146'
 
     repositories {
         jcenter()
diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt
new file mode 100644
index 000000000..a4a9a67ae
--- /dev/null
+++ b/kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt
@@ -0,0 +1,3 @@
+package scientifik.kmath.commons
+
+//val solver: DecompositionSolver
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
new file mode 100644
index 000000000..6d329809f
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
@@ -0,0 +1,321 @@
+package scientifik.kmath.linear
+
+import scientifik.kmath.operations.Field
+import scientifik.kmath.structures.MutableNDArray
+import scientifik.kmath.structures.NDArray
+import scientifik.kmath.structures.NDArrays
+
+/**
+ * Calculates the LUP-decomposition of a square matrix.
+ *
+ * The LUP-decomposition of a matrix A consists of three matrices L, U and
+ * P that satisfy: PA = LU. L is lower triangular (with unit
+ * diagonal terms), U is upper triangular and P is a permutation matrix. All
+ * matrices are mm.
+ *
+ * As shown by the presence of the P matrix, this decomposition is
+ * implemented using partial pivoting.
+ *
+ * This class is based on the class with similar name from the
+ * [JAMA](http://math.nist.gov/javanumerics/jama/) library.
+ *
+ *  * a [getP][.getP] method has been added,
+ *  * the `det` method has been renamed as [   getDeterminant][.getDeterminant],
+ *  * the `getDoublePivot` method has been removed (but the int based
+ * [getPivot][.getPivot] method has been kept),
+ *  * the `solve` and `isNonSingular` methods have been replaced
+ * by a [getSolver][.getSolver] method and the equivalent methods
+ * provided by the returned [DecompositionSolver].
+ *
+ *
+ * @see [MathWorld](http://mathworld.wolfram.com/LUDecomposition.html)
+ *
+ * @see [Wikipedia](http://en.wikipedia.org/wiki/LU_decomposition)
+ *
+ * @since 2.0 (changed to concrete class in 3.0)
+ *
+ * @param matrix The matrix to decompose.
+ * @param singularityThreshold threshold (based on partial row norm)
+ * under which a matrix is considered singular
+ * @throws NonSquareMatrixException if matrix is not square
+
+ */
+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 = matrix.context.field.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 Field<T>.calculateLU(): Pair<NDArray<T>, IntArray> {
+        if (matrix.rows != matrix.columns) {
+            error("LU decomposition supports only square matrices")
+        }
+
+        fun T.abs() = if (this > zero) this else -this
+
+        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]] }
+
+        // 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
+
+                sum.abs()
+            } ?: col
+
+            // Singularity check
+            if (isSingular(lu[max, col].abs())) {
+                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] /= luDiag
+            }
+        }
+        return Pair(lu, pivot)
+    }
+
+    /**
+     * Returns the pivot permutation vector.
+     * @return the pivot permutation vector
+     * @see .getP
+     */
+    fun getPivot(): IntArray {
+        return pivot.copyOf()
+    }
+
+    companion object {
+        /** Default bound to determine effective singularity in LU decomposition.  */
+        private const val DEFAULT_TOO_SMALL = 1e-11
+    }
+}
+
+class RealLUDecomposition(matrix: Matrix<Double>, private val singularityThreshold: Double = 1e-11) : LUDecomposition<Double>(matrix) {
+    override fun isSingular(value: Double): Boolean {
+        return value < singularityThreshold
+    }
+}
+
+
+/** Specialized solver.  */
+class 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] }
+    }
+}
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
new file mode 100644
index 000000000..07295e0aa
--- /dev/null
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
@@ -0,0 +1,250 @@
+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 LinearSpace<T : Any, V : Matrix<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. The space produced could be raised from cache since [LinearSpace] does not have mutable elements
+     */
+    abstract fun produceSpace(rows: Int, columns: Int): LinearSpace<T, V>
+
+    override val zero: V by lazy {
+        produce { _, _ -> field.zero }
+    }
+
+    val one: V by lazy {
+        produce { i, j -> if (i == j) field.one else 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. Throws exception on dimension mismatch
+     */
+    fun multiply(a: V, b: V): V {
+        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..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 specialized [LinearSpace] which works with vectors
+ */
+abstract class VectorSpace<T : Any, V : Vector<T>>(size: Int, field: Field<T>) : LinearSpace<T, V>(size, 1, field)
+
+/**
+ * A matrix-like structure
+ */
+interface Matrix<T : Any> {
+    val context: LinearSpace<T, out Matrix<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
+
+    fun transpose(): Matrix<T> {
+        return object : Matrix<T> {
+            override val context: LinearSpace<T, out Matrix<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]
+        }
+    }
+}
+
+interface Vector<T : Any> : Matrix<T> {
+    override val context: VectorSpace<T, Vector<T>>
+    override val columns: Int
+        get() = 1
+
+    operator fun get(i: Int) = get(i, 0)
+}
+
+
+/**
+ * NDArray-based implementation of vector space. By default uses slow [SimpleNDField], but could be overridden with custom [NDField] factory.
+ */
+class ArraySpace<T : Any>(
+        rows: Int,
+        columns: Int,
+        field: Field<T>,
+        val ndFactory: NDFieldFactory<T> = createFactory(field)
+) : LinearSpace<T, Matrix<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): ArraySpace<T> {
+        return ArraySpace(rows, columns, field, ndFactory)
+    }
+}
+
+class ArrayVectorSpace<T : Any>(
+        size: Int,
+        field: Field<T>,
+        val ndFactory: NDFieldFactory<T> = createFactory(field)
+) : VectorSpace<T, Vector<T>>(size, field) {
+    val ndField by lazy {
+        ndFactory(listOf(size))
+    }
+
+    override fun produce(initializer: (Int, Int) -> T): Vector<T> = produceVector { i -> initializer(i, 0) }
+
+    fun produceVector(initializer: (Int) -> T): Vector<T> = ArrayVector(this, initializer)
+
+    override fun produceSpace(rows: Int, columns: Int): LinearSpace<T, Vector<T>> {
+        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
+    }
+}
+
+/**
+ * Member of [ArraySpace] which wraps 2-D array
+ */
+class ArrayMatrix<T : Any> internal constructor(override val context: ArraySpace<T>, val array: NDArray<T>) : Matrix<T>, SpaceElement<Matrix<T>, ArraySpace<T>> {
+
+    constructor(context: ArraySpace<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>, SpaceElement<Vector<T>, ArrayVectorSpace<T>> {
+
+    constructor(context: ArrayVectorSpace<T>, initializer: (Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0]) })
+
+    init {
+        if (context.columns != 1) {
+            error("Vector must have single column")
+        }
+        if (context.rows != array.shape[0]) {
+            error("Array dimension mismatch")
+        }
+    }
+
+    //private val array = context.ndField.produce { list -> initializer(list[0]) }
+
+
+    override val rows: Int get() = context.rows
+
+    override val columns: Int = 1
+
+    override fun get(i: Int, j: 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 as Matrix<T>).toVector()
+    fun inverse(a: Matrix<T>): Matrix<T> = solve(a, a.context.one)
+}
+
+/**
+ * Create vector with custom field
+ */
+fun <T : Any> vector(size: Int, field: Field<T>, initializer: (Int) -> T) =
+        ArrayVector(ArrayVectorSpace(size, field), initializer)
+
+/**
+ * Create vector of [Double]
+ */
+fun realVector(size: Int, initializer: (Int) -> Double) =
+        ArrayVector(ArrayVectorSpace(size, DoubleField, realNDFieldFactory), initializer)
+
+/**
+ * Convert vector to array (copying content of array)
+ */
+fun <T : Any> Array<T>.asVector(field: Field<T>) = vector(size, field) { this[it] }
+
+fun DoubleArray.asVector() = realVector(this.size) { this[it] }
+
+/**
+ * Create [ArrayMatrix] with custom field
+ */
+fun <T : Any> matrix(rows: Int, columns: Int, field: Field<T>, initializer: (Int, Int) -> T) =
+        ArrayMatrix(ArraySpace(rows, columns, field), initializer)
+
+/**
+ * Create [ArrayMatrix] of doubles.
+ */
+fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
+        ArrayMatrix(ArraySpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
+
+
+/**
+ * Convert matrix to vector if it is possible
+ */
+fun <T : Any> Matrix<T>.toVector(): Vector<T> {
+    return when {
+        this is Vector -> return this
+        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) }
+            }
+        }
+        else -> error("Can't convert matrix with more than one column to vector")
+    }
+}
+
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Complex.kt
new file mode 100644
index 000000000..dd3be256e
--- /dev/null
+++ b/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)
+
+}
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
index 0d2600d82..eeba76222 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Fields.kt
+++ b/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
@@ -51,46 +50,6 @@ data 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
  */
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt
index 9fe7efe87..aab0f1b97 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/BufferNDField.kt
@@ -9,6 +9,11 @@ import scientifik.kmath.operations.Field
 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>
 }
 
 /**
@@ -63,8 +68,16 @@ abstract class BufferNDField<T>(shape: List<Int>, field: Field<T>) : NDField<T>(
         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)
+    }
 
-    class BufferNDArray<T>(override val context: BufferNDField<T>, val data: Buffer<T>) : NDArray<T> {
+
+    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())]
@@ -88,6 +101,12 @@ abstract class BufferNDField<T>(shape: List<Int>, field: Field<T>) : NDField<T>(
 
         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
+        }
+    }
 }
 
 
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt
deleted file mode 100644
index a5e6f2ce1..000000000
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/LinearAlgrebra.kt
+++ /dev/null
@@ -1,169 +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
-import scientifik.kmath.structures.NDArrays.createSimpleNDFieldFactory
-
-/**
- * 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. The space produced could be raised from cache since [LinearSpace] does not have mutable elements
-     */
-    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. Throws exception on dimension mismatch
-     */
-    fun multiply(a: V, b: V): V {
-        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..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> {
-    /**
-     * 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
-
-    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[j, i]
-        }
-    }
-}
-
-interface Vector<T : Any> : LinearStructure<T> {
-    override val columns: Int
-        get() = 1
-
-    operator fun get(i: Int) = get(i, 0)
-}
-
-
-/**
- * NDArray-based implementation of vector space. By default uses slow [SimpleNDField], but could be overridden with custom [NDField] factory.
- */
-class ArraySpace<T : Any>(
-        rows: Int,
-        columns: Int,
-        field: Field<T>,
-        val ndFactory: NDFieldFactory<T> = createSimpleNDFieldFactory(field)
-) : LinearSpace<T, LinearStructure<out T>>(rows, columns, field) {
-
-    val ndField by lazy {
-        ndFactory(listOf(rows, columns))
-    }
-
-    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, ndFactory)
-    }
-}
-
-/**
- * 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>> {
-
-    private val array = context.ndField.produce { list -> initializer(list[0], list[1]) }
-
-    //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 array[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")
-        }
-    }
-
-    private val array = context.ndField.produce { list -> initializer(list[0]) }
-
-
-    override val rows: Int get() = context.rows
-
-    override val columns: Int = 1
-
-    override fun get(i: Int, j: Int): T {
-        return array[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)
-
-fun realVector(size: Int, initializer: (Int) -> Double) =
-        ArrayVector(ArraySpace(size, 1, DoubleField, realNDFieldFactory), initializer)
-
-fun <T : Any> Array<T>.asVector(field: Field<T>) = vector(size, field) { this[it] }
-
-fun DoubleArray.asVector() = realVector(this.size) { this[it] }
-
-fun <T : Any> matrix(rows: Int, columns: Int, field: Field<T>, initializer: (Int, Int) -> T) =
-        ArrayMatrix(ArraySpace(rows, columns, field), initializer)
-
-fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
-        ArrayMatrix(ArraySpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArray.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArray.kt
index b1ceed5b1..85f60704f 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArray.kt
+++ b/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
  */
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt
index ee5ff0a74..739e91832 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/structures/NDArrays.kt
@@ -9,6 +9,28 @@ typealias NDFieldFactory<T> = (shape: List<Int>) -> NDField<T>
  */
 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
@@ -29,26 +51,22 @@ object NDArrays {
         return realNDArray(listOf(dim1, dim2, dim3)) { initializer(it[0], it[1], it[2]) }
     }
 
-    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))
-            }
+    /**
+     * Simple boxing NDField
+     */
+    fun <T : Any> createFactory(field: Field<T>): NDFieldFactory<T> = { shape -> SimpleNDField(field, shape) }
 
-            return object : Buffer<T> {
-                override fun get(index: Int): T = array[index]
-
-                override fun set(index: Int, value: T) {
-                    array[index] = initializer(index)
-                }
-            }
-        }
-    }
-
-    fun <T : Any> createSimpleNDFieldFactory(field: Field<T>): NDFieldFactory<T> = { list -> SimpleNDField(field, list) }
-
-    fun <T : Any> simpleNDArray(field: Field<T>, shape: List<Int>, initializer: (List<Int>) -> T): NDArray<T> {
+    /**
+     * 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) }
+    }
 }
\ No newline at end of file
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt
index 10d383a22..77a436626 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt
@@ -1,5 +1,6 @@
 package scientifik.kmath.structures
 
+import scientifik.kmath.linear.realVector
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
index 21431dad5..94c1e15cc 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/SimpleNDFieldTest.kt
@@ -1,7 +1,7 @@
 package scientifik.kmath.structures
 
 import scientifik.kmath.operations.DoubleField
-import scientifik.kmath.structures.NDArrays.simpleNDArray
+import scientifik.kmath.structures.NDArrays.create
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
@@ -9,7 +9,7 @@ import kotlin.test.assertEquals
 class SimpleNDFieldTest{
     @Test
     fun testStrides(){
-        val ndArray = simpleNDArray(DoubleField, listOf(10,10)){(it[0]+it[1]).toDouble()}
+        val ndArray = create(DoubleField, listOf(10,10)){(it[0]+it[1]).toDouble()}
         assertEquals(ndArray[5,5], 10.0)
     }
 
diff --git a/kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt b/kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
index 29e64a575..32d66a04b 100644
--- a/kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
+++ b/kmath-core/src/jsMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
@@ -5,4 +5,4 @@ import scientifik.kmath.operations.DoubleField
 /**
  * Using boxing implementation for js
  */
-actual val realNDFieldFactory: NDFieldFactory<Double> = NDArrays.createSimpleNDFieldFactory(DoubleField)
\ No newline at end of file
+actual val realNDFieldFactory: NDFieldFactory<Double> = NDArrays.createFactory(DoubleField)
\ No newline at end of file
diff --git a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
index c46d4210a..19d6c4041 100644
--- a/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
+++ b/kmath-core/src/jvmMain/kotlin/scientifik/kmath/structures/_NDArrays.kt
@@ -7,12 +7,18 @@ private class RealNDField(shape: List<Int>) : BufferNDField<Double>(shape, Doubl
     override fun createBuffer(capacity: Int, initializer: (Int) -> Double): Buffer<Double> {
         val array = DoubleArray(capacity, initializer)
         val buffer = DoubleBuffer.wrap(array)
-        return object : Buffer<Double> {
-            override fun get(index: Int): Double = buffer.get(index)
+        return BufferOfDoubles(buffer)
+    }
 
-            override fun set(index: Int, value: Double) {
-                buffer.put(index, value)
-            }
+    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)
         }
     }
 }
diff --git a/settings.gradle b/settings.gradle
index 03760221a..988a94843 100644
--- a/settings.gradle
+++ b/settings.gradle
@@ -10,4 +10,5 @@ enableFeaturePreview('GRADLE_METADATA')
 
 rootProject.name = 'kmath'
 include ':kmath-core'
+include ':kmath-commons'
 
-- 
2.34.1


From 0cde0cfea5d793c31b49390da1680b2ea546995f Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Thu, 11 Oct 2018 16:23:33 +0300
Subject: [PATCH 16/17] Basic matrix inversion

---
 .../scientifik/kmath/commons/CommonsMatrix.kt |   3 -
 .../kmath/linear/LUDecomposition.kt           | 153 ++++++++----------
 .../scientifik/kmath/linear/LinearAlgrebra.kt | 146 +++++++++--------
 .../scientifik/kmath/operations/Algebra.kt    |  11 +-
 .../{structures => linear}/ArrayMatrixTest.kt |  21 ++-
 .../kmath/linear/RealLUSolverTest.kt          |  14 ++
 settings.gradle                               |   1 -
 7 files changed, 175 insertions(+), 174 deletions(-)
 delete mode 100644 kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt
 rename kmath-core/src/commonTest/kotlin/scientifik/kmath/{structures => linear}/ArrayMatrixTest.kt (50%)
 create mode 100644 kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt

diff --git a/kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt b/kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt
deleted file mode 100644
index a4a9a67ae..000000000
--- a/kmath-commons/src/main/kotlin/scientifik/kmath/commons/CommonsMatrix.kt
+++ /dev/null
@@ -1,3 +0,0 @@
-package scientifik.kmath.commons
-
-//val solver: DecompositionSolver
\ No newline at end of file
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
index 6d329809f..e51d574f8 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LUDecomposition.kt
@@ -1,44 +1,12 @@
 package scientifik.kmath.linear
 
-import scientifik.kmath.operations.Field
 import scientifik.kmath.structures.MutableNDArray
 import scientifik.kmath.structures.NDArray
 import scientifik.kmath.structures.NDArrays
+import kotlin.math.absoluteValue
 
 /**
- * Calculates the LUP-decomposition of a square matrix.
- *
- * The LUP-decomposition of a matrix A consists of three matrices L, U and
- * P that satisfy: PA = LU. L is lower triangular (with unit
- * diagonal terms), U is upper triangular and P is a permutation matrix. All
- * matrices are mm.
- *
- * As shown by the presence of the P matrix, this decomposition is
- * implemented using partial pivoting.
- *
- * This class is based on the class with similar name from the
- * [JAMA](http://math.nist.gov/javanumerics/jama/) library.
- *
- *  * a [getP][.getP] method has been added,
- *  * the `det` method has been renamed as [   getDeterminant][.getDeterminant],
- *  * the `getDoublePivot` method has been removed (but the int based
- * [getPivot][.getPivot] method has been kept),
- *  * the `solve` and `isNonSingular` methods have been replaced
- * by a [getSolver][.getSolver] method and the equivalent methods
- * provided by the returned [DecompositionSolver].
- *
- *
- * @see [MathWorld](http://mathworld.wolfram.com/LUDecomposition.html)
- *
- * @see [Wikipedia](http://en.wikipedia.org/wiki/LU_decomposition)
- *
- * @since 2.0 (changed to concrete class in 3.0)
- *
- * @param matrix The matrix to decompose.
- * @param singularityThreshold threshold (based on partial row norm)
- * under which a matrix is considered singular
- * @throws NonSquareMatrixException if matrix is not square
-
+ * Implementation copier from Apache common-maths
  */
 abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
 
@@ -51,7 +19,7 @@ abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
     private var even: Boolean = false
 
     init {
-        val pair = matrix.context.field.calculateLU()
+        val pair = calculateLU()
         lu = pair.first
         pivot = pair.second
     }
@@ -134,72 +102,76 @@ abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
 
     abstract fun isSingular(value: T): Boolean
 
-    private fun Field<T>.calculateLU(): Pair<NDArray<T>, IntArray> {
+    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")
         }
 
-        fun T.abs() = if (this > zero) this else -this
-
         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]] }
 
-        // 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
+        with(matrix.context.field) {
+            // Initialize permutation array and parity
+            for (row in 0 until m) {
+                pivot[row] = row
             }
+            even = true
 
-            // 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]
+            // 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
                 }
-                //luRow[col] = sum
-                lu[row, col] = sum
 
-                sum.abs()
-            } ?: col
+                // 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
 
-            // Singularity check
-            if (isSingular(lu[max, col].abs())) {
-                error("Singular matrix")
-            }
+                    abs(sum)
+                } ?: col
 
-            // 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]
+                // Singularity check
+                if (isSingular(lu[max, col])) {
+                    error("Singular matrix")
                 }
-                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] /= luDiag
+                // 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)
@@ -214,21 +186,22 @@ abstract class LUDecomposition<T : Comparable<T>>(val matrix: Matrix<T>) {
         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
     }
 }
 
-class RealLUDecomposition(matrix: Matrix<Double>, private val singularityThreshold: Double = 1e-11) : LUDecomposition<Double>(matrix) {
-    override fun isSingular(value: Double): Boolean {
-        return value < singularityThreshold
-    }
-}
-
 
 /** Specialized solver.  */
-class RealLUSolver : LinearSolver<Double> {
+object RealLUSolver : LinearSolver<Double> {
 
 //
 //    /** {@inheritDoc}  */
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
index 07295e0aa..e7345163f 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
@@ -14,31 +14,31 @@ import scientifik.kmath.structures.realNDFieldFactory
  * @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 : Matrix<T>>(val rows: Int, val columns: Int, val field: Field<T>) : Space<V> {
+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): V
+    abstract fun produce(initializer: (Int, Int) -> T): Matrix<T>
 
     /**
-     * Produce new linear space with given dimensions. The space produced could be raised from cache since [LinearSpace] does not have mutable elements
+     * 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): LinearSpace<T, V>
+    abstract fun produceSpace(rows: Int, columns: Int): MatrixSpace<T>
 
-    override val zero: V by lazy {
+    override val zero: Matrix<T> by lazy {
         produce { _, _ -> field.zero }
     }
 
-    val one: V by lazy {
-        produce { i, j -> if (i == j) field.one else field.zero }
-    }
+//    val one: Matrix<T> by lazy {
+//        produce { i, j -> if (i == j) field.one else field.zero }
+//    }
 
-    override fun add(a: V, b: V): V {
+    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: V, k: Double): V {
+    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 } }
     }
@@ -46,29 +46,23 @@ abstract class LinearSpace<T : Any, V : Matrix<T>>(val rows: Int, val columns: I
     /**
      * Dot product. Throws exception on dimension mismatch
      */
-    fun multiply(a: V, b: V): V {
+    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..a.columns).asSequence().map { k -> field.multiply(a[i, k], b[k, j]) }.reduce { first, second -> field.add(first, second) }
+            (0 until 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 specialized [LinearSpace] which works with vectors
- */
-abstract class VectorSpace<T : Any, V : Vector<T>>(size: Int, field: Field<T>) : LinearSpace<T, V>(size, 1, field)
+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> {
-    val context: LinearSpace<T, out Matrix<T>>
+interface Matrix<T : Any> : SpaceElement<Matrix<T>, MatrixSpace<T>> {
     /**
      * Number of rows
      */
@@ -83,34 +77,58 @@ interface Matrix<T : Any> {
      */
     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: LinearSpace<T, out Matrix<T>> = this@Matrix.context
+            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 {
+        fun <T : Any> one(rows: Int, columns: Int, field: Field<T>): Matrix<T> {
+            return matrix(rows, columns, field) { i, j -> if (i == j) field.one else field.zero }
+        }
+    }
 }
 
-interface Vector<T : Any> : Matrix<T> {
-    override val context: VectorSpace<T, Vector<T>>
-    override val columns: Int
-        get() = 1
 
-    operator fun get(i: Int) = get(i, 0)
+/**
+ * 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
 }
 
 
 /**
  * NDArray-based implementation of vector space. By default uses slow [SimpleNDField], but could be overridden with custom [NDField] factory.
  */
-class ArraySpace<T : Any>(
+class ArrayMatrixSpace<T : Any>(
         rows: Int,
         columns: Int,
         field: Field<T>,
         val ndFactory: NDFieldFactory<T> = createFactory(field)
-) : LinearSpace<T, Matrix<T>>(rows, columns, field) {
+) : MatrixSpace<T>(rows, columns, field) {
 
     val ndField by lazy {
         ndFactory(listOf(rows, columns))
@@ -118,8 +136,8 @@ class ArraySpace<T : Any>(
 
     override fun produce(initializer: (Int, Int) -> T): Matrix<T> = ArrayMatrix(this, initializer)
 
-    override fun produceSpace(rows: Int, columns: Int): ArraySpace<T> {
-        return ArraySpace(rows, columns, field, ndFactory)
+    override fun produceSpace(rows: Int, columns: Int): ArrayMatrixSpace<T> {
+        return ArrayMatrixSpace(rows, columns, field, ndFactory)
     }
 }
 
@@ -127,26 +145,20 @@ class ArrayVectorSpace<T : Any>(
         size: Int,
         field: Field<T>,
         val ndFactory: NDFieldFactory<T> = createFactory(field)
-) : VectorSpace<T, Vector<T>>(size, field) {
+) : VectorSpace<T>(size, field) {
     val ndField by lazy {
         ndFactory(listOf(size))
     }
 
-    override fun produce(initializer: (Int, Int) -> T): Vector<T> = produceVector { i -> initializer(i, 0) }
-
-    fun produceVector(initializer: (Int) -> T): Vector<T> = ArrayVector(this, initializer)
-
-    override fun produceSpace(rows: Int, columns: Int): LinearSpace<T, Vector<T>> {
-        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
-    }
+    override fun produce(initializer: (Int) -> T): Vector<T> = ArrayVector(this, initializer)
 }
 
 /**
- * Member of [ArraySpace] which wraps 2-D array
+ * Member of [ArrayMatrixSpace] which wraps 2-D array
  */
-class ArrayMatrix<T : Any> internal constructor(override val context: ArraySpace<T>, val array: NDArray<T>) : Matrix<T>, SpaceElement<Matrix<T>, ArraySpace<T>> {
+class ArrayMatrix<T : Any> internal constructor(override val context: ArrayMatrixSpace<T>, val array: NDArray<T>) : Matrix<T> {
 
-    constructor(context: ArraySpace<T>, initializer: (Int, Int) -> T) : this(context, context.ndField.produce { list -> initializer(list[0], list[1]) })
+    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
 
@@ -160,32 +172,21 @@ class ArrayMatrix<T : Any> internal constructor(override val context: ArraySpace
 }
 
 
-class ArrayVector<T : Any> internal constructor(override val context: ArrayVectorSpace<T>, val array: NDArray<T>) : Vector<T>, SpaceElement<Vector<T>, ArrayVectorSpace<T>> {
+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.columns != 1) {
-            error("Vector must have single column")
-        }
-        if (context.rows != array.shape[0]) {
+        if (context.size != array.shape[0]) {
             error("Array dimension mismatch")
         }
     }
 
-    //private val array = context.ndField.produce { list -> initializer(list[0]) }
-
-
-    override val rows: Int get() = context.rows
-
-    override val columns: Int = 1
-
-    override fun get(i: Int, j: Int): T {
+    override fun get(i: Int): T {
         return array[i]
     }
 
     override val self: ArrayVector<T> get() = this
-
 }
 
 /**
@@ -193,8 +194,8 @@ class ArrayVector<T : Any> internal constructor(override val context: ArrayVecto
  */
 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 as Matrix<T>).toVector()
-    fun inverse(a: Matrix<T>): Matrix<T> = solve(a, a.context.one)
+    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.one(a.rows, a.columns, a.context.field))
 }
 
 /**
@@ -220,13 +221,13 @@ fun DoubleArray.asVector() = realVector(this.size) { this[it] }
  * Create [ArrayMatrix] with custom field
  */
 fun <T : Any> matrix(rows: Int, columns: Int, field: Field<T>, initializer: (Int, Int) -> T) =
-        ArrayMatrix(ArraySpace(rows, columns, field), initializer)
+        ArrayMatrix(ArrayMatrixSpace(rows, columns, field), initializer)
 
 /**
  * Create [ArrayMatrix] of doubles.
  */
 fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
-        ArrayMatrix(ArraySpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
+        ArrayMatrix(ArrayMatrixSpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
 
 
 /**
@@ -234,17 +235,28 @@ fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
  */
 fun <T : Any> Matrix<T>.toVector(): Vector<T> {
     return when {
-        this is Vector -> return this
         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) }
-            }
+//            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(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(size, 1, context.field) { i, j -> get(i) }
+}
+
diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt
index 9938f4856..9a56a4e91 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/operations/Algebra.kt
@@ -7,6 +7,11 @@ package scientifik.kmath.operations
  * @param S the type of mathematical context for this element
  */
 interface MathElement<T, S> {
+    /**
+     * Self value. Needed for static type checking.
+     */
+    val self: T
+
     /**
      * The context this element belongs to
      */
@@ -54,12 +59,6 @@ interface Space<T> {
  * @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
-
     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())
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/ArrayMatrixTest.kt
similarity index 50%
rename from kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt
rename to kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/ArrayMatrixTest.kt
index 77a436626..ecc2ca28c 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/structures/ArrayMatrixTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/ArrayMatrixTest.kt
@@ -1,6 +1,5 @@
-package scientifik.kmath.structures
+package scientifik.kmath.linear
 
-import scientifik.kmath.linear.realVector
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
@@ -11,17 +10,25 @@ class ArrayMatrixTest {
         val vector1 = realVector(5) { it.toDouble() }
         val vector2 = realVector(5) { 5 - it.toDouble() }
         val sum = vector1 + vector2
-        assertEquals(5.0, sum[2, 0])
+        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 = with(vector1.context) {
-            vector1 dot (vector2.transpose())
-        }
+        val product = vector1.toMatrix() dot (vector2.toMatrix().transpose())
 
-        assertEquals(10.0, product[1, 0])
+
+        assertEquals(5.0, product[1, 0])
+        assertEquals(6.0, product[2, 2])
     }
 }
\ No newline at end of file
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
new file mode 100644
index 000000000..d2237b80a
--- /dev/null
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
@@ -0,0 +1,14 @@
+package scientifik.kmath.linear
+
+import scientifik.kmath.operations.DoubleField
+import kotlin.test.Test
+import kotlin.test.assertEquals
+
+class RealLUSolverTest {
+    @Test
+    fun testInvert() {
+        val matrix = Matrix.one(2, 2, DoubleField)
+        val inverted = RealLUSolver.inverse(matrix)
+        assertEquals(1.0, inverted[0, 0])
+    }
+}
\ No newline at end of file
diff --git a/settings.gradle b/settings.gradle
index 988a94843..03760221a 100644
--- a/settings.gradle
+++ b/settings.gradle
@@ -10,5 +10,4 @@ enableFeaturePreview('GRADLE_METADATA')
 
 rootProject.name = 'kmath'
 include ':kmath-core'
-include ':kmath-commons'
 
-- 
2.34.1


From f894175897864fac8f204e19da64401e1d5fd5c3 Mon Sep 17 00:00:00 2001
From: Alexander Nozik <altavir@gmail.com>
Date: Fri, 12 Oct 2018 10:49:11 +0300
Subject: [PATCH 17/17] Basic matrix inversion

---
 .../scientifik/kmath/linear/LinearAlgrebra.kt | 113 +++++++++++++-----
 .../kmath/linear/RealLUSolverTest.kt          |  15 ++-
 2 files changed, 92 insertions(+), 36 deletions(-)

diff --git a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
index e7345163f..68961a2db 100644
--- a/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
+++ b/kmath-core/src/commonMain/kotlin/scientifik/kmath/linear/LinearAlgrebra.kt
@@ -55,6 +55,26 @@ abstract class MatrixSpace<T : Any>(val rows: Int, val columns: Int, val field:
             (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)
@@ -90,8 +110,38 @@ interface Matrix<T : Any> : SpaceElement<Matrix<T>, MatrixSpace<T>> {
     }
 
     companion object {
-        fun <T : Any> one(rows: Int, columns: Int, field: Field<T>): Matrix<T> {
-            return matrix(rows, columns, field) { i, j -> if (i == j) field.one else field.zero }
+
+        /**
+         * 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
         }
     }
 }
@@ -117,6 +167,30 @@ interface Vector<T : Any> : SpaceElement<Vector<T>, VectorSpace<T>> {
         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
+        }
+    }
 }
 
 
@@ -195,40 +269,15 @@ class ArrayVector<T : Any> internal constructor(override val context: ArrayVecto
 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.one(a.rows, a.columns, a.context.field))
+    fun inverse(a: Matrix<T>): Matrix<T> = solve(a, Matrix.diagonal(a.rows, a.columns, a.context.field))
 }
 
-/**
- * Create vector with custom field
- */
-fun <T : Any> vector(size: Int, field: Field<T>, initializer: (Int) -> T) =
-        ArrayVector(ArrayVectorSpace(size, field), initializer)
-
-/**
- * Create vector of [Double]
- */
-fun realVector(size: Int, initializer: (Int) -> Double) =
-        ArrayVector(ArrayVectorSpace(size, DoubleField, realNDFieldFactory), initializer)
-
 /**
  * Convert vector to array (copying content of array)
  */
-fun <T : Any> Array<T>.asVector(field: Field<T>) = vector(size, field) { this[it] }
-
-fun DoubleArray.asVector() = realVector(this.size) { this[it] }
-
-/**
- * Create [ArrayMatrix] with custom field
- */
-fun <T : Any> matrix(rows: Int, columns: Int, field: Field<T>, initializer: (Int, Int) -> T) =
-        ArrayMatrix(ArrayMatrixSpace(rows, columns, field), initializer)
-
-/**
- * Create [ArrayMatrix] of doubles.
- */
-fun realMatrix(rows: Int, columns: Int, initializer: (Int, Int) -> Double) =
-        ArrayMatrix(ArrayMatrixSpace(rows, columns, DoubleField, realNDFieldFactory), initializer)
+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
@@ -243,7 +292,7 @@ fun <T : Any> Matrix<T>.toVector(): Vector<T> {
 //                //Generic vector
 //                vector(rows, context.field) { get(it, 0) }
 //            }
-            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")
     }
@@ -257,6 +306,6 @@ fun <T : Any> Vector<T>.toMatrix(): Matrix<T> {
 //        //Generic vector
 //        matrix(size, 1, context.field) { i, j -> get(i) }
 //    }
-    return matrix(size, 1, context.field) { i, j -> get(i) }
+    return Matrix.of(size, 1, context.field) { i, j -> get(i) }
 }
 
diff --git a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
index d2237b80a..472d5262f 100644
--- a/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
+++ b/kmath-core/src/commonTest/kotlin/scientifik/kmath/linear/RealLUSolverTest.kt
@@ -2,13 +2,20 @@ package scientifik.kmath.linear
 
 import scientifik.kmath.operations.DoubleField
 import kotlin.test.Test
-import kotlin.test.assertEquals
+import kotlin.test.assertTrue
 
 class RealLUSolverTest {
     @Test
-    fun testInvert() {
-        val matrix = Matrix.one(2, 2, DoubleField)
+    fun testInvertOne() {
+        val matrix = Matrix.diagonal(2, 2, DoubleField)
         val inverted = RealLUSolver.inverse(matrix)
-        assertEquals(1.0, inverted[0, 0])
+        assertTrue { Matrix.equals(matrix,inverted) }
     }
+
+//    @Test
+//    fun testInvert() {
+//        val matrix = realMatrix(2,2){}
+//        val inverted = RealLUSolver.inverse(matrix)
+//        assertTrue { Matrix.equals(matrix,inverted) }
+//    }
 }
\ No newline at end of file
-- 
2.34.1