Merge branch 'dev' into mp-samplers

# Conflicts:
#	examples/src/main/kotlin/kscience/kmath/commons/prob/DistributionBenchmark.kt
#	examples/src/main/kotlin/kscience/kmath/commons/prob/DistributionDemo.kt
#	kmath-commons/src/main/kotlin/kscience/kmath/commons/expressions/DiffExpression.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/distributions/NormalDistribution.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/internal/InternalErf.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/internal/InternalGamma.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/internal/InternalUtils.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/AhrensDieterExponentialSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/AhrensDieterMarsagliaTsangGammaSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/AliasMethodDiscreteSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/BoxMullerNormalizedGaussianSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/GaussianSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/KempSmallMeanPoissonSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/LargeMeanPoissonSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/MarsagliaNormalizedGaussianSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/NormalizedGaussianSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/PoissonSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/SmallMeanPoissonSampler.kt
#	kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/samplers/ZigguratNormalizedGaussianSampler.kt
#	kmath-stat/src/jvmMain/kotlin/kscience/kmath/stat/distributions.kt

.space.kts

@@ -1 +1,3 @@
-job("Build") { gradlew("openjdk:11", "build") }
+job("Build") {
+    gradlew("openjdk:11", "build")
+}

CHANGELOG.md

@@ -7,7 +7,13 @@
 - Better trigonometric and hyperbolic functions for `AutoDiffField` (https://github.com/mipt-npm/kmath/pull/140).
 - Automatic README generation for features (#139)
 - Native support for `memory`, `core` and `dimensions`
-- `kmath-ejml` to supply EJML SimpleMatrix wrapper.
+- `kmath-ejml` to supply EJML SimpleMatrix wrapper (https://github.com/mipt-npm/kmath/pull/136).
+- A separate `Symbol` entity, which is used for global unbound symbol.
+- A `Symbol` indexing scope.
+- Basic optimization API for Commons-math.
+- Chi squared optimization for array-like data in CM
+- `Fitting` utility object in prob/stat
+- ND4J support module submitting `NDStructure` and `NDAlgebra` over `INDArray`.
 
 ### Changed
 - Package changed from `scientifik` to `kscience.kmath`.
@@ -16,6 +22,8 @@
 - `Polynomial` secondary constructor made function.
 - Kotlin version: 1.3.72 -> 1.4.20-M1
 - `kmath-ast` doesn't depend on heavy `kotlin-reflect` library.
+- Full autodiff refactoring based on `Symbol`
+- `kmath-prob` renamed to `kmath-stat`
 
 ### Deprecated
 

README.md

@@ -8,41 +8,50 @@ Bintray:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience
 Bintray-dev:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-core/_latestVersion)
 
 # KMath
-Could be pronounced as `key-math`.
+
-The Kotlin MATHematics library  was initially intended as a Kotlin-based analog to Python's `numpy` library. Later we found that kotlin is much more flexible language and allows superior architecture designs. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. The `numpy`-like experience could be achieved with [kmath-for-real](/kmath-for-real) extension module.
+Could be pronounced as `key-math`. The Kotlin MATHematics library was initially intended as a Kotlin-based analog to 
+Python's NumPy library. Later we found that kotlin is much more flexible language and allows superior architecture 
+designs. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. The `numpy`-like experience could 
+be achieved with [kmath-for-real](/kmath-for-real) extension module.
 
 ## Publications and talks
+
 * [A conceptual article about context-oriented design](https://proandroiddev.com/an-introduction-context-oriented-programming-in-kotlin-2e79d316b0a2)
 * [Another article about context-oriented design](https://proandroiddev.com/diving-deeper-into-context-oriented-programming-in-kotlin-3ecb4ec38814)
 * [ACAT 2019 conference paper](https://aip.scitation.org/doi/abs/10.1063/1.5130103)
 
 # Goal
-* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM and JS for now and Native in future). 
+
+* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM, JS and Native). 
 * Provide basic multiplatform implementations for those abstractions (without significant performance optimization).
 * Provide bindings and wrappers with those abstractions for popular optimized platform libraries.
 
 ## Non-goals
-* Be like Numpy. It was the idea at the beginning, but we decided that we can do better in terms of API.
+
-* Provide best performance out of the box. We have specialized libraries for that. Need only API wrappers for them.
+* Be like NumPy. It was the idea at the beginning, but we decided that we can do better in terms of API.
+* Provide the best performance out of the box. We have specialized libraries for that. Need only API wrappers for them.
 * Cover all cases as immediately and in one bundle. We will modularize everything and add new features gradually.
-* Provide specialized behavior in the core. API is made generic on purpose, so one needs to specialize for types, like for `Double` in the core. For that we will have specialization modules like `for-real`, which will give better experience for those, who want to work with specific types.
+* Provide specialized behavior in the core. API is made generic on purpose, so one needs to specialize for types, like 
+for `Double` in the core. For that we will have specialization modules like `for-real`, which will give better 
+experience for those, who want to work with specific types.
 
 ## Features
 
-Actual feature list is [here](/docs/features.md)
+Current feature list is [here](/docs/features.md)
 
 * **Algebra**
-    * Algebraic structures like rings, spaces and field (**TODO** add example to wiki)
+    * Algebraic structures like rings, spaces and fields (**TODO** add example to wiki)
     * Basic linear algebra operations (sums, products, etc.), backed by the `Space` API.
-    * Complex numbers backed by the `Field` API (meaning that they will be usable in any structure like vectors and N-dimensional arrays).
+    * Complex numbers backed by the `Field` API (meaning they will be usable in any structure like vectors and 
+    N-dimensional arrays).
     * Advanced linear algebra operations like matrix inversion and LU decomposition.
 
 * **Array-like structures** Full support of many-dimensional array-like structures 
 including mixed arithmetic operations and function operations over arrays and numbers (with the added benefit of static type checking).
 
-* **Expressions** By writing a single mathematical expression
+* **Expressions** By writing a single mathematical expression once, users will be able to apply different types of 
-once, users will be able to apply different types of objects to the expression by providing a context. Expressions
+objects to the expression by providing a context. Expressions can be used for a wide variety of purposes from high 
-can be used for a wide variety of purposes from high performance calculations to code generation.
+performance calculations to code generation.
 
 * **Histograms** Fast multi-dimensional histograms.
 
@@ -50,12 +59,11 @@ can be used for a wide variety of purposes from high performance calculations to
 
 * **Type-safe dimensions** Type-safe dimensions for matrix operations.
 
-* **Commons-math wrapper** It is planned to gradually wrap most parts of [Apache commons-math](http://commons.apache.org/proper/commons-math/)
+* **Commons-math wrapper** It is planned to gradually wrap most parts of 
-                           library in Kotlin code and maybe rewrite some parts to better suit the Kotlin programming paradigm, however there is no fixed roadmap for that. Feel free
+[Apache commons-math](http://commons.apache.org/proper/commons-math/) library in Kotlin code and maybe rewrite some 
-                           to submit a feature request if you want something to be done first.
+parts to better suit the Kotlin programming paradigm, however there is no established roadmap for that. Feel free to 
-
+submit a feature request if you want something to be implemented first.
-* **EJML wrapper** Provides EJML `SimpleMatrix` wrapper consistent with the core matrix structures.
+                           
-
 ## Planned features
 
 * **Messaging** A mathematical notation to support multi-language and multi-node communication for mathematical tasks.
@@ -101,7 +109,7 @@ can be used for a wide variety of purposes from high performance calculations to
 > - [buffers](kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : One-dimensional structure
 > - [expressions](kmath-core/src/commonMain/kotlin/kscience/kmath/expressions) : Functional Expressions
 > - [domains](kmath-core/src/commonMain/kotlin/kscience/kmath/domains) : Domains
-> - [autodif](kmath-core/src/commonMain/kotlin/kscience/kmath/misc/AutoDiff.kt) : Automatic differentiation
+> - [autodif](kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/SimpleAutoDiff.kt) : Automatic differentiation
 
 <hr/>
 
@@ -117,6 +125,12 @@ can be used for a wide variety of purposes from high performance calculations to
 > **Maturity**: EXPERIMENTAL
 <hr/>
 
+* ### [kmath-ejml](kmath-ejml)
+> 
+>
+> **Maturity**: EXPERIMENTAL
+<hr/>
+
 * ### [kmath-for-real](kmath-for-real)
 > 
 >
@@ -147,7 +161,19 @@ can be used for a wide variety of purposes from high performance calculations to
 > **Maturity**: EXPERIMENTAL
 <hr/>
 
-* ### [kmath-prob](kmath-prob)
+* ### [kmath-nd4j](kmath-nd4j)
+> ND4J NDStructure implementation and according NDAlgebra classes
+>
+> **Maturity**: EXPERIMENTAL
+>
+> **Features:**
+> - [nd4jarraystrucure](kmath-nd4j/src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt) : NDStructure wrapper for INDArray
+> - [nd4jarrayrings](kmath-nd4j/src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt) : Rings over Nd4jArrayStructure of Int and Long
+> - [nd4jarrayfields](kmath-nd4j/src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : Fields over Nd4jArrayStructure of Float and Double
+
+<hr/>
+
+* ### [kmath-stat](kmath-stat)
 > 
 >
 > **Maturity**: EXPERIMENTAL
@@ -162,39 +188,69 @@ can be used for a wide variety of purposes from high performance calculations to
 
 ## Multi-platform support
 
-KMath is developed as a multi-platform library, which means that most of the interfaces are declared in the [common module](/kmath-core/src/commonMain). Implementation is also done in the common module wherever possible. In some cases, features are delegated to platform-specific implementations even if they could be done in the common module for performance reasons. Currently, the JVM is the main focus of development, however Kotlin/Native and Kotlin/JS contributions are also welcome.
+KMath is developed as a multi-platform library, which means that most of the interfaces are declared in the 
+[common source sets](/kmath-core/src/commonMain) and implemented there wherever it is possible. In some cases, features 
+are delegated to platform-specific implementations even if they could be provided in the common module for performance 
+reasons. Currently, the Kotlin/JVM is the primary platform, however Kotlin/Native and Kotlin/JS contributions and 
+feedback are also welcome.
 
 ## Performance
 
-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 focus on creating convenient universal API first and then work on increasing performance for specific cases. We expect the worst KMath benchmarks will perform better than native Python, but worse than optimized native/SciPy (mostly due to boxing operations on primitive numbers). The best performance of optimized parts could be better than SciPy.
+Calculation performance is one of major goals of KMath in the future, but in some cases it is impossible to achieve 
+both performance and flexibility. 
 
-### Dependency
+We expect to focus on creating convenient universal API first and then work on increasing performance for specific 
+cases. We expect the worst KMath benchmarks will perform better than native Python, but worse than optimized 
+native/SciPy (mostly due to boxing operations on primitive numbers). The best performance of optimized parts could be 
+better than SciPy.
 
-Release artifacts are accessible from bintray with following configuration (see documentation for [kotlin-multiplatform](https://kotlinlang.org/docs/reference/multiplatform.html) form more details):
+### Repositories
+
+Release artifacts are accessible from bintray with following configuration (see documentation of 
+[Kotlin Multiplatform](https://kotlinlang.org/docs/reference/multiplatform.html) for more details):
 
 ```kotlin
-repositories{
+repositories {
+    jcenter()
+    maven("https://clojars.org/repo")
+    maven("https://dl.bintray.com/egor-bogomolov/astminer/")
+    maven("https://dl.bintray.com/hotkeytlt/maven")
+    maven("https://dl.bintray.com/kotlin/kotlin-eap")
+    maven("https://dl.bintray.com/kotlin/kotlinx")
     maven("https://dl.bintray.com/mipt-npm/kscience")
+    maven("https://jitpack.io")
+    mavenCentral()
 }
 
-dependencies{
+dependencies {
-    api("kscience.kmath:kmath-core:${kmathVersion}")
+    api("kscience.kmath:kmath-core:0.2.0-dev-3")
-    //api("scientifik:kmath-core:${kmathVersion}") for 0.1.3 and earlier
+    // api("kscience.kmath:kmath-core-jvm:0.2.0-dev-3") for jvm-specific version
 }
 ```
 
 Gradle `6.0+` is required for multiplatform artifacts.
 
-### Development
+#### Development
+
+Development builds are uploaded to the separate repository:
 
-Development builds are accessible from the reposirtory 
 ```kotlin
-repositories{
+repositories {
+    jcenter()
+    maven("https://clojars.org/repo")
+    maven("https://dl.bintray.com/egor-bogomolov/astminer/")
+    maven("https://dl.bintray.com/hotkeytlt/maven")
+    maven("https://dl.bintray.com/kotlin/kotlin-eap")
+    maven("https://dl.bintray.com/kotlin/kotlinx")
     maven("https://dl.bintray.com/mipt-npm/dev")
+    maven("https://jitpack.io")
+    mavenCentral()
 }
 ```
-with the same artifact names.
 
 ## Contributing
 
-The project requires a lot of additional work. Please feel free to contribute in any way and propose new features.
+The project requires a lot of additional work. The most important thing we need is a feedback about what features are 
+required the most. Feel free to create feature requests. We are also welcome to code contributions, 
+especially in issues marked with 
+[waiting for a hero](https://github.com/mipt-npm/kmath/labels/waiting%20for%20a%20hero) label.

build.gradle.kts

@@ -1,17 +1,26 @@
+import ru.mipt.npm.gradle.KSciencePublishPlugin
+
 plugins {
     id("ru.mipt.npm.project")
 }
 
-val kmathVersion: String by extra("0.2.0-dev-2")
+internal val kmathVersion: String by extra("0.2.0-dev-3")
-val bintrayRepo: String by extra("kscience")
+internal val bintrayRepo: String by extra("kscience")
-val githubProject: String by extra("kmath")
+internal val githubProject: String by extra("kmath")
 
 allprojects {
     repositories {
         jcenter()
+        maven("https://clojars.org/repo")
+        maven("https://dl.bintray.com/egor-bogomolov/astminer/")
+        maven("https://dl.bintray.com/hotkeytlt/maven")
         maven("https://dl.bintray.com/kotlin/kotlin-eap")
         maven("https://dl.bintray.com/kotlin/kotlinx")
-        maven("https://dl.bintray.com/hotkeytlt/maven")
+        maven("https://dl.bintray.com/mipt-npm/dev")
+        maven("https://dl.bintray.com/mipt-npm/kscience")
+        maven("https://jitpack.io")
+        maven("http://logicrunch.research.it.uu.se/maven/")
+        mavenCentral()
     }
 
     group = "kscience.kmath"
@@ -19,9 +28,13 @@ allprojects {
 }
 
 subprojects {
-    if (name.startsWith("kmath")) apply<ru.mipt.npm.gradle.KSciencePublishPlugin>()
+    if (name.startsWith("kmath")) apply<KSciencePublishPlugin>()
 }
 
 readme {
     readmeTemplate = file("docs/templates/README-TEMPLATE.md")
 }
+
+apiValidation {
+    validationDisabled = true
+}

docs/templates/README-TEMPLATE.md

@@ -8,41 +8,50 @@ Bintray:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience
 Bintray-dev:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-core/_latestVersion)
 
 # KMath
-Could be pronounced as `key-math`.
+
-The Kotlin MATHematics library  was initially intended as a Kotlin-based analog to Python's `numpy` library. Later we found that kotlin is much more flexible language and allows superior architecture designs. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. The `numpy`-like experience could be achieved with [kmath-for-real](/kmath-for-real) extension module.
+Could be pronounced as `key-math`. The Kotlin MATHematics library was initially intended as a Kotlin-based analog to 
+Python's NumPy library. Later we found that kotlin is much more flexible language and allows superior architecture 
+designs. In contrast to `numpy` and `scipy` it is modular and has a lightweight core. The `numpy`-like experience could 
+be achieved with [kmath-for-real](/kmath-for-real) extension module.
 
 ## Publications and talks
+
 * [A conceptual article about context-oriented design](https://proandroiddev.com/an-introduction-context-oriented-programming-in-kotlin-2e79d316b0a2)
 * [Another article about context-oriented design](https://proandroiddev.com/diving-deeper-into-context-oriented-programming-in-kotlin-3ecb4ec38814)
 * [ACAT 2019 conference paper](https://aip.scitation.org/doi/abs/10.1063/1.5130103)
 
 # Goal
-* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM and JS for now and Native in future). 
+
+* Provide a flexible and powerful API to work with mathematics abstractions in Kotlin-multiplatform (JVM, JS and Native). 
 * Provide basic multiplatform implementations for those abstractions (without significant performance optimization).
 * Provide bindings and wrappers with those abstractions for popular optimized platform libraries.
 
 ## Non-goals
-* Be like Numpy. It was the idea at the beginning, but we decided that we can do better in terms of API.
+
-* Provide best performance out of the box. We have specialized libraries for that. Need only API wrappers for them.
+* Be like NumPy. It was the idea at the beginning, but we decided that we can do better in terms of API.
+* Provide the best performance out of the box. We have specialized libraries for that. Need only API wrappers for them.
 * Cover all cases as immediately and in one bundle. We will modularize everything and add new features gradually.
-* Provide specialized behavior in the core. API is made generic on purpose, so one needs to specialize for types, like for `Double` in the core. For that we will have specialization modules like `for-real`, which will give better experience for those, who want to work with specific types.
+* Provide specialized behavior in the core. API is made generic on purpose, so one needs to specialize for types, like 
+for `Double` in the core. For that we will have specialization modules like `for-real`, which will give better 
+experience for those, who want to work with specific types.
 
 ## Features
 
-Actual feature list is [here](/docs/features.md)
+Current feature list is [here](/docs/features.md)
 
 * **Algebra**
-    * Algebraic structures like rings, spaces and field (**TODO** add example to wiki)
+    * Algebraic structures like rings, spaces and fields (**TODO** add example to wiki)
     * Basic linear algebra operations (sums, products, etc.), backed by the `Space` API.
-    * Complex numbers backed by the `Field` API (meaning that they will be usable in any structure like vectors and N-dimensional arrays).
+    * Complex numbers backed by the `Field` API (meaning they will be usable in any structure like vectors and 
+    N-dimensional arrays).
     * Advanced linear algebra operations like matrix inversion and LU decomposition.
 
 * **Array-like structures** Full support of many-dimensional array-like structures 
 including mixed arithmetic operations and function operations over arrays and numbers (with the added benefit of static type checking).
 
-* **Expressions** By writing a single mathematical expression
+* **Expressions** By writing a single mathematical expression once, users will be able to apply different types of 
-once, users will be able to apply different types of objects to the expression by providing a context. Expressions
+objects to the expression by providing a context. Expressions can be used for a wide variety of purposes from high 
-can be used for a wide variety of purposes from high performance calculations to code generation.
+performance calculations to code generation.
 
 * **Histograms** Fast multi-dimensional histograms.
 
@@ -50,9 +59,10 @@ can be used for a wide variety of purposes from high performance calculations to
 
 * **Type-safe dimensions** Type-safe dimensions for matrix operations.
 
-* **Commons-math wrapper** It is planned to gradually wrap most parts of [Apache commons-math](http://commons.apache.org/proper/commons-math/)
+* **Commons-math wrapper** It is planned to gradually wrap most parts of 
-                           library in Kotlin code and maybe rewrite some parts to better suit the Kotlin programming paradigm, however there is no fixed roadmap for that. Feel free
+[Apache commons-math](http://commons.apache.org/proper/commons-math/) library in Kotlin code and maybe rewrite some 
-                           to submit a feature request if you want something to be done first.
+parts to better suit the Kotlin programming paradigm, however there is no established roadmap for that. Feel free to 
+submit a feature request if you want something to be implemented first.
                            
 ## Planned features
 
@@ -72,39 +82,53 @@ $modules
 
 ## Multi-platform support
 
-KMath is developed as a multi-platform library, which means that most of the interfaces are declared in the [common module](/kmath-core/src/commonMain). Implementation is also done in the common module wherever possible. In some cases, features are delegated to platform-specific implementations even if they could be done in the common module for performance reasons. Currently, the JVM is the main focus of development, however Kotlin/Native and Kotlin/JS contributions are also welcome.
+KMath is developed as a multi-platform library, which means that most of the interfaces are declared in the 
+[common source sets](/kmath-core/src/commonMain) and implemented there wherever it is possible. In some cases, features 
+are delegated to platform-specific implementations even if they could be provided in the common module for performance 
+reasons. Currently, the Kotlin/JVM is the primary platform, however Kotlin/Native and Kotlin/JS contributions and 
+feedback are also welcome.
 
 ## Performance
 
-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 focus on creating convenient universal API first and then work on increasing performance for specific cases. We expect the worst KMath benchmarks will perform better than native Python, but worse than optimized native/SciPy (mostly due to boxing operations on primitive numbers). The best performance of optimized parts could be better than SciPy.
+Calculation performance is one of major goals of KMath in the future, but in some cases it is impossible to achieve 
+both performance and flexibility. 
 
-### Dependency
+We expect to focus on creating convenient universal API first and then work on increasing performance for specific 
+cases. We expect the worst KMath benchmarks will perform better than native Python, but worse than optimized 
+native/SciPy (mostly due to boxing operations on primitive numbers). The best performance of optimized parts could be 
+better than SciPy.
 
-Release artifacts are accessible from bintray with following configuration (see documentation for [kotlin-multiplatform](https://kotlinlang.org/docs/reference/multiplatform.html) form more details):
+### Repositories
+
+Release artifacts are accessible from bintray with following configuration (see documentation of 
+[Kotlin Multiplatform](https://kotlinlang.org/docs/reference/multiplatform.html) for more details):
 
 ```kotlin
-repositories{
+repositories {
     maven("https://dl.bintray.com/mipt-npm/kscience")
 }
 
-dependencies{
+dependencies {
     api("kscience.kmath:kmath-core:$version")
-    //api("kscience.kmath:kmath-core-jvm:$version") for jvm-specific version
+    // api("kscience.kmath:kmath-core-jvm:$version") for jvm-specific version
 }
 ```
 
 Gradle `6.0+` is required for multiplatform artifacts.
 
-### Development
+#### Development
+
+Development builds are uploaded to the separate repository:
 
-Development builds are accessible from the reposirtory 
 ```kotlin
-repositories{
+repositories {
     maven("https://dl.bintray.com/mipt-npm/dev")
 }
 ```
-with the same artifact names.
 
 ## Contributing
 
-The project requires a lot of additional work. Please feel free to contribute in any way and propose new features.
+The project requires a lot of additional work. The most important thing we need is a feedback about what features are 
+required the most. Feel free to create feature requests. We are also welcome to code contributions, 
+especially in issues marked with 
+[waiting for a hero](https://github.com/mipt-npm/kmath/labels/waiting%20for%20a%20hero) label.

examples/build.gradle.kts

@@ -8,25 +8,46 @@ plugins {
 }
 
 allOpen.annotation("org.openjdk.jmh.annotations.State")
+sourceSets.register("benchmarks")
 
 repositories {
-    maven("https://dl.bintray.com/mipt-npm/kscience")
+    jcenter()
+    maven("https://clojars.org/repo")
+    maven("https://dl.bintray.com/egor-bogomolov/astminer/")
+    maven("https://dl.bintray.com/hotkeytlt/maven")
+    maven("https://dl.bintray.com/kotlin/kotlin-eap")
+    maven("https://dl.bintray.com/kotlin/kotlinx")
     maven("https://dl.bintray.com/mipt-npm/dev")
-    maven("https://dl.bintray.com/kotlin/kotlin-dev/")
+    maven("https://dl.bintray.com/mipt-npm/kscience")
+    maven("https://jitpack.io")
+    maven("http://logicrunch.research.it.uu.se/maven/")
     mavenCentral()
 }
 
-sourceSets.register("benchmarks")
-
 dependencies {
-//    implementation(project(":kmath-ast"))
+    implementation(project(":kmath-ast"))
+    implementation(project(":kmath-kotlingrad"))
     implementation(project(":kmath-core"))
     implementation(project(":kmath-coroutines"))
     implementation(project(":kmath-commons"))
-    implementation(project(":kmath-prob"))
+    implementation(project(":kmath-stat"))
     implementation(project(":kmath-viktor"))
     implementation(project(":kmath-dimensions"))
     implementation(project(":kmath-ejml"))
+    implementation(project(":kmath-nd4j"))
+    implementation("org.deeplearning4j:deeplearning4j-core:1.0.0-beta7")
+    implementation("org.nd4j:nd4j-native:1.0.0-beta7")
+
+//    uncomment if your system supports AVX2
+//    val os = System.getProperty("os.name")
+//
+//    if (System.getProperty("os.arch") in arrayOf("x86_64", "amd64")) when {
+//        os.startsWith("Windows") -> implementation("org.nd4j:nd4j-native:1.0.0-beta7:windows-x86_64-avx2")
+//        os == "Linux" -> implementation("org.nd4j:nd4j-native:1.0.0-beta7:linux-x86_64-avx2")
+//        os == "Mac OS X" -> implementation("org.nd4j:nd4j-native:1.0.0-beta7:macosx-x86_64-avx2")
+//    } else
+    implementation("org.nd4j:nd4j-native-platform:1.0.0-beta7")
+
     implementation("org.jetbrains.kotlinx:kotlinx-io:0.2.0-npm-dev-11")
     implementation("org.jetbrains.kotlinx:kotlinx.benchmark.runtime:0.2.0-dev-20")
     implementation("org.slf4j:slf4j-simple:1.7.30")
@@ -55,4 +76,6 @@ kotlin.sourceSets.all {
     }
 }
 
-tasks.withType<KotlinCompile> { kotlinOptions.jvmTarget = "11" }
+tasks.withType<KotlinCompile> {
+    kotlinOptions.jvmTarget = "11"
+}

examples/src/main/kotlin/kscience/kmath/ast/ExpressionsInterpretersBenchmark.kt

@@ -1,70 +1,80 @@
 package kscience.kmath.ast
-//
+
-//import kscience.kmath.asm.compile
+import kscience.kmath.asm.compile
-//import kscience.kmath.expressions.Expression
+import kscience.kmath.expressions.Expression
-//import kscience.kmath.expressions.expressionInField
+import kscience.kmath.expressions.expressionInField
-//import kscience.kmath.expressions.invoke
+import kscience.kmath.expressions.invoke
-//import kscience.kmath.operations.Field
+import kscience.kmath.operations.Field
-//import kscience.kmath.operations.RealField
+import kscience.kmath.operations.RealField
-//import kotlin.random.Random
+import kotlin.random.Random
-//import kotlin.system.measureTimeMillis
+import kotlin.system.measureTimeMillis
-//
+
-//class ExpressionsInterpretersBenchmark {
+internal class ExpressionsInterpretersBenchmark {
-//    private val algebra: Field<Double> = RealField
+    private val algebra: Field<Double> = RealField
-//    fun functionalExpression() {
+    fun functionalExpression() {
-//        val expr = algebra.expressionInField {
+        val expr = algebra.expressionInField {
-//            variable("x") * const(2.0) + const(2.0) / variable("x") - const(16.0)
+            symbol("x") * const(2.0) + const(2.0) / symbol("x") - const(16.0)
-//        }
+        }
-//
+
-//        invokeAndSum(expr)
+        invokeAndSum(expr)
-//    }
+    }
-//
+
-//    fun mstExpression() {
+    fun mstExpression() {
-//        val expr = algebra.mstInField {
+        val expr = algebra.mstInField {
-//            symbol("x") * number(2.0) + number(2.0) / symbol("x") - number(16.0)
+            symbol("x") * number(2.0) + number(2.0) / symbol("x") - number(16.0)
-//        }
+        }
-//
+
-//        invokeAndSum(expr)
+        invokeAndSum(expr)
-//    }
+    }
-//
+
-//    fun asmExpression() {
+    fun asmExpression() {
-//        val expr = algebra.mstInField {
+        val expr = algebra.mstInField {
-//            symbol("x") * number(2.0) + number(2.0) / symbol("x") - number(16.0)
+            symbol("x") * number(2.0) + number(2.0) / symbol("x") - number(16.0)
-//        }.compile()
+        }.compile()
-//
+
-//        invokeAndSum(expr)
+        invokeAndSum(expr)
-//    }
+    }
-//
+
-//    private fun invokeAndSum(expr: Expression<Double>) {
+    private fun invokeAndSum(expr: Expression<Double>) {
-//        val random = Random(0)
+        val random = Random(0)
-//        var sum = 0.0
+        var sum = 0.0
-//
+
-//        repeat(1000000) {
+        repeat(1000000) {
-//            sum += expr("x" to random.nextDouble())
+            sum += expr("x" to random.nextDouble())
-//        }
+        }
-//
+
-//        println(sum)
+        println(sum)
-//    }
+    }
-//}
+}
-//
+
-//fun main() {
+/**
-//    val benchmark = ExpressionsInterpretersBenchmark()
+ * This benchmark compares basically evaluation of simple function with MstExpression interpreter, ASM backend and
-//
+ * core FunctionalExpressions API.
-//    val fe = measureTimeMillis {
+ *
-//        benchmark.functionalExpression()
+ * The expected rating is:
-//    }
+ *
-//
+ * 1. ASM.
-//    println("fe=$fe")
+ * 2. MST.
-//
+ * 3. FE.
-//    val mst = measureTimeMillis {
+ */
-//        benchmark.mstExpression()
+fun main() {
-//    }
+    val benchmark = ExpressionsInterpretersBenchmark()
-//
+
-//    println("mst=$mst")
+    val fe = measureTimeMillis {
-//
+        benchmark.functionalExpression()
-//    val asm = measureTimeMillis {
+    }
-//        benchmark.asmExpression()
+
-//    }
+    println("fe=$fe")
-//
+
-//    println("asm=$asm")
+    val mst = measureTimeMillis {
-//}
+        benchmark.mstExpression()
+    }
+
+    println("mst=$mst")
+
+    val asm = measureTimeMillis {
+        benchmark.asmExpression()
+    }
+
+    println("asm=$asm")
+}

examples/src/main/kotlin/kscience/kmath/ast/KotlingradSupport.kt

@@ -0,0 +1,24 @@
+package kscience.kmath.ast
+
+import kscience.kmath.asm.compile
+import kscience.kmath.expressions.derivative
+import kscience.kmath.expressions.invoke
+import kscience.kmath.expressions.symbol
+import kscience.kmath.kotlingrad.differentiable
+import kscience.kmath.operations.RealField
+
+/**
+ * In this example, x^2-4*x-44 function is differentiated with Kotlin∇, and the autodiff result is compared with
+ * valid derivative.
+ */
+fun main() {
+    val x by symbol
+
+    val actualDerivative = MstExpression(RealField, "x^2-4*x-44".parseMath())
+        .differentiable()
+        .derivative(x)
+        .compile()
+
+    val expectedDerivative = MstExpression(RealField, "2*x-4".parseMath()).compile()
+    assert(actualDerivative("x" to 123.0) == expectedDerivative("x" to 123.0))
+}

examples/src/main/kotlin/kscience/kmath/operations/ComplexDemo.kt

@@ -6,8 +6,8 @@ import kscience.kmath.structures.complex
 
 fun main() {
     // 2d element
-    val element = NDElement.complex(2, 2) { index: IntArray ->
+    val element = NDElement.complex(2, 2) { (i,j) ->
-        Complex(index[0].toDouble() - index[1].toDouble(), index[0].toDouble() + index[1].toDouble())
+        Complex(i.toDouble() - j.toDouble(), i.toDouble() + j.toDouble())
     }
     println(element)
 

examples/src/main/kotlin/kscience/kmath/structures/NDField.kt

@@ -1,8 +1,10 @@
 package kscience.kmath.structures
 
 import kotlinx.coroutines.GlobalScope
+import kscience.kmath.nd4j.Nd4jArrayField
 import kscience.kmath.operations.RealField
 import kscience.kmath.operations.invoke
+import org.nd4j.linalg.factory.Nd4j
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract
 import kotlin.system.measureTimeMillis
@@ -14,6 +16,8 @@ internal inline fun measureAndPrint(title: String, block: () -> Unit) {
 }
 
 fun main() {
+    // initializing Nd4j
+    Nd4j.zeros(0)
     val dim = 1000
     val n = 1000
 
@@ -23,6 +27,8 @@ fun main() {
     val specializedField = NDField.real(dim, dim)
     //A generic boxing field. It should be used for objects, not primitives.
     val genericField = NDField.boxing(RealField, dim, dim)
+    // Nd4j specialized field.
+    val nd4jField = Nd4jArrayField.real(dim, dim)
 
     measureAndPrint("Automatic field addition") {
         autoField {
@@ -43,6 +49,13 @@ fun main() {
         }
     }
 
+    measureAndPrint("Nd4j specialized addition") {
+        nd4jField {
+            var res = one
+            repeat(n) { res += 1.0 as Number }
+        }
+    }
+
     measureAndPrint("Lazy addition") {
         val res = specializedField.one.mapAsync(GlobalScope) {
             var c = 0.0

kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MstAlgebra.kt

@@ -6,14 +6,14 @@ import kscience.kmath.operations.*
  * [Algebra] over [MST] nodes.
  */
 public object MstAlgebra : NumericAlgebra<MST> {
-    override fun number(value: Number): MST = MST.Numeric(value)
+    override fun number(value: Number): MST.Numeric = MST.Numeric(value)
 
-    override fun symbol(value: String): MST = MST.Symbolic(value)
+    override fun symbol(value: String): MST.Symbolic = MST.Symbolic(value)
 
-    override fun unaryOperation(operation: String, arg: MST): MST =
+    override fun unaryOperation(operation: String, arg: MST): MST.Unary =
         MST.Unary(operation, arg)
 
-    override fun binaryOperation(operation: String, left: MST, right: MST): MST =
+    override fun binaryOperation(operation: String, left: MST, right: MST): MST.Binary =
         MST.Binary(operation, left, right)
 }
 
@@ -21,97 +21,100 @@ public object MstAlgebra : NumericAlgebra<MST> {
  * [Space] over [MST] nodes.
  */
 public object MstSpace : Space<MST>, NumericAlgebra<MST> {
-    override val zero: MST = number(0.0)
+    override val zero: MST.Numeric by lazy { number(0.0) }
 
-    override fun number(value: Number): MST = MstAlgebra.number(value)
+    override fun number(value: Number): MST.Numeric = MstAlgebra.number(value)
-    override fun symbol(value: String): MST = MstAlgebra.symbol(value)
+    override fun symbol(value: String): MST.Symbolic = MstAlgebra.symbol(value)
-    override fun add(a: MST, b: MST): MST = binaryOperation(SpaceOperations.PLUS_OPERATION, a, b)
+    override fun add(a: MST, b: MST): MST.Binary = binaryOperation(SpaceOperations.PLUS_OPERATION, a, b)
-    override fun multiply(a: MST, k: Number): MST = binaryOperation(RingOperations.TIMES_OPERATION, a, number(k))
+    override fun multiply(a: MST, k: Number): MST.Binary = binaryOperation(RingOperations.TIMES_OPERATION, a, number(k))
 
-    override fun binaryOperation(operation: String, left: MST, right: MST): MST =
+    override fun binaryOperation(operation: String, left: MST, right: MST): MST.Binary =
         MstAlgebra.binaryOperation(operation, left, right)
 
-    override fun unaryOperation(operation: String, arg: MST): MST = MstAlgebra.unaryOperation(operation, arg)
+    override fun unaryOperation(operation: String, arg: MST): MST.Unary = MstAlgebra.unaryOperation(operation, arg)
 }
 
 /**
  * [Ring] over [MST] nodes.
  */
 public object MstRing : Ring<MST>, NumericAlgebra<MST> {
-    override val zero: MST
+    override val zero: MST.Numeric
         get() = MstSpace.zero
-    override val one: MST = number(1.0)
 
-    override fun number(value: Number): MST = MstSpace.number(value)
+    override val one: MST.Numeric by lazy { number(1.0) }
-    override fun symbol(value: String): MST = MstSpace.symbol(value)
-    override fun add(a: MST, b: MST): MST = MstSpace.add(a, b)
 
-    override fun multiply(a: MST, k: Number): MST = MstSpace.multiply(a, k)
+    override fun number(value: Number): MST.Numeric = MstSpace.number(value)
+    override fun symbol(value: String): MST.Symbolic = MstSpace.symbol(value)
+    override fun add(a: MST, b: MST): MST.Binary = MstSpace.add(a, b)
+    override fun multiply(a: MST, k: Number): MST.Binary = MstSpace.multiply(a, k)
+    override fun multiply(a: MST, b: MST): MST.Binary = binaryOperation(RingOperations.TIMES_OPERATION, a, b)
 
-    override fun multiply(a: MST, b: MST): MST = binaryOperation(RingOperations.TIMES_OPERATION, a, b)
+    override fun binaryOperation(operation: String, left: MST, right: MST): MST.Binary =
-
-    override fun binaryOperation(operation: String, left: MST, right: MST): MST =
         MstSpace.binaryOperation(operation, left, right)
 
-    override fun unaryOperation(operation: String, arg: MST): MST = MstAlgebra.unaryOperation(operation, arg)
+    override fun unaryOperation(operation: String, arg: MST): MST.Unary = MstSpace.unaryOperation(operation, arg)
 }
 
 /**
  * [Field] over [MST] nodes.
  */
 public object MstField : Field<MST> {
-    public override val zero: MST
+    public override val zero: MST.Numeric
         get() = MstRing.zero
 
-    public override val one: MST
+    public override val one: MST.Numeric
         get() = MstRing.one
 
-    public override fun symbol(value: String): MST = MstRing.symbol(value)
+    public override fun symbol(value: String): MST.Symbolic = MstRing.symbol(value)
-    public override fun number(value: Number): MST = MstRing.number(value)
+    public override fun number(value: Number): MST.Numeric = MstRing.number(value)
-    public override fun add(a: MST, b: MST): MST = MstRing.add(a, b)
+    public override fun add(a: MST, b: MST): MST.Binary = MstRing.add(a, b)
-    public override fun multiply(a: MST, k: Number): MST = MstRing.multiply(a, k)
+    public override fun multiply(a: MST, k: Number): MST.Binary = MstRing.multiply(a, k)
-    public override fun multiply(a: MST, b: MST): MST = MstRing.multiply(a, b)
+    public override fun multiply(a: MST, b: MST): MST.Binary = MstRing.multiply(a, b)
-    public override fun divide(a: MST, b: MST): MST = binaryOperation(FieldOperations.DIV_OPERATION, a, b)
+    public override fun divide(a: MST, b: MST): MST.Binary = binaryOperation(FieldOperations.DIV_OPERATION, a, b)
 
-    public override fun binaryOperation(operation: String, left: MST, right: MST): MST =
+    public override fun binaryOperation(operation: String, left: MST, right: MST): MST.Binary =
         MstRing.binaryOperation(operation, left, right)
 
-    override fun unaryOperation(operation: String, arg: MST): MST = MstRing.unaryOperation(operation, arg)
+    override fun unaryOperation(operation: String, arg: MST): MST.Unary = MstRing.unaryOperation(operation, arg)
 }
 
 /**
  * [ExtendedField] over [MST] nodes.
  */
 public object MstExtendedField : ExtendedField<MST> {
-    override val zero: MST
+    override val zero: MST.Numeric
         get() = MstField.zero
 
-    override val one: MST
+    override val one: MST.Numeric
         get() = MstField.one
 
-    override fun symbol(value: String): MST = MstField.symbol(value)
+    override fun symbol(value: String): MST.Symbolic = MstField.symbol(value)
-    override fun sin(arg: MST): MST = unaryOperation(TrigonometricOperations.SIN_OPERATION, arg)
+    override fun number(value: Number): MST.Numeric = MstField.number(value)
-    override fun cos(arg: MST): MST = unaryOperation(TrigonometricOperations.COS_OPERATION, arg)
+    override fun sin(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.SIN_OPERATION, arg)
-    override fun tan(arg: MST): MST = unaryOperation(TrigonometricOperations.TAN_OPERATION, arg)
+    override fun cos(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.COS_OPERATION, arg)
-    override fun asin(arg: MST): MST = unaryOperation(TrigonometricOperations.ASIN_OPERATION, arg)
+    override fun tan(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.TAN_OPERATION, arg)
-    override fun acos(arg: MST): MST = unaryOperation(TrigonometricOperations.ACOS_OPERATION, arg)
+    override fun asin(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.ASIN_OPERATION, arg)
-    override fun atan(arg: MST): MST = unaryOperation(TrigonometricOperations.ATAN_OPERATION, arg)
+    override fun acos(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.ACOS_OPERATION, arg)
-    override fun sinh(arg: MST): MST = unaryOperation(HyperbolicOperations.SINH_OPERATION, arg)
+    override fun atan(arg: MST): MST.Unary = unaryOperation(TrigonometricOperations.ATAN_OPERATION, arg)
-    override fun cosh(arg: MST): MST = unaryOperation(HyperbolicOperations.COSH_OPERATION, arg)
+    override fun sinh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.SINH_OPERATION, arg)
-    override fun tanh(arg: MST): MST = unaryOperation(HyperbolicOperations.TANH_OPERATION, arg)
+    override fun cosh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.COSH_OPERATION, arg)
-    override fun asinh(arg: MST): MST = unaryOperation(HyperbolicOperations.ASINH_OPERATION, arg)
+    override fun tanh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.TANH_OPERATION, arg)
-    override fun acosh(arg: MST): MST = unaryOperation(HyperbolicOperations.ACOSH_OPERATION, arg)
+    override fun asinh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.ASINH_OPERATION, arg)
-    override fun atanh(arg: MST): MST = unaryOperation(HyperbolicOperations.ATANH_OPERATION, arg)
+    override fun acosh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.ACOSH_OPERATION, arg)
-    override fun add(a: MST, b: MST): MST = MstField.add(a, b)
+    override fun atanh(arg: MST): MST.Unary = unaryOperation(HyperbolicOperations.ATANH_OPERATION, arg)
-    override fun multiply(a: MST, k: Number): MST = MstField.multiply(a, k)
+    override fun add(a: MST, b: MST): MST.Binary = MstField.add(a, b)
-    override fun multiply(a: MST, b: MST): MST = MstField.multiply(a, b)
+    override fun multiply(a: MST, k: Number): MST.Binary = MstField.multiply(a, k)
-    override fun divide(a: MST, b: MST): MST = MstField.divide(a, b)
+    override fun multiply(a: MST, b: MST): MST.Binary = MstField.multiply(a, b)
-    override fun power(arg: MST, pow: Number): MST = binaryOperation(PowerOperations.POW_OPERATION, arg, number(pow))
+    override fun divide(a: MST, b: MST): MST.Binary = MstField.divide(a, b)
-    override fun exp(arg: MST): MST = unaryOperation(ExponentialOperations.EXP_OPERATION, arg)
-    override fun ln(arg: MST): MST = unaryOperation(ExponentialOperations.LN_OPERATION, arg)
 
-    override fun binaryOperation(operation: String, left: MST, right: MST): MST =
+    override fun power(arg: MST, pow: Number): MST.Binary =
+        binaryOperation(PowerOperations.POW_OPERATION, arg, number(pow))
+
+    override fun exp(arg: MST): MST.Unary = unaryOperation(ExponentialOperations.EXP_OPERATION, arg)
+    override fun ln(arg: MST): MST.Unary = unaryOperation(ExponentialOperations.LN_OPERATION, arg)
+
+    override fun binaryOperation(operation: String, left: MST, right: MST): MST.Binary =
         MstField.binaryOperation(operation, left, right)
 
-    override fun unaryOperation(operation: String, arg: MST): MST = MstField.unaryOperation(operation, arg)
+    override fun unaryOperation(operation: String, arg: MST): MST.Unary = MstField.unaryOperation(operation, arg)
 }

kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MstExpression.kt

@@ -13,21 +13,22 @@ import kotlin.contracts.contract
  * @property mst the [MST] node.
  * @author Alexander Nozik
  */
-public class MstExpression<T>(public val algebra: Algebra<T>, public val mst: MST) : Expression<T> {
+public class MstExpression<T, out A : Algebra<T>>(public val algebra: A, public val mst: MST) : Expression<T> {
-    private inner class InnerAlgebra(val arguments: Map<String, T>) : NumericAlgebra<T> {
+    private inner class InnerAlgebra(val arguments: Map<Symbol, T>) : NumericAlgebra<T> {
-        override fun symbol(value: String): T = arguments[value] ?: algebra.symbol(value)
+        override fun symbol(value: String): T = arguments[StringSymbol(value)] ?: algebra.symbol(value)
         override fun unaryOperation(operation: String, arg: T): T = algebra.unaryOperation(operation, arg)
 
         override fun binaryOperation(operation: String, left: T, right: T): T =
             algebra.binaryOperation(operation, left, right)
 
-        override fun number(value: Number): T = if (algebra is NumericAlgebra)
+        @Suppress("UNCHECKED_CAST")
-            algebra.number(value)
+        override fun number(value: Number): T = if (algebra is NumericAlgebra<*>)
+            (algebra as NumericAlgebra<T>).number(value)
         else
             error("Numeric nodes are not supported by $this")
     }
 
-    override operator fun invoke(arguments: Map<String, T>): T = InnerAlgebra(arguments).evaluate(mst)
+    override operator fun invoke(arguments: Map<Symbol, T>): T = InnerAlgebra(arguments).evaluate(mst)
 }
 
 /**
@@ -37,15 +38,15 @@ public class MstExpression<T>(public val algebra: Algebra<T>, public val mst: MS
  */
 public inline fun <reified T : Any, A : Algebra<T>, E : Algebra<MST>> A.mst(
     mstAlgebra: E,
-    block: E.() -> MST
+    block: E.() -> MST,
-): MstExpression<T> = MstExpression(this, mstAlgebra.block())
+): MstExpression<T, A> = MstExpression(this, mstAlgebra.block())
 
 /**
  * Builds [MstExpression] over [Space].
  *
  * @author Alexander Nozik
  */
-public inline fun <reified T : Any> Space<T>.mstInSpace(block: MstSpace.() -> MST): MstExpression<T> {
+public inline fun <reified T : Any, A : Space<T>> A.mstInSpace(block: MstSpace.() -> MST): MstExpression<T, A> {
     contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
     return MstExpression(this, MstSpace.block())
 }
@@ -55,7 +56,7 @@ public inline fun <reified T : Any> Space<T>.mstInSpace(block: MstSpace.() -> MS
  *
  * @author Alexander Nozik
  */
-public inline fun <reified T : Any> Ring<T>.mstInRing(block: MstRing.() -> MST): MstExpression<T> {
+public inline fun <reified T : Any, A : Ring<T>> A.mstInRing(block: MstRing.() -> MST): MstExpression<T, A> {
     contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
     return MstExpression(this, MstRing.block())
 }
@@ -65,7 +66,7 @@ public inline fun <reified T : Any> Ring<T>.mstInRing(block: MstRing.() -> MST):
  *
  * @author Alexander Nozik
  */
-public inline fun <reified T : Any> Field<T>.mstInField(block: MstField.() -> MST): MstExpression<T> {
+public inline fun <reified T : Any, A : Field<T>> A.mstInField(block: MstField.() -> MST): MstExpression<T, A> {
     contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
     return MstExpression(this, MstField.block())
 }
@@ -75,7 +76,7 @@ public inline fun <reified T : Any> Field<T>.mstInField(block: MstField.() -> MS
  *
  * @author Iaroslav Postovalov
  */
-public inline fun <reified T : Any> Field<T>.mstInExtendedField(block: MstExtendedField.() -> MST): MstExpression<T> {
+public inline fun <reified T : Any, A : ExtendedField<T>> A.mstInExtendedField(block: MstExtendedField.() -> MST): MstExpression<T, A> {
     contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
     return MstExpression(this, MstExtendedField.block())
 }
@@ -85,7 +86,7 @@ public inline fun <reified T : Any> Field<T>.mstInExtendedField(block: MstExtend
  *
  * @author Alexander Nozik
  */
-public inline fun <reified T : Any, A : Space<T>> FunctionalExpressionSpace<T, A>.mstInSpace(block: MstSpace.() -> MST): MstExpression<T> {
+public inline fun <reified T : Any, A : Space<T>> FunctionalExpressionSpace<T, A>.mstInSpace(block: MstSpace.() -> MST): MstExpression<T, A> {
     contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
     return algebra.mstInSpace(block)
 }
@@ -95,7 +96,7 @@ public inline fun <reified T : Any, A : Space<T>> FunctionalExpressionSpace<T, A
  *
  * @author Alexander Nozik
  */
-public inline fun <reified T : Any, A : Ring<T>> FunctionalExpressionRing<T, A>.mstInRing(block: MstRing.() -> MST): MstExpression<T> {
+public inline fun <reified T : Any, A : Ring<T>> FunctionalExpressionRing<T, A>.mstInRing(block: MstRing.() -> MST): MstExpression<T, A> {
     contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
     return algebra.mstInRing(block)
 }
@@ -105,7 +106,7 @@ public inline fun <reified T : Any, A : Ring<T>> FunctionalExpressionRing<T, A>.
  *
  * @author Alexander Nozik
  */
-public inline fun <reified T : Any, A : Field<T>> FunctionalExpressionField<T, A>.mstInField(block: MstField.() -> MST): MstExpression<T> {
+public inline fun <reified T : Any, A : Field<T>> FunctionalExpressionField<T, A>.mstInField(block: MstField.() -> MST): MstExpression<T, A> {
     contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
     return algebra.mstInField(block)
 }
@@ -116,8 +117,8 @@ public inline fun <reified T : Any, A : Field<T>> FunctionalExpressionField<T, A
  * @author Iaroslav Postovalov
  */
 public inline fun <reified T : Any, A : ExtendedField<T>> FunctionalExpressionExtendedField<T, A>.mstInExtendedField(
-    block: MstExtendedField.() -> MST
+    block: MstExtendedField.() -> MST,
-): MstExpression<T> {
+): MstExpression<T, A> {
     contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
     return algebra.mstInExtendedField(block)
 }

kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/asm.kt

@@ -69,4 +69,5 @@ public inline fun <reified T : Any> Algebra<T>.expression(mst: MST): Expression<
  *
  * @author Alexander Nozik.
  */
-public inline fun <reified T : Any> MstExpression<T>.compile(): Expression<T> = mst.compileWith(T::class.java, algebra)
+public inline fun <reified T : Any> MstExpression<T, Algebra<T>>.compile(): Expression<T> =
+    mst.compileWith(T::class.java, algebra)

kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/internal/AsmBuilder.kt

@@ -25,7 +25,7 @@ internal class AsmBuilder<T> internal constructor(
     private val classOfT: Class<*>,
     private val algebra: Algebra<T>,
     private val className: String,
-    private val invokeLabel0Visitor: AsmBuilder<T>.() -> Unit
+    private val invokeLabel0Visitor: AsmBuilder<T>.() -> Unit,
 ) {
     /**
      * Internal classloader of [AsmBuilder] with alias to define class from byte array.
@@ -379,22 +379,14 @@ internal class AsmBuilder<T> internal constructor(
      * Loads a variable [name] from arguments [Map] parameter of [Expression.invoke]. The [defaultValue] may be
      * provided.
      */
-    internal fun loadVariable(name: String, defaultValue: T? = null): Unit = invokeMethodVisitor.run {
+    internal fun loadVariable(name: String): Unit = invokeMethodVisitor.run {
         load(invokeArgumentsVar, MAP_TYPE)
         aconst(name)
 
-        if (defaultValue != null)
-            loadTConstant(defaultValue)
-
         invokestatic(
             MAP_INTRINSICS_TYPE.internalName,
             "getOrFail",
-
+            Type.getMethodDescriptor(OBJECT_TYPE, MAP_TYPE, STRING_TYPE),
-            Type.getMethodDescriptor(
-                OBJECT_TYPE,
-                MAP_TYPE,
-                OBJECT_TYPE,
-                *OBJECT_TYPE.wrapToArrayIf { defaultValue != null }),
             false
         )
 
@@ -429,7 +421,7 @@ internal class AsmBuilder<T> internal constructor(
         method: String,
         descriptor: String,
         expectedArity: Int,
-        opcode: Int = INVOKEINTERFACE
+        opcode: Int = INVOKEINTERFACE,
     ) {
         run loop@{
             repeat(expectedArity) {

kmath-ast/src/jvmMain/kotlin/kscience/kmath/asm/internal/mapIntrinsics.kt

@@ -2,11 +2,12 @@
 
 package kscience.kmath.asm.internal
 
+import kscience.kmath.expressions.StringSymbol
+import kscience.kmath.expressions.Symbol
+
 /**
- * Gets value with given [key] or throws [IllegalStateException] whenever it is not present.
+ * Gets value with given [key] or throws [NoSuchElementException] whenever it is not present.
  *
  * @author Iaroslav Postovalov
  */
-@JvmOverloads
+internal fun <V> Map<Symbol, V>.getOrFail(key: String): V = getValue(StringSymbol(key))
-internal fun <K, V> Map<K, V>.getOrFail(key: K, default: V? = null): V =
-    this[key] ?: default ?: error("Parameter not found: $key")

kmath-ast/src/jvmTest/kotlin/kscience/kmath/asm/TestAsmAlgebras.kt

@@ -1,6 +1,5 @@
 package kscience.kmath.asm
 
-import kscience.kmath.asm.compile
 import kscience.kmath.ast.mstInField
 import kscience.kmath.ast.mstInRing
 import kscience.kmath.ast.mstInSpace
@@ -11,6 +10,7 @@ import kotlin.test.Test
 import kotlin.test.assertEquals
 
 internal class TestAsmAlgebras {
+
     @Test
     fun space() {
         val res1 = ByteRing.mstInSpace {

kmath-commons/build.gradle.kts

@@ -6,7 +6,7 @@ description = "Commons math binding for kmath"
 dependencies {
     api(project(":kmath-core"))
     api(project(":kmath-coroutines"))
-    api(project(":kmath-prob"))
+    api(project(":kmath-stat"))
-//    api(project(":kmath-functions"))
+    api(project(":kmath-functions"))
     api("org.apache.commons:commons-math3:3.6.1")
 }

kmath-commons/src/main/kotlin/kscience/kmath/commons/expressions/DerivativeStructureExpression.kt

@@ -0,0 +1,121 @@
+package kscience.kmath.commons.expressions
+
+import kscience.kmath.expressions.*
+import kscience.kmath.operations.ExtendedField
+import org.apache.commons.math3.analysis.differentiation.DerivativeStructure
+
+/**
+ * A field over commons-math [DerivativeStructure].
+ *
+ * @property order The derivation order.
+ * @property bindings The map of bindings values. All bindings are considered free parameters
+ */
+public class DerivativeStructureField(
+    public val order: Int,
+    bindings: Map<Symbol, Double>,
+) : ExtendedField<DerivativeStructure>, ExpressionAlgebra<Double, DerivativeStructure> {
+    public val numberOfVariables: Int = bindings.size
+
+    public override val zero: DerivativeStructure by lazy { DerivativeStructure(numberOfVariables, order) }
+    public override val one: DerivativeStructure by lazy { DerivativeStructure(numberOfVariables, order, 1.0) }
+
+    /**
+     * A class that implements both [DerivativeStructure] and a [Symbol]
+     */
+    public inner class DerivativeStructureSymbol(
+        size: Int,
+        index: Int,
+        symbol: Symbol,
+        value: Double,
+    ) : DerivativeStructure(size, order, index, value), Symbol {
+        override val identity: String = symbol.identity
+        override fun toString(): String = identity
+        override fun equals(other: Any?): Boolean = this.identity == (other as? Symbol)?.identity
+        override fun hashCode(): Int = identity.hashCode()
+    }
+
+    /**
+     * Identity-based symbol bindings map
+     */
+    private val variables: Map<String, DerivativeStructureSymbol> = bindings.entries.mapIndexed { index, (key, value) ->
+        key.identity to DerivativeStructureSymbol(numberOfVariables, index, key, value)
+    }.toMap()
+
+    override fun const(value: Double): DerivativeStructure = DerivativeStructure(numberOfVariables, order, value)
+
+    public override fun bindOrNull(symbol: Symbol): DerivativeStructureSymbol? = variables[symbol.identity]
+
+    public fun bind(symbol: Symbol): DerivativeStructureSymbol = variables.getValue(symbol.identity)
+
+    override fun symbol(value: String): DerivativeStructureSymbol = bind(StringSymbol(value))
+
+    public fun DerivativeStructure.derivative(symbols: List<Symbol>): Double {
+        require(symbols.size <= order) { "The order of derivative ${symbols.size} exceeds computed order $order" }
+        val ordersCount = symbols.map { it.identity }.groupBy { it }.mapValues { it.value.size }
+        return getPartialDerivative(*variables.keys.map { ordersCount[it] ?: 0 }.toIntArray())
+    }
+
+    public fun DerivativeStructure.derivative(vararg symbols: Symbol): Double = derivative(symbols.toList())
+
+    public override fun add(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.add(b)
+
+    public override fun multiply(a: DerivativeStructure, k: Number): DerivativeStructure = when (k) {
+        is Double -> a.multiply(k)
+        is Int -> a.multiply(k)
+        else -> a.multiply(k.toDouble())
+    }
+
+    public override fun multiply(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.multiply(b)
+    public override fun divide(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.divide(b)
+    public override fun sin(arg: DerivativeStructure): DerivativeStructure = arg.sin()
+    public override fun cos(arg: DerivativeStructure): DerivativeStructure = arg.cos()
+    public override fun tan(arg: DerivativeStructure): DerivativeStructure = arg.tan()
+    public override fun asin(arg: DerivativeStructure): DerivativeStructure = arg.asin()
+    public override fun acos(arg: DerivativeStructure): DerivativeStructure = arg.acos()
+    public override fun atan(arg: DerivativeStructure): DerivativeStructure = arg.atan()
+    public override fun sinh(arg: DerivativeStructure): DerivativeStructure = arg.sinh()
+    public override fun cosh(arg: DerivativeStructure): DerivativeStructure = arg.cosh()
+    public override fun tanh(arg: DerivativeStructure): DerivativeStructure = arg.tanh()
+    public override fun asinh(arg: DerivativeStructure): DerivativeStructure = arg.asinh()
+    public override fun acosh(arg: DerivativeStructure): DerivativeStructure = arg.acosh()
+    public override fun atanh(arg: DerivativeStructure): DerivativeStructure = arg.atanh()
+
+    public override fun power(arg: DerivativeStructure, pow: Number): DerivativeStructure = when (pow) {
+        is Double -> arg.pow(pow)
+        is Int -> arg.pow(pow)
+        else -> arg.pow(pow.toDouble())
+    }
+
+    public fun power(arg: DerivativeStructure, pow: DerivativeStructure): DerivativeStructure = arg.pow(pow)
+    public override fun exp(arg: DerivativeStructure): DerivativeStructure = arg.exp()
+    public override fun ln(arg: DerivativeStructure): DerivativeStructure = arg.log()
+
+    public override operator fun DerivativeStructure.plus(b: Number): DerivativeStructure = add(b.toDouble())
+    public override operator fun DerivativeStructure.minus(b: Number): DerivativeStructure = subtract(b.toDouble())
+    public override operator fun Number.plus(b: DerivativeStructure): DerivativeStructure = b + this
+    public override operator fun Number.minus(b: DerivativeStructure): DerivativeStructure = b - this
+
+    public companion object :
+        AutoDiffProcessor<Double, DerivativeStructure, DerivativeStructureField, Expression<Double>> {
+        public override fun process(function: DerivativeStructureField.() -> DerivativeStructure): DifferentiableExpression<Double, Expression<Double>> =
+            DerivativeStructureExpression(function)
+    }
+}
+
+
+/**
+ * A constructs that creates a derivative structure with required order on-demand
+ */
+public class DerivativeStructureExpression(
+    public val function: DerivativeStructureField.() -> DerivativeStructure,
+) : DifferentiableExpression<Double, Expression<Double>> {
+    public override operator fun invoke(arguments: Map<Symbol, Double>): Double =
+        DerivativeStructureField(0, arguments).function().value
+
+    /**
+     * Get the derivative expression with given orders
+     */
+    public override fun derivativeOrNull(symbols: List<Symbol>): Expression<Double> = Expression { arguments ->
+        with(DerivativeStructureField(symbols.size, arguments)) { function().derivative(symbols) }
+    }
+}

kmath-commons/src/main/kotlin/kscience/kmath/commons/expressions/DiffExpression.kt

@@ -1,127 +0,0 @@
-package kscience.kmath.commons.expressions
-
-import kscience.kmath.expressions.Expression
-import kscience.kmath.expressions.ExpressionAlgebra
-import kscience.kmath.operations.ExtendedField
-import kscience.kmath.operations.Field
-import kscience.kmath.operations.invoke
-import org.apache.commons.math3.analysis.differentiation.DerivativeStructure
-import kotlin.properties.ReadOnlyProperty
-
-/**
- * A field over commons-math [DerivativeStructure].
- *
- * @property order The derivation order.
- * @property parameters The map of free parameters.
- */
-public class DerivativeStructureField(
-    public val order: Int,
-    public val parameters: Map<String, Double>
-) : ExtendedField<DerivativeStructure> {
-    public override val zero: DerivativeStructure by lazy { DerivativeStructure(parameters.size, order) }
-    public override val one: DerivativeStructure by lazy { DerivativeStructure(parameters.size, order, 1.0) }
-
-    private val variables: Map<String, DerivativeStructure> = parameters.mapValues { (key, value) ->
-        DerivativeStructure(parameters.size, order, parameters.keys.indexOf(key), value)
-    }
-
-    public val variable: ReadOnlyProperty<Any?, DerivativeStructure> = ReadOnlyProperty { _, property ->
-        variables[property.name] ?: error("A variable with name ${property.name} does not exist")
-    }
-
-    public fun variable(name: String, default: DerivativeStructure? = null): DerivativeStructure =
-        variables[name] ?: default ?: error("A variable with name $name does not exist")
-
-    public fun Number.const(): DerivativeStructure = DerivativeStructure(order, parameters.size, toDouble())
-    public fun DerivativeStructure.deriv(parName: String, order: Int = 1): Double = deriv(mapOf(parName to order))
-
-    public fun DerivativeStructure.deriv(orders: Map<String, Int>): Double {
-        return getPartialDerivative(*parameters.keys.map { orders[it] ?: 0 }.toIntArray())
-    }
-
-    public fun DerivativeStructure.deriv(vararg orders: Pair<String, Int>): Double = deriv(mapOf(*orders))
-    public override fun add(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.add(b)
-
-    public override fun multiply(a: DerivativeStructure, k: Number): DerivativeStructure = when (k) {
-        is Double -> a.multiply(k)
-        is Int -> a.multiply(k)
-        else -> a.multiply(k.toDouble())
-    }
-
-    public override fun multiply(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.multiply(b)
-    public override fun divide(a: DerivativeStructure, b: DerivativeStructure): DerivativeStructure = a.divide(b)
-    public override fun sin(arg: DerivativeStructure): DerivativeStructure = arg.sin()
-    public override fun cos(arg: DerivativeStructure): DerivativeStructure = arg.cos()
-    public override fun tan(arg: DerivativeStructure): DerivativeStructure = arg.tan()
-    public override fun asin(arg: DerivativeStructure): DerivativeStructure = arg.asin()
-    public override fun acos(arg: DerivativeStructure): DerivativeStructure = arg.acos()
-    public override fun atan(arg: DerivativeStructure): DerivativeStructure = arg.atan()
-    public override fun sinh(arg: DerivativeStructure): DerivativeStructure = arg.sinh()
-    public override fun cosh(arg: DerivativeStructure): DerivativeStructure = arg.cosh()
-    public override fun tanh(arg: DerivativeStructure): DerivativeStructure = arg.tanh()
-    public override fun asinh(arg: DerivativeStructure): DerivativeStructure = arg.asinh()
-    public override fun acosh(arg: DerivativeStructure): DerivativeStructure = arg.acosh()
-    public override fun atanh(arg: DerivativeStructure): DerivativeStructure = arg.atanh()
-
-    public override fun power(arg: DerivativeStructure, pow: Number): DerivativeStructure = when (pow) {
-        is Double -> arg.pow(pow)
-        is Int -> arg.pow(pow)
-        else -> arg.pow(pow.toDouble())
-    }
-
-    public fun power(arg: DerivativeStructure, pow: DerivativeStructure): DerivativeStructure = arg.pow(pow)
-    public override fun exp(arg: DerivativeStructure): DerivativeStructure = arg.exp()
-    public override fun ln(arg: DerivativeStructure): DerivativeStructure = arg.log()
-    public override operator fun DerivativeStructure.plus(b: Number): DerivativeStructure = add(b.toDouble())
-    public override operator fun DerivativeStructure.minus(b: Number): DerivativeStructure = subtract(b.toDouble())
-    public override operator fun Number.plus(b: DerivativeStructure): DerivativeStructure = b + this
-    public override operator fun Number.minus(b: DerivativeStructure): DerivativeStructure = b - this
-}
-
-/**
- * A constructs that creates a derivative structure with required order on-demand
- */
-public class DiffExpression(public val function: DerivativeStructureField.() -> DerivativeStructure) :
-    Expression<Double> {
-    public override operator fun invoke(arguments: Map<String, Double>): Double = DerivativeStructureField(
-        0,
-        arguments
-    ).function().value
-
-    /**
-     * Get the derivative expression with given orders
-     * TODO make result [DiffExpression]
-     */
-    public fun derivative(orders: Map<String, Int>): Expression<Double> = Expression { arguments ->
-        (DerivativeStructureField(orders.values.maxOrNull() ?: 0, arguments)) { function().deriv(orders) }
-    }
-
-    //TODO add gradient and maybe other vector operators
-}
-
-public fun DiffExpression.derivative(vararg orders: Pair<String, Int>): Expression<Double> = derivative(mapOf(*orders))
-public fun DiffExpression.derivative(name: String): Expression<Double> = derivative(name to 1)
-
-/**
- * A context for [DiffExpression] (not to be confused with [DerivativeStructure])
- */
-public object DiffExpressionAlgebra : ExpressionAlgebra<Double, DiffExpression>, Field<DiffExpression> {
-    public override val zero: DiffExpression = DiffExpression { 0.0.const() }
-    public override val one: DiffExpression = DiffExpression { 1.0.const() }
-
-    public override fun variable(name: String, default: Double?): DiffExpression =
-        DiffExpression { variable(name, default?.const()) }
-
-    public override fun const(value: Double): DiffExpression = DiffExpression { value.const() }
-
-    public override fun add(a: DiffExpression, b: DiffExpression): DiffExpression =
-        DiffExpression { a.function(this) + b.function(this) }
-
-    public override fun multiply(a: DiffExpression, k: Number): DiffExpression = DiffExpression { a.function(this) * k }
-
-    public override fun multiply(a: DiffExpression, b: DiffExpression): DiffExpression =
-        DiffExpression { a.function(this) * b.function(this) }
-
-    public override fun divide(a: DiffExpression, b: DiffExpression): DiffExpression =
-        DiffExpression { a.function(this) / b.function(this) }
-}

kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/CMOptimizationProblem.kt

@@ -0,0 +1,110 @@
+package kscience.kmath.commons.optimization
+
+import kscience.kmath.expressions.*
+import kscience.kmath.stat.OptimizationFeature
+import kscience.kmath.stat.OptimizationProblem
+import kscience.kmath.stat.OptimizationProblemFactory
+import kscience.kmath.stat.OptimizationResult
+import org.apache.commons.math3.optim.*
+import org.apache.commons.math3.optim.nonlinear.scalar.GoalType
+import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer
+import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction
+import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunctionGradient
+import org.apache.commons.math3.optim.nonlinear.scalar.gradient.NonLinearConjugateGradientOptimizer
+import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.AbstractSimplex
+import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.NelderMeadSimplex
+import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.SimplexOptimizer
+import kotlin.reflect.KClass
+
+public operator fun PointValuePair.component1(): DoubleArray = point
+public operator fun PointValuePair.component2(): Double = value
+
+public class CMOptimizationProblem(override val symbols: List<Symbol>, ) :
+    OptimizationProblem<Double>, SymbolIndexer, OptimizationFeature {
+    private val optimizationData: HashMap<KClass<out OptimizationData>, OptimizationData> = HashMap()
+    private var optimizatorBuilder: (() -> MultivariateOptimizer)? = null
+    public var convergenceChecker: ConvergenceChecker<PointValuePair> = SimpleValueChecker(DEFAULT_RELATIVE_TOLERANCE,
+        DEFAULT_ABSOLUTE_TOLERANCE, DEFAULT_MAX_ITER)
+
+    public fun addOptimizationData(data: OptimizationData) {
+        optimizationData[data::class] = data
+    }
+
+    init {
+        addOptimizationData(MaxEval.unlimited())
+    }
+
+    public fun exportOptimizationData(): List<OptimizationData> = optimizationData.values.toList()
+
+    public override fun initialGuess(map: Map<Symbol, Double>): Unit {
+        addOptimizationData(InitialGuess(map.toDoubleArray()))
+    }
+
+    public override fun expression(expression: Expression<Double>): Unit {
+        val objectiveFunction = ObjectiveFunction {
+            val args = it.toMap()
+            expression(args)
+        }
+        addOptimizationData(objectiveFunction)
+    }
+
+    public override fun diffExpression(expression: DifferentiableExpression<Double, Expression<Double>>) {
+        expression(expression)
+        val gradientFunction = ObjectiveFunctionGradient {
+            val args = it.toMap()
+            DoubleArray(symbols.size) { index ->
+                expression.derivative(symbols[index])(args)
+            }
+        }
+        addOptimizationData(gradientFunction)
+        if (optimizatorBuilder == null) {
+            optimizatorBuilder = {
+                NonLinearConjugateGradientOptimizer(
+                    NonLinearConjugateGradientOptimizer.Formula.FLETCHER_REEVES,
+                    convergenceChecker
+                )
+            }
+        }
+    }
+
+    public fun simplex(simplex: AbstractSimplex) {
+        addOptimizationData(simplex)
+        //Set optimization builder to simplex if it is not present
+        if (optimizatorBuilder == null) {
+            optimizatorBuilder = { SimplexOptimizer(convergenceChecker) }
+        }
+    }
+
+    public fun simplexSteps(steps: Map<Symbol, Double>) {
+        simplex(NelderMeadSimplex(steps.toDoubleArray()))
+    }
+
+    public fun goal(goalType: GoalType) {
+        addOptimizationData(goalType)
+    }
+
+    public fun optimizer(block: () -> MultivariateOptimizer) {
+        optimizatorBuilder = block
+    }
+
+    override fun update(result: OptimizationResult<Double>) {
+        initialGuess(result.point)
+    }
+
+    override fun optimize(): OptimizationResult<Double> {
+        val optimizer = optimizatorBuilder?.invoke() ?: error("Optimizer not defined")
+        val (point, value) = optimizer.optimize(*optimizationData.values.toTypedArray())
+        return OptimizationResult(point.toMap(), value, setOf(this))
+    }
+
+    public companion object : OptimizationProblemFactory<Double, CMOptimizationProblem> {
+        public const val DEFAULT_RELATIVE_TOLERANCE: Double = 1e-4
+        public const val DEFAULT_ABSOLUTE_TOLERANCE: Double = 1e-4
+        public const val DEFAULT_MAX_ITER: Int = 1000
+
+        override fun build(symbols: List<Symbol>): CMOptimizationProblem = CMOptimizationProblem(symbols)
+    }
+}
+
+public fun CMOptimizationProblem.initialGuess(vararg pairs: Pair<Symbol, Double>): Unit = initialGuess(pairs.toMap())
+public fun CMOptimizationProblem.simplexSteps(vararg pairs: Pair<Symbol, Double>): Unit = simplexSteps(pairs.toMap())

kmath-commons/src/main/kotlin/kscience/kmath/commons/optimization/cmFit.kt

@@ -0,0 +1,67 @@
+package kscience.kmath.commons.optimization
+
+import kscience.kmath.commons.expressions.DerivativeStructureField
+import kscience.kmath.expressions.DifferentiableExpression
+import kscience.kmath.expressions.Expression
+import kscience.kmath.expressions.Symbol
+import kscience.kmath.stat.Fitting
+import kscience.kmath.stat.OptimizationResult
+import kscience.kmath.stat.optimizeWith
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.asBuffer
+import org.apache.commons.math3.analysis.differentiation.DerivativeStructure
+import org.apache.commons.math3.optim.nonlinear.scalar.GoalType
+
+/**
+ * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation
+ */
+public fun Fitting.chiSquared(
+    x: Buffer<Double>,
+    y: Buffer<Double>,
+    yErr: Buffer<Double>,
+    model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure,
+): DifferentiableExpression<Double, Expression<Double>> = chiSquared(DerivativeStructureField, x, y, yErr, model)
+
+/**
+ * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation
+ */
+public fun Fitting.chiSquared(
+    x: Iterable<Double>,
+    y: Iterable<Double>,
+    yErr: Iterable<Double>,
+    model: DerivativeStructureField.(x: DerivativeStructure) -> DerivativeStructure,
+): DifferentiableExpression<Double, Expression<Double>> = chiSquared(
+    DerivativeStructureField,
+    x.toList().asBuffer(),
+    y.toList().asBuffer(),
+    yErr.toList().asBuffer(),
+    model
+)
+
+/**
+ * Optimize expression without derivatives
+ */
+public fun Expression<Double>.optimize(
+    vararg symbols: Symbol,
+    configuration: CMOptimizationProblem.() -> Unit,
+): OptimizationResult<Double> = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration)
+
+/**
+ * Optimize differentiable expression
+ */
+public fun DifferentiableExpression<Double, Expression<Double>>.optimize(
+    vararg symbols: Symbol,
+    configuration: CMOptimizationProblem.() -> Unit,
+): OptimizationResult<Double> = optimizeWith(CMOptimizationProblem, symbols = symbols, configuration)
+
+public fun DifferentiableExpression<Double, Expression<Double>>.minimize(
+    vararg startPoint: Pair<Symbol, Double>,
+    configuration: CMOptimizationProblem.() -> Unit = {},
+): OptimizationResult<Double> {
+    require(startPoint.isNotEmpty()) { "Must provide a list of symbols for optimization" }
+    val problem = CMOptimizationProblem(startPoint.map { it.first }).apply(configuration)
+    problem.diffExpression(this)
+    problem.initialGuess(startPoint.toMap())
+    problem.goal(GoalType.MINIMIZE)
+    return problem.optimize()
+}

kmath-commons/src/main/kotlin/kscience/kmath/commons/random/CMRandomGeneratorWrapper.kt

@@ -1,9 +1,10 @@
 package kscience.kmath.commons.random
 
-import kscience.kmath.prob.RandomGenerator
+import kscience.kmath.stat.RandomGenerator
 
-public class CMRandomGeneratorWrapper(public val factory: (IntArray) -> RandomGenerator) :
+public class CMRandomGeneratorWrapper(
-    org.apache.commons.math3.random.RandomGenerator {
+    public val factory: (IntArray) -> RandomGenerator,
+) : org.apache.commons.math3.random.RandomGenerator {
     private var generator: RandomGenerator = factory(intArrayOf())
 
     public override fun nextBoolean(): Boolean = generator.nextBoolean()

kmath-commons/src/test/kotlin/kscience/kmath/commons/expressions/AutoDiffTest.kt

@@ -1,40 +0,0 @@
-package kscience.kmath.commons.expressions
-
-import kscience.kmath.expressions.invoke
-import kotlin.contracts.InvocationKind
-import kotlin.contracts.contract
-import kotlin.test.Test
-import kotlin.test.assertEquals
-
-internal inline fun <R> diff(
-    order: Int,
-    vararg parameters: Pair<String, Double>,
-    block: DerivativeStructureField.() -> R
-): R {
-    contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
-    return DerivativeStructureField(order, mapOf(*parameters)).run(block)
-}
-
-internal class AutoDiffTest {
-    @Test
-    fun derivativeStructureFieldTest() {
-        val res = diff(3, "x" to 1.0, "y" to 1.0) {
-            val x by variable
-            val y = variable("y")
-            val z = x * (-sin(x * y) + y)
-            z.deriv("x")
-        }
-    }
-
-    @Test
-    fun autoDifTest() {
-        val f = DiffExpression {
-            val x by variable
-            val y by variable
-            x.pow(2) + 2 * x * y + y.pow(2) + 1
-        }
-
-        assertEquals(10.0, f("x" to 1.0, "y" to 2.0))
-        assertEquals(6.0, f.derivative("x")("x" to 1.0, "y" to 2.0))
-    }
-}

kmath-commons/src/test/kotlin/kscience/kmath/commons/expressions/DerivativeStructureExpressionTest.kt

@@ -0,0 +1,50 @@
+package kscience.kmath.commons.expressions
+
+import kscience.kmath.expressions.*
+import kotlin.contracts.InvocationKind
+import kotlin.contracts.contract
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertFails
+
+internal inline fun diff(
+    order: Int,
+    vararg parameters: Pair<Symbol, Double>,
+    block: DerivativeStructureField.() -> Unit,
+): Unit {
+    contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
+    DerivativeStructureField(order, mapOf(*parameters)).run(block)
+}
+
+internal class AutoDiffTest {
+    private val x by symbol
+    private val y by symbol
+
+    @Test
+    fun derivativeStructureFieldTest() {
+        diff(2, x to 1.0, y to 1.0) {
+            val x = bind(x)//by binding()
+            val y = symbol("y")
+            val z = x * (-sin(x * y) + y) + 2.0
+            println(z.derivative(x))
+            println(z.derivative(y,x))
+            assertEquals(z.derivative(x, y), z.derivative(y, x))
+            //check that improper order cause failure
+            assertFails { z.derivative(x,x,y) }
+        }
+    }
+
+    @Test
+    fun autoDifTest() {
+        val f = DerivativeStructureExpression {
+            val x by binding()
+            val y by binding()
+            x.pow(2) + 2 * x * y + y.pow(2) + 1
+        }
+
+        assertEquals(10.0, f(x to 1.0, y to 2.0))
+        assertEquals(6.0, f.derivative(x)(x to 1.0, y to 2.0))
+        assertEquals(2.0, f.derivative(x, x)(x to 1.234, y to -2.0))
+        assertEquals(2.0, f.derivative(x, y)(x to 1.0, y to 2.0))
+    }
+}

kmath-commons/src/test/kotlin/kscience/kmath/commons/optimization/OptimizeTest.kt

@@ -0,0 +1,68 @@
+package kscience.kmath.commons.optimization
+
+import kscience.kmath.commons.expressions.DerivativeStructureExpression
+import kscience.kmath.expressions.symbol
+import kscience.kmath.stat.Distribution
+import kscience.kmath.stat.Fitting
+import kscience.kmath.stat.RandomGenerator
+import kscience.kmath.stat.normal
+import org.junit.jupiter.api.Test
+import kotlin.math.pow
+
+internal class OptimizeTest {
+    val x by symbol
+    val y by symbol
+
+    val normal = DerivativeStructureExpression {
+        exp(-bind(x).pow(2) / 2) + exp(-bind(y).pow(2) / 2)
+    }
+
+    @Test
+    fun testGradientOptimization() {
+        val result = normal.optimize(x, y) {
+            initialGuess(x to 1.0, y to 1.0)
+            //no need to select optimizer. Gradient optimizer is used by default because gradients are provided by function
+        }
+        println(result.point)
+        println(result.value)
+    }
+
+    @Test
+    fun testSimplexOptimization() {
+        val result = normal.optimize(x, y) {
+            initialGuess(x to 1.0, y to 1.0)
+            simplexSteps(x to 2.0, y to 0.5)
+            //this sets simplex optimizer
+        }
+        println(result.point)
+        println(result.value)
+    }
+
+    @Test
+    fun testCmFit() {
+        val a by symbol
+        val b by symbol
+        val c by symbol
+
+        val sigma = 1.0
+        val generator = Distribution.normal(0.0, sigma)
+        val chain = generator.sample(RandomGenerator.default(112667))
+        val x = (1..100).map(Int::toDouble)
+
+        val y = x.map {
+            it.pow(2) + it + 1 + chain.nextDouble()
+        }
+
+        val yErr = List(x.size) { sigma }
+
+        val chi2 = Fitting.chiSquared(x, y, yErr) { x1 ->
+            val cWithDefault = bindOrNull(c) ?: one
+            bind(a) * x1.pow(2) + bind(b) * x1 + cWithDefault
+        }
+
+        val result = chi2.minimize(a to 1.5, b to 0.9, c to 1.0)
+        println(result)
+        println("Chi2/dof = ${result.value / (x.size - 3)}")
+    }
+
+}

kmath-core/README.md

@@ -7,12 +7,12 @@ The core features of KMath:
  - [buffers](src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : One-dimensional structure
  - [expressions](src/commonMain/kotlin/kscience/kmath/expressions) : Functional Expressions
  - [domains](src/commonMain/kotlin/kscience/kmath/domains) : Domains
- - [autodif](src/commonMain/kotlin/kscience/kmath/misc/AutoDiff.kt) : Automatic differentiation
+ - [autodif](src/commonMain/kotlin/kscience/kmath/expressions/SimpleAutoDiff.kt) : Automatic differentiation
 
 
 > #### Artifact:
 >
-> This module artifact: `kscience.kmath:kmath-core:0.2.0-dev-1`.
+> This module artifact: `kscience.kmath:kmath-core:0.2.0-dev-3`.
 >
 > Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-core/_latestVersion)
 >
@@ -22,25 +22,28 @@ The core features of KMath:
 >
 > ```gradle
 > repositories {
+>     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" }
 >     maven { url 'https://dl.bintray.com/mipt-npm/kscience' }
 >     maven { url 'https://dl.bintray.com/mipt-npm/dev' }
 >     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+
 > }
 > 
 > dependencies {
->     implementation 'kscience.kmath:kmath-core:0.2.0-dev-1'
+>     implementation 'kscience.kmath:kmath-core:0.2.0-dev-3'
 > }
 > ```
 > **Gradle Kotlin DSL:**
 >
 > ```kotlin
 > repositories {
+>     maven("https://dl.bintray.com/kotlin/kotlin-eap")
 >     maven("https://dl.bintray.com/mipt-npm/kscience")
 >     maven("https://dl.bintray.com/mipt-npm/dev")
 >     maven("https://dl.bintray.com/hotkeytlt/maven")
 > }
 > 
 > dependencies {
->     implementation("kscience.kmath:kmath-core:0.2.0-dev-1")
+>     implementation("kscience.kmath:kmath-core:0.2.0-dev-3")
 > }
 > ```

kmath-core/build.gradle.kts

@@ -49,6 +49,6 @@ readme {
     feature(
         id = "autodif",
         description = "Automatic differentiation",
-        ref = "src/commonMain/kotlin/kscience/kmath/misc/AutoDiff.kt"
+        ref = "src/commonMain/kotlin/kscience/kmath/expressions/SimpleAutoDiff.kt"
     )
 }

kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/DifferentiableExpression.kt

@@ -0,0 +1,48 @@
+package kscience.kmath.expressions
+
+/**
+ * Represents expression which structure can be differentiated.
+ *
+ * @param T the type this expression takes as argument and returns.
+ * @param R the type of expression this expression can be differentiated to.
+ */
+public interface DifferentiableExpression<T, out R : Expression<T>> : Expression<T> {
+    /**
+     * Differentiates this expression by ordered collection of [symbols].
+     *
+     * @param symbols the symbols.
+     * @return the derivative or `null`.
+     */
+    public fun derivativeOrNull(symbols: List<Symbol>): R?
+}
+
+public fun <T, R : Expression<T>> DifferentiableExpression<T, R>.derivative(symbols: List<Symbol>): R =
+    derivativeOrNull(symbols) ?: error("Derivative by symbols $symbols not provided")
+
+public fun <T, R : Expression<T>> DifferentiableExpression<T, R>.derivative(vararg symbols: Symbol): R =
+    derivative(symbols.toList())
+
+public fun <T, R : Expression<T>> DifferentiableExpression<T, R>.derivative(name: String): R =
+    derivative(StringSymbol(name))
+
+/**
+ * A [DifferentiableExpression] that defines only first derivatives
+ */
+public abstract class FirstDerivativeExpression<T, R : Expression<T>> : DifferentiableExpression<T,R> {
+    /**
+     * Returns first derivative of this expression by given [symbol].
+     */
+    public abstract fun derivativeOrNull(symbol: Symbol): R?
+
+    public final override fun derivativeOrNull(symbols: List<Symbol>): R? {
+        val dSymbol = symbols.firstOrNull() ?: return null
+        return derivativeOrNull(dSymbol)
+    }
+}
+
+/**
+ * A factory that converts an expression in autodiff variables to a [DifferentiableExpression]
+ */
+public fun interface AutoDiffProcessor<T : Any, I : Any, A : ExpressionAlgebra<T, I>, out R : Expression<T>> {
+    public fun process(function: A.() -> I): DifferentiableExpression<T, R>
+}

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

@@ -1,9 +1,30 @@
 package kscience.kmath.expressions
 
 import kscience.kmath.operations.Algebra
+import kotlin.jvm.JvmName
+import kotlin.properties.ReadOnlyProperty
 
 /**
- * An elementary function that could be invoked on a map of arguments
+ * A marker interface for a symbol. A symbol mus have an identity
+ */
+public interface Symbol {
+    /**
+     * Identity object for the symbol. Two symbols with the same identity are considered to be the same symbol.
+     */
+    public val identity: String
+}
+
+/**
+ * A [Symbol] with a [String] identity
+ */
+public inline class StringSymbol(override val identity: String) : Symbol {
+    override fun toString(): String = identity
+}
+
+/**
+ * An elementary function that could be invoked on a map of arguments.
+ *
+ * @param T the type this expression takes as argument and returns.
  */
 public fun interface Expression<T> {
     /**
@@ -12,30 +33,75 @@ public fun interface Expression<T> {
      * @param arguments the map of arguments.
      * @return the value.
      */
-    public operator fun invoke(arguments: Map<String, T>): T
+    public operator fun invoke(arguments: Map<Symbol, T>): T
-
-    public companion object
 }
 
+/**
+ * Calls this expression without providing any arguments.
+ *
+ * @return a value.
+ */
+public operator fun <T> Expression<T>.invoke(): T = invoke(emptyMap())
+
 /**
  * Calls this expression from arguments.
  *
- * @param pairs the pair of arguments' names to values.
+ * @param pairs the pairs of arguments to values.
- * @return the value.
+ * @return a value.
  */
-public operator fun <T> Expression<T>.invoke(vararg pairs: Pair<String, T>): T = invoke(mapOf(*pairs))
+@JvmName("callBySymbol")
+public operator fun <T> Expression<T>.invoke(vararg pairs: Pair<Symbol, T>): T = invoke(mapOf(*pairs))
+
+/**
+ * Calls this expression from arguments.
+ *
+ * @param pairs the pairs of arguments' names to values.
+ * @return a value.
+ */
+@JvmName("callByString")
+public operator fun <T> Expression<T>.invoke(vararg pairs: Pair<String, T>): T =
+    invoke(mapOf(*pairs).mapKeys { StringSymbol(it.key) })
+
 
 /**
  * A context for expression construction
+ *
+ * @param T type of the constants for the expression
+ * @param E type of the actual expression state
  */
-public interface ExpressionAlgebra<T, E> : Algebra<E> {
+public interface ExpressionAlgebra<in T, E> : Algebra<E> {
     /**
-     * Introduce a variable into expression context
+     * Bind a given [Symbol] to this context variable and produce context-specific object. Return null if symbol could not be bound in current context.
      */
-    public fun variable(name: String, default: T? = null): E
+    public fun bindOrNull(symbol: Symbol): E?
+
+    /**
+     * Bind a string to a context using [StringSymbol]
+     */
+    override fun symbol(value: String): E = bind(StringSymbol(value))
 
     /**
      * A constant expression which does not depend on arguments
      */
     public fun const(value: T): E
 }
+
+/**
+ * Bind a given [Symbol] to this context variable and produce context-specific object.
+ */
+public fun <T, E> ExpressionAlgebra<T, E>.bind(symbol: Symbol): E =
+    bindOrNull(symbol) ?: error("Symbol $symbol could not be bound to $this")
+
+/**
+ * A delegate to create a symbol with a string identity in this scope
+ */
+public val symbol: ReadOnlyProperty<Any?, StringSymbol> = ReadOnlyProperty { _, property ->
+    StringSymbol(property.name)
+}
+
+/**
+ * Bind a symbol by name inside the [ExpressionAlgebra]
+ */
+public fun <T, E> ExpressionAlgebra<T, E>.binding(): ReadOnlyProperty<Any?, E> = ReadOnlyProperty { _, property ->
+    bind(StringSymbol(property.name)) ?: error("A variable with name ${property.name} does not exist")
+}

kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/FunctionalExpressionAlgebra.kt

@@ -2,67 +2,43 @@ package kscience.kmath.expressions
 
 import kscience.kmath.operations.*
 
-internal class FunctionalUnaryOperation<T>(val context: Algebra<T>, val name: String, private val expr: Expression<T>) :
-    Expression<T> {
-    override operator fun invoke(arguments: Map<String, T>): T =
-        context.unaryOperation(name, expr.invoke(arguments))
-}
-
-internal class FunctionalBinaryOperation<T>(
-    val context: Algebra<T>,
-    val name: String,
-    val first: Expression<T>,
-    val second: Expression<T>
-) : Expression<T> {
-    override operator fun invoke(arguments: Map<String, T>): T =
-        context.binaryOperation(name, first.invoke(arguments), second.invoke(arguments))
-}
-
-internal class FunctionalVariableExpression<T>(val name: String, val default: T? = null) : Expression<T> {
-    override operator fun invoke(arguments: Map<String, T>): T =
-        arguments[name] ?: default ?: error("Parameter not found: $name")
-}
-
-internal class FunctionalConstantExpression<T>(val value: T) : Expression<T> {
-    override operator fun invoke(arguments: Map<String, T>): T = value
-}
-
-internal class FunctionalConstProductExpression<T>(
-    val context: Space<T>,
-    private val expr: Expression<T>,
-    val const: Number
-) : Expression<T> {
-    override operator fun invoke(arguments: Map<String, T>): T = context.multiply(expr.invoke(arguments), const)
-}
-
 /**
  * A context class for [Expression] construction.
  *
  * @param algebra The algebra to provide for Expressions built.
  */
-public abstract class FunctionalExpressionAlgebra<T, A : Algebra<T>>(public val algebra: A) :
+public abstract class FunctionalExpressionAlgebra<T, A : Algebra<T>>(
-    ExpressionAlgebra<T, Expression<T>> {
+    public val algebra: A,
+) : ExpressionAlgebra<T, Expression<T>> {
     /**
      * Builds an Expression of constant expression which does not depend on arguments.
      */
-    public override fun const(value: T): Expression<T> = FunctionalConstantExpression(value)
+    public override fun const(value: T): Expression<T> = Expression { value }
 
     /**
      * Builds an Expression to access a variable.
      */
-    public override fun variable(name: String, default: T?): Expression<T> = FunctionalVariableExpression(name, default)
+    public override fun bindOrNull(symbol: Symbol): Expression<T>? = Expression { arguments ->
+        arguments[symbol] ?: error("Argument not found: $symbol")
+    }
 
     /**
      * Builds an Expression of dynamic call of binary operation [operation] on [left] and [right].
      */
-    public override fun binaryOperation(operation: String, left: Expression<T>, right: Expression<T>): Expression<T> =
+    public override fun binaryOperation(
-        FunctionalBinaryOperation(algebra, operation, left, right)
+        operation: String,
+        left: Expression<T>,
+        right: Expression<T>,
+    ): Expression<T> = Expression { arguments ->
+        algebra.binaryOperation(operation, left.invoke(arguments), right.invoke(arguments))
+    }
 
     /**
      * Builds an Expression of dynamic call of unary operation with name [operation] on [arg].
      */
-    public override fun unaryOperation(operation: String, arg: Expression<T>): Expression<T> =
+    public override fun unaryOperation(operation: String, arg: Expression<T>): Expression<T> = Expression { arguments ->
-        FunctionalUnaryOperation(algebra, operation, arg)
+        algebra.unaryOperation(operation, arg.invoke(arguments))
+    }
 }
 
 /**
@@ -81,8 +57,9 @@ public open class FunctionalExpressionSpace<T, A : Space<T>>(algebra: A) :
     /**
      * Builds an Expression of multiplication of expression by number.
      */
-    public override fun multiply(a: Expression<T>, k: Number): Expression<T> =
+    public override fun multiply(a: Expression<T>, k: Number): Expression<T> = Expression { arguments ->
-        FunctionalConstProductExpression(algebra, a, k)
+        algebra.multiply(a.invoke(arguments), k)
+    }
 
     public operator fun Expression<T>.plus(arg: T): Expression<T> = this + const(arg)
     public operator fun Expression<T>.minus(arg: T): Expression<T> = this - const(arg)
@@ -118,8 +95,8 @@ public open class FunctionalExpressionRing<T, A>(algebra: A) : FunctionalExpress
 }
 
 public open class FunctionalExpressionField<T, A>(algebra: A) :
-    FunctionalExpressionRing<T, A>(algebra),
+    FunctionalExpressionRing<T, A>(algebra), Field<Expression<T>>
-    Field<Expression<T>> where A : Field<T>, A : NumericAlgebra<T> {
+        where A : Field<T>, A : NumericAlgebra<T> {
     /**
      * Builds an Expression of division an expression by another one.
      */

kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/SimpleAutoDiff.kt

@@ -0,0 +1,393 @@
+package kscience.kmath.expressions
+
+import kscience.kmath.linear.Point
+import kscience.kmath.operations.*
+import kscience.kmath.structures.asBuffer
+import kotlin.contracts.InvocationKind
+import kotlin.contracts.contract
+
+/*
+ * Implementation of backward-mode automatic differentiation.
+ * Initial gist by Roman Elizarov: https://gist.github.com/elizarov/1ad3a8583e88cb6ea7a0ad09bb591d3d
+ */
+
+
+public open class AutoDiffValue<out T>(public val value: T)
+
+/**
+ * Represents result of [simpleAutoDiff] call.
+ *
+ * @param T the non-nullable type of value.
+ * @param value the value of result.
+ * @property simpleAutoDiff The mapping of differentiated variables to their derivatives.
+ * @property context The field over [T].
+ */
+public class DerivationResult<T : Any>(
+    public val value: T,
+    private val derivativeValues: Map<String, T>,
+    public val context: Field<T>,
+) {
+    /**
+     * Returns derivative of [variable] or returns [Ring.zero] in [context].
+     */
+    public fun derivative(variable: Symbol): T = derivativeValues[variable.identity] ?: context.zero
+
+    /**
+     * Computes the divergence.
+     */
+    public fun div(): T = context { sum(derivativeValues.values) }
+}
+
+/**
+ * Computes the gradient for variables in given order.
+ */
+public fun <T : Any> DerivationResult<T>.grad(vararg variables: Symbol): Point<T> {
+    check(variables.isNotEmpty()) { "Variable order is not provided for gradient construction" }
+    return variables.map(::derivative).asBuffer()
+}
+
+/**
+ * Runs differentiation and establishes [SimpleAutoDiffField] context inside the block of code.
+ *
+ * The partial derivatives are placed in argument `d` variable
+ *
+ * Example:
+ * ```
+ * val x by symbol // define variable(s) and their values
+ * val y = RealField.withAutoDiff() { sqr(x) + 5 * x + 3 } // write formulate in deriv context
+ * assertEquals(17.0, y.x) // the value of result (y)
+ * assertEquals(9.0, x.d)  // dy/dx
+ * ```
+ *
+ * @param body the action in [SimpleAutoDiffField] context returning [AutoDiffVariable] to differentiate with respect to.
+ * @return the result of differentiation.
+ */
+public fun <T : Any, F : Field<T>> F.simpleAutoDiff(
+    bindings: Map<Symbol, T>,
+    body: SimpleAutoDiffField<T, F>.() -> AutoDiffValue<T>,
+): DerivationResult<T> {
+    contract { callsInPlace(body, InvocationKind.EXACTLY_ONCE) }
+
+    return SimpleAutoDiffField(this, bindings).differentiate(body)
+}
+
+public fun <T : Any, F : Field<T>> F.simpleAutoDiff(
+    vararg bindings: Pair<Symbol, T>,
+    body: SimpleAutoDiffField<T, F>.() -> AutoDiffValue<T>,
+): DerivationResult<T> = simpleAutoDiff(bindings.toMap(), body)
+
+/**
+ * Represents field in context of which functions can be derived.
+ */
+public open class SimpleAutoDiffField<T : Any, F : Field<T>>(
+    public val context: F,
+    bindings: Map<Symbol, T>,
+) : Field<AutoDiffValue<T>>, ExpressionAlgebra<T, AutoDiffValue<T>> {
+    public override val zero: AutoDiffValue<T>
+        get() = const(context.zero)
+
+    public override val one: AutoDiffValue<T>
+        get() = const(context.one)
+
+    // this stack contains pairs of blocks and values to apply them to
+    private var stack: Array<Any?> = arrayOfNulls<Any?>(8)
+    private var sp: Int = 0
+    private val derivatives: MutableMap<AutoDiffValue<T>, T> = hashMapOf()
+
+    private val bindings: Map<String, AutoDiffVariableWithDerivative<T>> = bindings.entries.associate {
+        it.key.identity to AutoDiffVariableWithDerivative(it.key.identity, it.value, context.zero)
+    }
+
+    /**
+     * Differentiable variable with value and derivative of differentiation ([simpleAutoDiff]) result
+     * with respect to this variable.
+     *
+     * @param T the non-nullable type of value.
+     * @property value The value of this variable.
+     */
+    private class AutoDiffVariableWithDerivative<T : Any>(
+        override val identity: String,
+        value: T,
+        var d: T,
+    ) : AutoDiffValue<T>(value), Symbol {
+        override fun toString(): String = identity
+        override fun equals(other: Any?): Boolean = this.identity == (other as? Symbol)?.identity
+        override fun hashCode(): Int = identity.hashCode()
+    }
+
+    public override fun bindOrNull(symbol: Symbol): AutoDiffValue<T>? = bindings[symbol.identity]
+
+    private fun getDerivative(variable: AutoDiffValue<T>): T =
+        (variable as? AutoDiffVariableWithDerivative)?.d ?: derivatives[variable] ?: context.zero
+
+    private fun setDerivative(variable: AutoDiffValue<T>, value: T) {
+        if (variable is AutoDiffVariableWithDerivative) variable.d = value else derivatives[variable] = value
+    }
+
+    @Suppress("UNCHECKED_CAST")
+    private fun runBackwardPass() {
+        while (sp > 0) {
+            val value = stack[--sp]
+            val block = stack[--sp] as F.(Any?) -> Unit
+            context.block(value)
+        }
+    }
+
+    override fun const(value: T): AutoDiffValue<T> = AutoDiffValue(value)
+
+    /**
+     * A variable accessing inner state of derivatives.
+     * Use this value in inner builders to avoid creating additional derivative bindings.
+     */
+    public var AutoDiffValue<T>.d: T
+        get() = getDerivative(this)
+        set(value) = setDerivative(this, value)
+
+    public inline fun const(block: F.() -> T): AutoDiffValue<T> = const(context.block())
+
+    /**
+     * Performs update of derivative after the rest of the formula in the back-pass.
+     *
+     * For example, implementation of `sin` function is:
+     *
+     * ```
+     * fun AD.sin(x: Variable): Variable = derive(Variable(sin(x.x)) { z -> // call derive with function result
+     *     x.d += z.d * cos(x.x) // update derivative using chain rule and derivative of the function
+     * }
+     * ```
+     */
+    @Suppress("UNCHECKED_CAST")
+    public fun <R> derive(value: R, block: F.(R) -> Unit): R {
+        // save block to stack for backward pass
+        if (sp >= stack.size) stack = stack.copyOf(stack.size * 2)
+        stack[sp++] = block
+        stack[sp++] = value
+        return value
+    }
+
+
+    internal fun differentiate(function: SimpleAutoDiffField<T, F>.() -> AutoDiffValue<T>): DerivationResult<T> {
+        val result = function()
+        result.d = context.one // computing derivative w.r.t result
+        runBackwardPass()
+        return DerivationResult(result.value, bindings.mapValues { it.value.d }, context)
+    }
+
+    // Overloads for Double constants
+
+    public override operator fun Number.plus(b: AutoDiffValue<T>): AutoDiffValue<T> =
+        derive(const { this@plus.toDouble() * one + b.value }) { z ->
+            b.d += z.d
+        }
+
+    public override operator fun AutoDiffValue<T>.plus(b: Number): AutoDiffValue<T> = b.plus(this)
+
+    public override operator fun Number.minus(b: AutoDiffValue<T>): AutoDiffValue<T> =
+        derive(const { this@minus.toDouble() * one - b.value }) { z -> b.d -= z.d }
+
+    public override operator fun AutoDiffValue<T>.minus(b: Number): AutoDiffValue<T> =
+        derive(const { this@minus.value - one * b.toDouble() }) { z -> this@minus.d += z.d }
+
+
+    // Basic math (+, -, *, /)
+
+    public override fun add(a: AutoDiffValue<T>, b: AutoDiffValue<T>): AutoDiffValue<T> =
+        derive(const { a.value + b.value }) { z ->
+            a.d += z.d
+            b.d += z.d
+        }
+
+    public override fun multiply(a: AutoDiffValue<T>, b: AutoDiffValue<T>): AutoDiffValue<T> =
+        derive(const { a.value * b.value }) { z ->
+            a.d += z.d * b.value
+            b.d += z.d * a.value
+        }
+
+    public override fun divide(a: AutoDiffValue<T>, b: AutoDiffValue<T>): AutoDiffValue<T> =
+        derive(const { a.value / b.value }) { z ->
+            a.d += z.d / b.value
+            b.d -= z.d * a.value / (b.value * b.value)
+        }
+
+    public override fun multiply(a: AutoDiffValue<T>, k: Number): AutoDiffValue<T> =
+        derive(const { k.toDouble() * a.value }) { z ->
+            a.d += z.d * k.toDouble()
+        }
+}
+
+/**
+ * A constructs that creates a derivative structure with required order on-demand
+ */
+public class SimpleAutoDiffExpression<T : Any, F : Field<T>>(
+    public val field: F,
+    public val function: SimpleAutoDiffField<T, F>.() -> AutoDiffValue<T>,
+) : FirstDerivativeExpression<T, Expression<T>>() {
+    public override operator fun invoke(arguments: Map<Symbol, T>): T {
+        //val bindings = arguments.entries.map { it.key.bind(it.value) }
+        return SimpleAutoDiffField(field, arguments).function().value
+    }
+
+    public override fun derivativeOrNull(symbol: Symbol): Expression<T> = Expression { arguments ->
+        //val bindings = arguments.entries.map { it.key.bind(it.value) }
+        val derivationResult = SimpleAutoDiffField(field, arguments).differentiate(function)
+        derivationResult.derivative(symbol)
+    }
+}
+
+/**
+ * Generate [AutoDiffProcessor] for [SimpleAutoDiffExpression]
+ */
+public fun <T : Any, F : Field<T>> simpleAutoDiff(field: F): AutoDiffProcessor<T, AutoDiffValue<T>, SimpleAutoDiffField<T, F>, Expression<T>> =
+    AutoDiffProcessor { function ->
+        SimpleAutoDiffExpression(field, function)
+    }
+
+// Extensions for differentiation of various basic mathematical functions
+
+// x ^ 2
+public fun <T : Any, F : Field<T>> SimpleAutoDiffField<T, F>.sqr(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { x.value * x.value }) { z -> x.d += z.d * 2 * x.value }
+
+// x ^ 1/2
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.sqrt(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { sqrt(x.value) }) { z -> x.d += z.d * 0.5 / z.value }
+
+// x ^ y (const)
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.pow(
+    x: AutoDiffValue<T>,
+    y: Double,
+): AutoDiffValue<T> =
+    derive(const { power(x.value, y) }) { z -> x.d += z.d * y * power(x.value, y - 1) }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.pow(
+    x: AutoDiffValue<T>,
+    y: Int,
+): AutoDiffValue<T> = pow(x, y.toDouble())
+
+// exp(x)
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.exp(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { exp(x.value) }) { z -> x.d += z.d * z.value }
+
+// ln(x)
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.ln(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { ln(x.value) }) { z -> x.d += z.d / x.value }
+
+// x ^ y (any)
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.pow(
+    x: AutoDiffValue<T>,
+    y: AutoDiffValue<T>,
+): AutoDiffValue<T> =
+    exp(y * ln(x))
+
+// sin(x)
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.sin(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { sin(x.value) }) { z -> x.d += z.d * cos(x.value) }
+
+// cos(x)
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.cos(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { cos(x.value) }) { z -> x.d -= z.d * sin(x.value) }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.tan(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { tan(x.value) }) { z ->
+        val c = cos(x.value)
+        x.d += z.d / (c * c)
+    }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.asin(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { asin(x.value) }) { z -> x.d += z.d / sqrt(one - x.value * x.value) }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.acos(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { acos(x.value) }) { z -> x.d -= z.d / sqrt(one - x.value * x.value) }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.atan(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { atan(x.value) }) { z -> x.d += z.d / (one + x.value * x.value) }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.sinh(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { sinh(x.value) }) { z -> x.d += z.d * cosh(x.value) }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.cosh(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { cosh(x.value) }) { z -> x.d += z.d * sinh(x.value) }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.tanh(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { tanh(x.value) }) { z ->
+        val c = cosh(x.value)
+        x.d += z.d / (c * c)
+    }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.asinh(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { asinh(x.value) }) { z -> x.d += z.d / sqrt(one + x.value * x.value) }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.acosh(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { acosh(x.value) }) { z -> x.d += z.d / (sqrt((x.value - one) * (x.value + one))) }
+
+public fun <T : Any, F : ExtendedField<T>> SimpleAutoDiffField<T, F>.atanh(x: AutoDiffValue<T>): AutoDiffValue<T> =
+    derive(const { atanh(x.value) }) { z -> x.d += z.d / (one - x.value * x.value) }
+
+public class SimpleAutoDiffExtendedField<T : Any, F : ExtendedField<T>>(
+    context: F,
+    bindings: Map<Symbol, T>,
+) : ExtendedField<AutoDiffValue<T>>, SimpleAutoDiffField<T, F>(context, bindings) {
+    // x ^ 2
+    public fun sqr(x: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).sqr(x)
+
+    // x ^ 1/2
+    public override fun sqrt(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).sqrt(arg)
+
+    // x ^ y (const)
+    public override fun power(arg: AutoDiffValue<T>, pow: Number): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).pow(arg, pow.toDouble())
+
+    // exp(x)
+    public override fun exp(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).exp(arg)
+
+    // ln(x)
+    public override fun ln(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).ln(arg)
+
+    // x ^ y (any)
+    public fun pow(
+        x: AutoDiffValue<T>,
+        y: AutoDiffValue<T>,
+    ): AutoDiffValue<T> = exp(y * ln(x))
+
+    // sin(x)
+    public override fun sin(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).sin(arg)
+
+    // cos(x)
+    public override fun cos(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).cos(arg)
+
+    public override fun tan(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).tan(arg)
+
+    public override fun asin(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).asin(arg)
+
+    public override fun acos(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).acos(arg)
+
+    public override fun atan(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).atan(arg)
+
+    public override fun sinh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).sinh(arg)
+
+    public override fun cosh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).cosh(arg)
+
+    public override fun tanh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).tanh(arg)
+
+    public override fun asinh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).asinh(arg)
+
+    public override fun acosh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).acosh(arg)
+
+    public override fun atanh(arg: AutoDiffValue<T>): AutoDiffValue<T> =
+        (this as SimpleAutoDiffField<T, F>).atanh(arg)
+}

kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/SymbolIndexer.kt

@@ -0,0 +1,61 @@
+package kscience.kmath.expressions
+
+import kscience.kmath.linear.Point
+import kscience.kmath.structures.BufferFactory
+import kscience.kmath.structures.Structure2D
+
+/**
+ * An environment to easy transform indexed variables to symbols and back.
+ * TODO requires multi-receivers to be beutiful
+ */
+public interface SymbolIndexer {
+    public val symbols: List<Symbol>
+    public fun indexOf(symbol: Symbol): Int = symbols.indexOf(symbol)
+
+    public operator fun <T> List<T>.get(symbol: Symbol): T {
+        require(size == symbols.size) { "The input list size for indexer should be ${symbols.size} but $size found" }
+        return get(this@SymbolIndexer.indexOf(symbol))
+    }
+
+    public operator fun <T> Array<T>.get(symbol: Symbol): T {
+        require(size == symbols.size) { "The input array size for indexer should be ${symbols.size} but $size found" }
+        return get(this@SymbolIndexer.indexOf(symbol))
+    }
+
+    public operator fun DoubleArray.get(symbol: Symbol): Double {
+        require(size == symbols.size) { "The input array size for indexer should be ${symbols.size} but $size found" }
+        return get(this@SymbolIndexer.indexOf(symbol))
+    }
+
+    public operator fun <T> Point<T>.get(symbol: Symbol): T {
+        require(size == symbols.size) { "The input buffer size for indexer should be ${symbols.size} but $size found" }
+        return get(this@SymbolIndexer.indexOf(symbol))
+    }
+
+    public fun DoubleArray.toMap(): Map<Symbol, Double> {
+        require(size == symbols.size) { "The input array size for indexer should be ${symbols.size} but $size found" }
+        return symbols.indices.associate { symbols[it] to get(it) }
+    }
+
+    public operator fun <T> Structure2D<T>.get(rowSymbol: Symbol, columnSymbol: Symbol): T =
+        get(indexOf(rowSymbol), indexOf(columnSymbol))
+
+
+    public fun <T> Map<Symbol, T>.toList(): List<T> = symbols.map { getValue(it) }
+
+    public fun <T> Map<Symbol, T>.toPoint(bufferFactory: BufferFactory<T>): Point<T> =
+        bufferFactory(symbols.size) { getValue(symbols[it]) }
+
+    public fun Map<Symbol, Double>.toDoubleArray(): DoubleArray = DoubleArray(symbols.size) { getValue(symbols[it]) }
+}
+
+public inline class SimpleSymbolIndexer(override val symbols: List<Symbol>) : SymbolIndexer
+
+/**
+ * Execute the block with symbol indexer based on given symbol order
+ */
+public inline fun <R> withSymbols(vararg symbols: Symbol, block: SymbolIndexer.() -> R): R =
+    with(SimpleSymbolIndexer(symbols.toList()), block)
+
+public inline fun <R> withSymbols(symbols: Collection<Symbol>, block: SymbolIndexer.() -> R): R =
+    with(SimpleSymbolIndexer(symbols.toList()), block)

kmath-core/src/commonMain/kotlin/kscience/kmath/expressions/expressionBuilders.kt

@@ -7,6 +7,7 @@ import kscience.kmath.operations.Space
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract
 
+
 /**
  * Creates a functional expression with this [Space].
  */

kmath-core/src/commonMain/kotlin/kscience/kmath/misc/AutoDiff.kt

@@ -1,266 +0,0 @@
-package kscience.kmath.misc
-
-import kscience.kmath.linear.Point
-import kscience.kmath.operations.*
-import kscience.kmath.structures.asBuffer
-import kotlin.contracts.InvocationKind
-import kotlin.contracts.contract
-
-/*
- * Implementation of backward-mode automatic differentiation.
- * Initial gist by Roman Elizarov: https://gist.github.com/elizarov/1ad3a8583e88cb6ea7a0ad09bb591d3d
- */
-
-/**
- * Differentiable variable with value and derivative of differentiation ([deriv]) result
- * with respect to this variable.
- *
- * @param T the non-nullable type of value.
- * @property value The value of this variable.
- */
-public open class Variable<T : Any>(public val value: T)
-
-/**
- * Represents result of [deriv] call.
- *
- * @param T the non-nullable type of value.
- * @param value the value of result.
- * @property deriv The mapping of differentiated variables to their derivatives.
- * @property context The field over [T].
- */
-public class DerivationResult<T : Any>(
-    value: T,
-    public val deriv: Map<Variable<T>, T>,
-    public val context: Field<T>
-) : Variable<T>(value) {
-    /**
-     * Returns derivative of [variable] or returns [Ring.zero] in [context].
-     */
-    public fun deriv(variable: Variable<T>): T = deriv[variable] ?: context.zero
-
-    /**
-     * Computes the divergence.
-     */
-    public fun div(): T = context { sum(deriv.values) }
-
-    /**
-     * Computes the gradient for variables in given order.
-     */
-    public fun grad(vararg variables: Variable<T>): Point<T> {
-        check(variables.isNotEmpty()) { "Variable order is not provided for gradient construction" }
-        return variables.map(::deriv).asBuffer()
-    }
-}
-
-/**
- * Runs differentiation and establishes [AutoDiffField] context inside the block of code.
- *
- * The partial derivatives are placed in argument `d` variable
- *
- * Example:
- * ```
- * val x = Variable(2) // define variable(s) and their values
- * val y = deriv { sqr(x) + 5 * x + 3 } // write formulate in deriv context
- * assertEquals(17.0, y.x) // the value of result (y)
- * assertEquals(9.0, x.d)  // dy/dx
- * ```
- *
- * @param body the action in [AutoDiffField] context returning [Variable] to differentiate with respect to.
- * @return the result of differentiation.
- */
-public inline fun <T : Any, F : Field<T>> F.deriv(body: AutoDiffField<T, F>.() -> Variable<T>): DerivationResult<T> {
-    contract { callsInPlace(body, InvocationKind.EXACTLY_ONCE) }
-
-    return (AutoDiffContext(this)) {
-        val result = body()
-        result.d = context.one // computing derivative w.r.t result
-        runBackwardPass()
-        DerivationResult(result.value, derivatives, this@deriv)
-    }
-}
-
-/**
- * Represents field in context of which functions can be derived.
- */
-public abstract class AutoDiffField<T : Any, F : Field<T>> : Field<Variable<T>> {
-    public abstract val context: F
-
-    /**
-     * A variable accessing inner state of derivatives.
-     * Use this value in inner builders to avoid creating additional derivative bindings.
-     */
-    public abstract var Variable<T>.d: T
-
-    /**
-     * Performs update of derivative after the rest of the formula in the back-pass.
-     *
-     * For example, implementation of `sin` function is:
-     *
-     * ```
-     * fun AD.sin(x: Variable): Variable = derive(Variable(sin(x.x)) { z -> // call derive with function result
-     *     x.d += z.d * cos(x.x) // update derivative using chain rule and derivative of the function
-     * }
-     * ```
-     */
-    public abstract fun <R> derive(value: R, block: F.(R) -> Unit): R
-
-    /**
-     *
-     */
-    public abstract fun variable(value: T): Variable<T>
-
-    public inline fun variable(block: F.() -> T): Variable<T> = variable(context.block())
-
-    // Overloads for Double constants
-
-    override operator fun Number.plus(b: Variable<T>): Variable<T> =
-        derive(variable { this@plus.toDouble() * one + b.value }) { z ->
-            b.d += z.d
-        }
-
-    override operator fun Variable<T>.plus(b: Number): Variable<T> = b.plus(this)
-
-    override operator fun Number.minus(b: Variable<T>): Variable<T> =
-        derive(variable { this@minus.toDouble() * one - b.value }) { z -> b.d -= z.d }
-
-    override operator fun Variable<T>.minus(b: Number): Variable<T> =
-        derive(variable { this@minus.value - one * b.toDouble() }) { z -> this@minus.d += z.d }
-}
-
-/**
- * Automatic Differentiation context class.
- */
-@PublishedApi
-internal class AutoDiffContext<T : Any, F : Field<T>>(override val context: F) : AutoDiffField<T, F>() {
-    // this stack contains pairs of blocks and values to apply them to
-    private var stack: Array<Any?> = arrayOfNulls<Any?>(8)
-    private var sp: Int = 0
-    val derivatives: MutableMap<Variable<T>, T> = hashMapOf()
-    override val zero: Variable<T> get() = Variable(context.zero)
-    override val one: Variable<T> get() = Variable(context.one)
-
-    /**
-     * A variable coupled with its derivative. For internal use only
-     */
-    private class VariableWithDeriv<T : Any>(x: T, var d: T) : Variable<T>(x)
-
-    override fun variable(value: T): Variable<T> =
-        VariableWithDeriv(value, context.zero)
-
-    override var Variable<T>.d: T
-        get() = (this as? VariableWithDeriv)?.d ?: derivatives[this] ?: context.zero
-        set(value) = if (this is VariableWithDeriv) d = value else derivatives[this] = value
-
-    @Suppress("UNCHECKED_CAST")
-    override fun <R> derive(value: R, block: F.(R) -> Unit): R {
-        // save block to stack for backward pass
-        if (sp >= stack.size) stack = stack.copyOf(stack.size * 2)
-        stack[sp++] = block
-        stack[sp++] = value
-        return value
-    }
-
-    @Suppress("UNCHECKED_CAST")
-    fun runBackwardPass() {
-        while (sp > 0) {
-            val value = stack[--sp]
-            val block = stack[--sp] as F.(Any?) -> Unit
-            context.block(value)
-        }
-    }
-
-    // Basic math (+, -, *, /)
-
-    override fun add(a: Variable<T>, b: Variable<T>): Variable<T> = derive(variable { a.value + b.value }) { z ->
-        a.d += z.d
-        b.d += z.d
-    }
-
-    override fun multiply(a: Variable<T>, b: Variable<T>): Variable<T> = derive(variable { a.value * b.value }) { z ->
-        a.d += z.d * b.value
-        b.d += z.d * a.value
-    }
-
-    override fun divide(a: Variable<T>, b: Variable<T>): Variable<T> = derive(variable { a.value / b.value }) { z ->
-        a.d += z.d / b.value
-        b.d -= z.d * a.value / (b.value * b.value)
-    }
-
-    override fun multiply(a: Variable<T>, k: Number): Variable<T> = derive(variable { k.toDouble() * a.value }) { z ->
-        a.d += z.d * k.toDouble()
-    }
-}
-
-// Extensions for differentiation of various basic mathematical functions
-
-// x ^ 2
-public fun <T : Any, F : Field<T>> AutoDiffField<T, F>.sqr(x: Variable<T>): Variable<T> =
-    derive(variable { x.value * x.value }) { z -> x.d += z.d * 2 * x.value }
-
-// x ^ 1/2
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.sqrt(x: Variable<T>): Variable<T> =
-    derive(variable { sqrt(x.value) }) { z -> x.d += z.d * 0.5 / z.value }
-
-// x ^ y (const)
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.pow(x: Variable<T>, y: Double): Variable<T> =
-    derive(variable { power(x.value, y) }) { z -> x.d += z.d * y * power(x.value, y - 1) }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.pow(x: Variable<T>, y: Int): Variable<T> =
-    pow(x, y.toDouble())
-
-// exp(x)
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.exp(x: Variable<T>): Variable<T> =
-    derive(variable { exp(x.value) }) { z -> x.d += z.d * z.value }
-
-// ln(x)
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.ln(x: Variable<T>): Variable<T> =
-    derive(variable { ln(x.value) }) { z -> x.d += z.d / x.value }
-
-// x ^ y (any)
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.pow(x: Variable<T>, y: Variable<T>): Variable<T> =
-    exp(y * ln(x))
-
-// sin(x)
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.sin(x: Variable<T>): Variable<T> =
-    derive(variable { sin(x.value) }) { z -> x.d += z.d * cos(x.value) }
-
-// cos(x)
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.cos(x: Variable<T>): Variable<T> =
-    derive(variable { cos(x.value) }) { z -> x.d -= z.d * sin(x.value) }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.tan(x: Variable<T>): Variable<T> =
-    derive(variable { tan(x.value) }) { z ->
-        val c = cos(x.value)
-        x.d += z.d / (c * c)
-    }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.asin(x: Variable<T>): Variable<T> =
-    derive(variable { asin(x.value) }) { z -> x.d += z.d / sqrt(one - x.value * x.value) }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.acos(x: Variable<T>): Variable<T> =
-    derive(variable { acos(x.value) }) { z -> x.d -= z.d / sqrt(one - x.value * x.value) }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.atan(x: Variable<T>): Variable<T> =
-    derive(variable { atan(x.value) }) { z -> x.d += z.d / (one + x.value * x.value) }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.sinh(x: Variable<T>): Variable<T> =
-    derive(variable { sin(x.value) }) { z -> x.d += z.d * cosh(x.value) }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.cosh(x: Variable<T>): Variable<T> =
-    derive(variable { cos(x.value) }) { z -> x.d += z.d * sinh(x.value) }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.tanh(x: Variable<T>): Variable<T> =
-    derive(variable { tan(x.value) }) { z ->
-        val c = cosh(x.value)
-        x.d += z.d / (c * c)
-    }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.asinh(x: Variable<T>): Variable<T> =
-    derive(variable { asinh(x.value) }) { z -> x.d += z.d / sqrt(one + x.value * x.value) }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.acosh(x: Variable<T>): Variable<T> =
-    derive(variable { acosh(x.value) }) { z -> x.d += z.d / (sqrt((x.value - one) * (x.value + one))) }
-
-public fun <T : Any, F : ExtendedField<T>> AutoDiffField<T, F>.atanh(x: Variable<T>): Variable<T> =
-    derive(variable { atanh(x.value) }) { z -> x.d += z.d / (one - x.value * x.value) }
-

kmath-core/src/commonMain/kotlin/kscience/kmath/operations/AlgebraElements.kt

@@ -74,9 +74,9 @@ public interface SpaceElement<T, I : SpaceElement<T, I, S>, S : Space<T>> : Math
 /**
  * The element of [Ring].
  *
- * @param T the type of space operation results.
+ * @param T the type of ring operation results.
  * @param I self type of the element. Needed for static type checking.
- * @param R the type of space.
+ * @param R the type of ring.
  */
 public interface RingElement<T, I : RingElement<T, I, R>, R : Ring<T>> : SpaceElement<T, I, R> {
     /**
@@ -91,7 +91,7 @@ public interface RingElement<T, I : RingElement<T, I, R>, R : Ring<T>> : SpaceEl
 /**
  * The element of [Field].
  *
- * @param T the type of space operation results.
+ * @param T the type of field operation results.
  * @param I self type of the element. Needed for static type checking.
  * @param F the type of field.
  */

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

@@ -195,6 +195,7 @@ public data class Complex(val re: Double, val im: Double) : FieldElement<Complex
     }
 }
 
+
 /**
  * Creates a complex number with real part equal to this real.
  *

kmath-core/src/commonMain/kotlin/kscience/kmath/structures/NDAlgebra.kt

@@ -73,7 +73,7 @@ public interface NDAlgebra<T, C, N : NDStructure<T>> {
     public fun check(vararg elements: N): Array<out N> = elements
         .map(NDStructure<T>::shape)
         .singleOrNull { !shape.contentEquals(it) }
-        ?.let { throw ShapeMismatchException(shape, it) }
+        ?.let<IntArray, Array<out N>> { throw ShapeMismatchException(shape, it) }
         ?: elements
 
     /**

kmath-core/src/commonTest/kotlin/kscience/kmath/expressions/ExpressionFieldTest.kt

@@ -6,19 +6,21 @@ import kscience.kmath.operations.RealField
 import kscience.kmath.operations.invoke
 import kotlin.test.Test
 import kotlin.test.assertEquals
+import kotlin.test.assertFails
 
 class ExpressionFieldTest {
+    val x by symbol
     @Test
     fun testExpression() {
         val context = FunctionalExpressionField(RealField)
 
         val expression = context {
-            val x = variable("x", 2.0)
+            val x by binding()
             x * x + 2 * x + one
         }
 
-        assertEquals(expression("x" to 1.0), 4.0)
+        assertEquals(expression(x to 1.0), 4.0)
-        assertEquals(expression(), 9.0)
+        assertFails {  expression()}
     }
 
     @Test
@@ -26,33 +28,33 @@ class ExpressionFieldTest {
         val context = FunctionalExpressionField(ComplexField)
 
         val expression = context {
-            val x = variable("x", Complex(2.0, 0.0))
+            val  x = bind(x)
             x * x + 2 * x + one
         }
 
-        assertEquals(expression("x" to Complex(1.0, 0.0)), Complex(4.0, 0.0))
+        assertEquals(expression(x to Complex(1.0, 0.0)), Complex(4.0, 0.0))
-        assertEquals(expression(), Complex(9.0, 0.0))
+        //assertEquals(expression(), Complex(9.0, 0.0))
     }
 
     @Test
     fun separateContext() {
         fun <T> FunctionalExpressionField<T, *>.expression(): Expression<T> {
-            val x = variable("x")
+            val x by binding()
             return x * x + 2 * x + one
         }
 
         val expression = FunctionalExpressionField(RealField).expression()
-        assertEquals(expression("x" to 1.0), 4.0)
+        assertEquals(expression(x to 1.0), 4.0)
     }
 
     @Test
     fun valueExpression() {
         val expressionBuilder: FunctionalExpressionField<Double, *>.() -> Expression<Double> = {
-            val x = variable("x")
+            val x by binding()
             x * x + 2 * x + one
         }
 
         val expression = FunctionalExpressionField(RealField).expressionBuilder()
-        assertEquals(expression("x" to 1.0), 4.0)
+        assertEquals(expression(x to 1.0), 4.0)
     }
 }

kmath-core/src/commonTest/kotlin/kscience/kmath/expressions/SimpleAutoDiffTest.kt

@@ -0,0 +1,285 @@
+package kscience.kmath.expressions
+
+import kscience.kmath.operations.RealField
+import kscience.kmath.structures.asBuffer
+import kotlin.math.E
+import kotlin.math.PI
+import kotlin.math.pow
+import kotlin.math.sqrt
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertTrue
+
+class SimpleAutoDiffTest {
+
+    fun dx(
+        xBinding: Pair<Symbol, Double>,
+        body: SimpleAutoDiffField<Double, RealField>.(x: AutoDiffValue<Double>) -> AutoDiffValue<Double>,
+    ): DerivationResult<Double> = RealField.simpleAutoDiff(xBinding) { body(bind(xBinding.first)) }
+
+    fun dxy(
+        xBinding: Pair<Symbol, Double>,
+        yBinding: Pair<Symbol, Double>,
+        body: SimpleAutoDiffField<Double, RealField>.(x: AutoDiffValue<Double>, y: AutoDiffValue<Double>) -> AutoDiffValue<Double>,
+    ): DerivationResult<Double> = RealField.simpleAutoDiff(xBinding, yBinding) {
+        body(bind(xBinding.first), bind(yBinding.first))
+    }
+
+    fun diff(block: SimpleAutoDiffField<Double, RealField>.() -> AutoDiffValue<Double>): SimpleAutoDiffExpression<Double, RealField> {
+        return SimpleAutoDiffExpression(RealField, block)
+    }
+
+    val x by symbol
+    val y by symbol
+    val z by symbol
+
+    @Test
+    fun testPlusX2() {
+        val y = RealField.simpleAutoDiff(x to 3.0) {
+            // diff w.r.t this x at 3
+            val x = bind(x)
+            x + x
+        }
+        assertEquals(6.0, y.value) //    y  = x + x = 6
+        assertEquals(2.0, y.derivative(x)) // dy/dx = 2
+    }
+
+    @Test
+    fun testPlusX2Expr() {
+        val expr = diff {
+            val x = bind(x)
+            x + x
+        }
+        assertEquals(6.0, expr(x to 3.0)) //    y  = x + x = 6
+        assertEquals(2.0, expr.derivative(x)(x to 3.0)) // dy/dx = 2
+    }
+
+
+    @Test
+    fun testPlus() {
+        // two variables
+        val z = RealField.simpleAutoDiff(x to 2.0, y to 3.0) {
+            val x = bind(x)
+            val y = bind(y)
+            x + y
+        }
+        assertEquals(5.0, z.value) //    z  = x + y = 5
+        assertEquals(1.0, z.derivative(x)) // dz/dx = 1
+        assertEquals(1.0, z.derivative(y)) // dz/dy = 1
+    }
+
+    @Test
+    fun testMinus() {
+        // two variables
+        val z = RealField.simpleAutoDiff(x to 7.0, y to 3.0) {
+            val x = bind(x)
+            val y = bind(y)
+
+            x - y
+        }
+        assertEquals(4.0, z.value)  //    z  = x - y = 4
+        assertEquals(1.0, z.derivative(x))  // dz/dx = 1
+        assertEquals(-1.0, z.derivative(y)) // dz/dy = -1
+    }
+
+    @Test
+    fun testMulX2() {
+        val y = dx(x to 3.0) { x ->
+            // diff w.r.t this x at 3
+            x * x
+        }
+        assertEquals(9.0, y.value) //    y  = x * x = 9
+        assertEquals(6.0, y.derivative(x)) // dy/dx = 2 * x = 7
+    }
+
+    @Test
+    fun testSqr() {
+        val y = dx(x to 3.0) { x -> sqr(x) }
+        assertEquals(9.0, y.value) //    y  = x ^ 2 = 9
+        assertEquals(6.0, y.derivative(x)) // dy/dx = 2 * x = 7
+    }
+
+    @Test
+    fun testSqrSqr() {
+        val y = dx(x to 2.0) { x -> sqr(sqr(x)) }
+        assertEquals(16.0, y.value) //     y = x ^ 4   = 16
+        assertEquals(32.0, y.derivative(x)) // dy/dx = 4 * x^3 = 32
+    }
+
+    @Test
+    fun testX3() {
+        val y = dx(x to 2.0) { x ->
+            // diff w.r.t this x at 2
+            x * x * x
+        }
+        assertEquals(8.0, y.value)  //    y  = x * x * x = 8
+        assertEquals(12.0, y.derivative(x)) // dy/dx = 3 * x * x = 12
+    }
+
+    @Test
+    fun testDiv() {
+        val z = dxy(x to 5.0, y to 2.0) { x, y ->
+            x / y
+        }
+        assertEquals(2.5, z.value)   //     z =  x / y   = 2.5
+        assertEquals(0.5, z.derivative(x))   // dz/dx =  1 / y   = 0.5
+        assertEquals(-1.25, z.derivative(y)) // dz/dy = -x / y^2 = -1.25
+    }
+
+    @Test
+    fun testPow3() {
+        val y = dx(x to 2.0) { x ->
+            // diff w.r.t this x at 2
+            pow(x, 3)
+        }
+        assertEquals(8.0, y.value)  //    y  = x ^ 3     = 8
+        assertEquals(12.0, y.derivative(x)) // dy/dx = 3 * x ^ 2 = 12
+    }
+
+    @Test
+    fun testPowFull() {
+        val z = dxy(x to 2.0, y to 3.0) { x, y ->
+            pow(x, y)
+        }
+        assertApprox(8.0, z.value)           //     z = x ^ y = 8
+        assertApprox(12.0, z.derivative(x))          // dz/dx = y * x ^ (y - 1) = 12
+        assertApprox(8.0 * kotlin.math.ln(2.0), z.derivative(y)) // dz/dy = x ^ y * ln(x)
+    }
+
+    @Test
+    fun testFromPaper() {
+        val y = dx(x to 3.0) { x -> 2 * x + x * x * x }
+        assertEquals(33.0, y.value)  //     y = 2 * x + x * x * x = 33
+        assertEquals(29.0, y.derivative(x))  // dy/dx = 2 + 3 * x * x = 29
+    }
+
+    @Test
+    fun testInnerVariable() {
+        val y = dx(x to 1.0) { x ->
+            const(1.0) * x
+        }
+        assertEquals(1.0, y.value)          //     y = x ^ n = 1
+        assertEquals(1.0, y.derivative(x)) // dy/dx = n * x ^ (n - 1) = n - 1
+    }
+
+    @Test
+    fun testLongChain() {
+        val n = 10_000
+        val y = dx(x to 1.0) { x ->
+            var res = const(1.0)
+            for (i in 1..n) res *= x
+            res
+        }
+        assertEquals(1.0, y.value)          //     y = x ^ n = 1
+        assertEquals(n.toDouble(), y.derivative(x)) // dy/dx = n * x ^ (n - 1) = n - 1
+    }
+
+    @Test
+    fun testExample() {
+        val y = dx(x to 2.0) { x -> sqr(x) + 5 * x + 3 }
+        assertEquals(17.0, y.value) // the value of result (y)
+        assertEquals(9.0, y.derivative(x))  // dy/dx
+    }
+
+    @Test
+    fun testSqrt() {
+        val y = dx(x to 16.0) { x -> sqrt(x) }
+        assertEquals(4.0, y.value)     //     y = x ^ 1/2 = 4
+        assertEquals(1.0 / 8, y.derivative(x)) // dy/dx = 1/2 / x ^ 1/4 = 1/8
+    }
+
+    @Test
+    fun testSin() {
+        val y = dx(x to PI / 6.0) { x -> sin(x) }
+        assertApprox(0.5, y.value) // y = sin(PI/6) = 0.5
+        assertApprox(sqrt(3.0) / 2, y.derivative(x)) // dy/dx = cos(pi/6) = sqrt(3)/2
+    }
+
+    @Test
+    fun testCos() {
+        val y = dx(x to PI / 6) { x -> cos(x) }
+        assertApprox(sqrt(3.0) / 2, y.value) //y = cos(pi/6) = sqrt(3)/2
+        assertApprox(-0.5, y.derivative(x)) // dy/dx = -sin(pi/6) = -0.5
+    }
+
+    @Test
+    fun testTan() {
+        val y = dx(x to PI / 6) { x -> tan(x) }
+        assertApprox(1.0 / sqrt(3.0), y.value) // y = tan(pi/6) = 1/sqrt(3)
+        assertApprox(4.0 / 3.0, y.derivative(x)) // dy/dx = sec(pi/6)^2 = 4/3
+    }
+
+    @Test
+    fun testAsin() {
+        val y = dx(x to PI / 6) { x -> asin(x) }
+        assertApprox(kotlin.math.asin(PI / 6.0), y.value) // y = asin(pi/6)
+        assertApprox(6.0 / sqrt(36 - PI * PI), y.derivative(x)) // dy/dx = 6/sqrt(36-pi^2)
+    }
+
+    @Test
+    fun testAcos() {
+        val y = dx(x to PI / 6) { x -> acos(x) }
+        assertApprox(kotlin.math.acos(PI / 6.0), y.value) // y = acos(pi/6)
+        assertApprox(-6.0 / sqrt(36.0 - PI * PI), y.derivative(x)) // dy/dx = -6/sqrt(36-pi^2)
+    }
+
+    @Test
+    fun testAtan() {
+        val y = dx(x to PI / 6) { x -> atan(x) }
+        assertApprox(kotlin.math.atan(PI / 6.0), y.value) // y = atan(pi/6)
+        assertApprox(36.0 / (36.0 + PI * PI), y.derivative(x)) // dy/dx = 36/(36+pi^2)
+    }
+
+    @Test
+    fun testSinh() {
+        val y = dx(x to 0.0) { x -> sinh(x) }
+        assertApprox(kotlin.math.sinh(0.0), y.value) // y = sinh(0)
+        assertApprox(kotlin.math.cosh(0.0), y.derivative(x)) // dy/dx = cosh(0)
+    }
+
+    @Test
+    fun testCosh() {
+        val y = dx(x to 0.0) { x -> cosh(x) }
+        assertApprox(1.0, y.value) //y = cosh(0)
+        assertApprox(0.0, y.derivative(x)) // dy/dx = sinh(0)
+    }
+
+    @Test
+    fun testTanh() {
+        val y = dx(x to 1.0) { x -> tanh(x) }
+        assertApprox((E * E - 1) / (E * E + 1), y.value) // y = tanh(pi/6)
+        assertApprox(1.0 / kotlin.math.cosh(1.0).pow(2), y.derivative(x)) // dy/dx = sech(pi/6)^2
+    }
+
+    @Test
+    fun testAsinh() {
+        val y = dx(x to PI / 6) { x -> asinh(x) }
+        assertApprox(kotlin.math.asinh(PI / 6.0), y.value) // y = asinh(pi/6)
+        assertApprox(6.0 / sqrt(36 + PI * PI), y.derivative(x)) // dy/dx = 6/sqrt(pi^2+36)
+    }
+
+    @Test
+    fun testAcosh() {
+        val y = dx(x to PI / 6) { x -> acosh(x) }
+        assertApprox(kotlin.math.acosh(PI / 6.0), y.value) // y = acosh(pi/6)
+        assertApprox(-6.0 / sqrt(36.0 - PI * PI), y.derivative(x)) // dy/dx = -6/sqrt(36-pi^2)
+    }
+
+    @Test
+    fun testAtanh() {
+        val y = dx(x to PI / 6) { x -> atanh(x) }
+        assertApprox(kotlin.math.atanh(PI / 6.0), y.value) // y = atanh(pi/6)
+        assertApprox(-36.0 / (PI * PI - 36.0), y.derivative(x)) // dy/dx = -36/(pi^2-36)
+    }
+
+    @Test
+    fun testDivGrad() {
+        val res = dxy(x to 1.0, y to 2.0) { x, y -> x * x + y * y }
+        assertEquals(6.0, res.div())
+        assertTrue(res.grad(x, y).contentEquals(doubleArrayOf(2.0, 4.0).asBuffer()))
+    }
+
+    private fun assertApprox(a: Double, b: Double) {
+        if ((a - b) > 1e-10) assertEquals(a, b)
+    }
+}

kmath-core/src/commonTest/kotlin/kscience/kmath/misc/AutoDiffTest.kt

@@ -1,261 +0,0 @@
-package kscience.kmath.misc
-
-import kscience.kmath.operations.RealField
-import kscience.kmath.structures.asBuffer
-import kotlin.math.PI
-import kotlin.math.pow
-import kotlin.math.sqrt
-import kotlin.test.Test
-import kotlin.test.assertEquals
-import kotlin.test.assertTrue
-
-class AutoDiffTest {
-    inline fun deriv(body: AutoDiffField<Double, RealField>.() -> Variable<Double>): DerivationResult<Double> =
-        RealField.deriv(body)
-
-    @Test
-    fun testPlusX2() {
-        val x = Variable(3.0) // diff w.r.t this x at 3
-        val y = deriv { x + x }
-        assertEquals(6.0, y.value) //    y  = x + x = 6
-        assertEquals(2.0, y.deriv(x)) // dy/dx = 2
-    }
-
-    @Test
-    fun testPlus() {
-        // two variables
-        val x = Variable(2.0)
-        val y = Variable(3.0)
-        val z = deriv { x + y }
-        assertEquals(5.0, z.value) //    z  = x + y = 5
-        assertEquals(1.0, z.deriv(x)) // dz/dx = 1
-        assertEquals(1.0, z.deriv(y)) // dz/dy = 1
-    }
-
-    @Test
-    fun testMinus() {
-        // two variables
-        val x = Variable(7.0)
-        val y = Variable(3.0)
-        val z = deriv { x - y }
-        assertEquals(4.0, z.value)  //    z  = x - y = 4
-        assertEquals(1.0, z.deriv(x))  // dz/dx = 1
-        assertEquals(-1.0, z.deriv(y)) // dz/dy = -1
-    }
-
-    @Test
-    fun testMulX2() {
-        val x = Variable(3.0) // diff w.r.t this x at 3
-        val y = deriv { x * x }
-        assertEquals(9.0, y.value) //    y  = x * x = 9
-        assertEquals(6.0, y.deriv(x)) // dy/dx = 2 * x = 7
-    }
-
-    @Test
-    fun testSqr() {
-        val x = Variable(3.0)
-        val y = deriv { sqr(x) }
-        assertEquals(9.0, y.value) //    y  = x ^ 2 = 9
-        assertEquals(6.0, y.deriv(x)) // dy/dx = 2 * x = 7
-    }
-
-    @Test
-    fun testSqrSqr() {
-        val x = Variable(2.0)
-        val y = deriv { sqr(sqr(x)) }
-        assertEquals(16.0, y.value) //     y = x ^ 4   = 16
-        assertEquals(32.0, y.deriv(x)) // dy/dx = 4 * x^3 = 32
-    }
-
-    @Test
-    fun testX3() {
-        val x = Variable(2.0) // diff w.r.t this x at 2
-        val y = deriv { x * x * x }
-        assertEquals(8.0, y.value)  //    y  = x * x * x = 8
-        assertEquals(12.0, y.deriv(x)) // dy/dx = 3 * x * x = 12
-    }
-
-    @Test
-    fun testDiv() {
-        val x = Variable(5.0)
-        val y = Variable(2.0)
-        val z = deriv { x / y }
-        assertEquals(2.5, z.value)   //     z =  x / y   = 2.5
-        assertEquals(0.5, z.deriv(x))   // dz/dx =  1 / y   = 0.5
-        assertEquals(-1.25, z.deriv(y)) // dz/dy = -x / y^2 = -1.25
-    }
-
-    @Test
-    fun testPow3() {
-        val x = Variable(2.0) // diff w.r.t this x at 2
-        val y = deriv { pow(x, 3) }
-        assertEquals(8.0, y.value)  //    y  = x ^ 3     = 8
-        assertEquals(12.0, y.deriv(x)) // dy/dx = 3 * x ^ 2 = 12
-    }
-
-    @Test
-    fun testPowFull() {
-        val x = Variable(2.0)
-        val y = Variable(3.0)
-        val z = deriv { pow(x, y) }
-        assertApprox(8.0, z.value)           //     z = x ^ y = 8
-        assertApprox(12.0, z.deriv(x))          // dz/dx = y * x ^ (y - 1) = 12
-        assertApprox(8.0 * kotlin.math.ln(2.0), z.deriv(y)) // dz/dy = x ^ y * ln(x)
-    }
-
-    @Test
-    fun testFromPaper() {
-        val x = Variable(3.0)
-        val y = deriv { 2 * x + x * x * x }
-        assertEquals(33.0, y.value)  //     y = 2 * x + x * x * x = 33
-        assertEquals(29.0, y.deriv(x))  // dy/dx = 2 + 3 * x * x = 29
-    }
-
-    @Test
-    fun testInnerVariable() {
-        val x = Variable(1.0)
-        val y = deriv {
-            Variable(1.0) * x
-        }
-        assertEquals(1.0, y.value)          //     y = x ^ n = 1
-        assertEquals(1.0, y.deriv(x)) // dy/dx = n * x ^ (n - 1) = n - 1
-    }
-
-    @Test
-    fun testLongChain() {
-        val n = 10_000
-        val x = Variable(1.0)
-        val y = deriv {
-            var res = Variable(1.0)
-            for (i in 1..n) res *= x
-            res
-        }
-        assertEquals(1.0, y.value)          //     y = x ^ n = 1
-        assertEquals(n.toDouble(), y.deriv(x)) // dy/dx = n * x ^ (n - 1) = n - 1
-    }
-
-    @Test
-    fun testExample() {
-        val x = Variable(2.0)
-        val y = deriv { sqr(x) + 5 * x + 3 }
-        assertEquals(17.0, y.value) // the value of result (y)
-        assertEquals(9.0, y.deriv(x))  // dy/dx
-    }
-
-    @Test
-    fun testSqrt() {
-        val x = Variable(16.0)
-        val y = deriv { sqrt(x) }
-        assertEquals(4.0, y.value)     //     y = x ^ 1/2 = 4
-        assertEquals(1.0 / 8, y.deriv(x)) // dy/dx = 1/2 / x ^ 1/4 = 1/8
-    }
-
-    @Test
-    fun testSin() {
-        val x = Variable(PI / 6.0)
-        val y = deriv { sin(x) }
-        assertApprox(0.5, y.value) // y = sin(PI/6) = 0.5
-        assertApprox(sqrt(3.0) / 2, y.deriv(x)) // dy/dx = cos(pi/6) = sqrt(3)/2
-    }
-
-    @Test
-    fun testCos() {
-        val x = Variable(PI / 6)
-        val y = deriv { cos(x) }
-        assertApprox(sqrt(3.0) / 2, y.value) //y = cos(pi/6) = sqrt(3)/2
-        assertApprox(-0.5, y.deriv(x)) // dy/dx = -sin(pi/6) = -0.5
-    }
-
-    @Test
-    fun testTan() {
-        val x = Variable(PI / 6)
-        val y = deriv { tan(x) }
-        assertApprox(1.0 / sqrt(3.0), y.value) // y = tan(pi/6) = 1/sqrt(3)
-        assertApprox(4.0 / 3.0, y.deriv(x)) // dy/dx = sec(pi/6)^2 = 4/3
-    }
-
-    @Test
-    fun testAsin() {
-        val x = Variable(PI / 6)
-        val y = deriv { asin(x) }
-        assertApprox(kotlin.math.asin(PI / 6.0), y.value) // y = asin(pi/6)
-        assertApprox(6.0 / sqrt(36 - PI * PI), y.deriv(x)) // dy/dx = 6/sqrt(36-pi^2)
-    }
-
-    @Test
-    fun testAcos() {
-        val x = Variable(PI / 6)
-        val y = deriv { acos(x) }
-        assertApprox(kotlin.math.acos(PI / 6.0), y.value) // y = acos(pi/6)
-        assertApprox(-6.0 / sqrt(36.0 - PI * PI), y.deriv(x)) // dy/dx = -6/sqrt(36-pi^2)
-    }
-
-    @Test
-    fun testAtan() {
-        val x = Variable(PI / 6)
-        val y = deriv { atan(x) }
-        assertApprox(kotlin.math.atan(PI / 6.0), y.value) // y = atan(pi/6)
-        assertApprox(36.0 / (36.0 + PI * PI), y.deriv(x)) // dy/dx = 36/(36+pi^2)
-    }
-
-    @Test
-    fun testSinh() {
-        val x = Variable(0.0)
-        val y = deriv { sinh(x) }
-        assertApprox(kotlin.math.sinh(0.0), y.value) // y = sinh(0)
-        assertApprox(kotlin.math.cosh(0.0), y.deriv(x)) // dy/dx = cosh(0)
-    }
-
-    @Test
-    fun testCosh() {
-        val x = Variable(0.0)
-        val y = deriv { cosh(x) }
-        assertApprox(1.0, y.value) //y = cosh(0)
-        assertApprox(0.0, y.deriv(x)) // dy/dx = sinh(0)
-    }
-
-    @Test
-    fun testTanh() {
-        val x = Variable(PI / 6)
-        val y = deriv { tanh(x) }
-        assertApprox(1.0 / sqrt(3.0), y.value) // y = tanh(pi/6)
-        assertApprox(1.0 / kotlin.math.cosh(PI / 6.0).pow(2), y.deriv(x)) // dy/dx = sech(pi/6)^2
-    }
-
-    @Test
-    fun testAsinh() {
-        val x = Variable(PI / 6)
-        val y = deriv { asinh(x) }
-        assertApprox(kotlin.math.asinh(PI / 6.0), y.value) // y = asinh(pi/6)
-        assertApprox(6.0 / sqrt(36 + PI * PI), y.deriv(x)) // dy/dx = 6/sqrt(pi^2+36)
-    }
-
-    @Test
-    fun testAcosh() {
-        val x = Variable(PI / 6)
-        val y = deriv { acosh(x) }
-        assertApprox(kotlin.math.acosh(PI / 6.0), y.value) // y = acosh(pi/6)
-        assertApprox(-6.0 / sqrt(36.0 - PI * PI), y.deriv(x)) // dy/dx = -6/sqrt(36-pi^2)
-    }
-
-    @Test
-    fun testAtanh() {
-        val x = Variable(PI / 6.0)
-        val y = deriv { atanh(x) }
-        assertApprox(kotlin.math.atanh(PI / 6.0), y.value) // y = atanh(pi/6)
-        assertApprox(-36.0 / (PI * PI - 36.0), y.deriv(x)) // dy/dx = -36/(pi^2-36)
-    }
-
-    @Test
-    fun testDivGrad() {
-        val x = Variable(1.0)
-        val y = Variable(2.0)
-        val res = deriv { x * x + y * y }
-        assertEquals(6.0, res.div())
-        assertTrue(res.grad(x, y).contentEquals(doubleArrayOf(2.0, 4.0).asBuffer()))
-    }
-
-    private fun assertApprox(a: Double, b: Double) {
-        if ((a - b) > 1e-10) assertEquals(a, b)
-    }
-}

kmath-functions/src/commonMain/kotlin/kscience/kmath/functions/Polynomial.kt

@@ -8,13 +8,13 @@ import kotlin.contracts.contract
 import kotlin.math.max
 import kotlin.math.pow
 
-// TODO make `inline`, when KT-41771 gets fixed
 /**
  * Polynomial coefficients without fixation on specific context they are applied to
  * @param coefficients constant is the leftmost coefficient
  */
 public inline class Polynomial<T : Any>(public val coefficients: List<T>)
 
+@Suppress("FunctionName")
 public fun <T : Any> Polynomial(vararg coefficients: T): Polynomial<T> = Polynomial(coefficients.toList())
 
 public fun Polynomial<Double>.value(): Double = coefficients.reduceIndexed { index, acc, d -> acc + d.pow(index) }
@@ -33,14 +33,6 @@ public fun <T : Any, C : Ring<T>> Polynomial<T>.value(ring: C, arg: T): T = ring
     res
 }
 
-/**
- * Represent a polynomial as a context-dependent function
- */
-public fun <T : Any, C : Ring<T>> Polynomial<T>.asMathFunction(): MathFunction<T, C, T> =
-    object : MathFunction<T, C, T> {
-        override fun C.invoke(arg: T): T = value(this, arg)
-    }
-
 /**
  * Represent the polynomial as a regular context-less function
  */
@@ -49,7 +41,7 @@ public fun <T : Any, C : Ring<T>> Polynomial<T>.asFunction(ring: C): (T) -> T =
 /**
  * An algebra for polynomials
  */
-public class PolynomialSpace<T : Any, C : Ring<T>>(public val ring: C) : Space<Polynomial<T>> {
+public class PolynomialSpace<T : Any, C : Ring<T>>(private val ring: C) : Space<Polynomial<T>> {
     public override val zero: Polynomial<T> = Polynomial(emptyList())
 
     public override fun add(a: Polynomial<T>, b: Polynomial<T>): Polynomial<T> {

kmath-functions/src/commonMain/kotlin/kscience/kmath/functions/functions.kt

@@ -1,34 +0,0 @@
-package kscience.kmath.functions
-
-import kscience.kmath.operations.Algebra
-import kscience.kmath.operations.RealField
-
-// TODO make fun interface when KT-41770 is fixed
-/**
- * A regular function that could be called only inside specific algebra context
- * @param T source type
- * @param C source algebra constraint
- * @param R result type
- */
-public /*fun*/ interface MathFunction<T, C : Algebra<T>, R> {
-    public operator fun C.invoke(arg: T): R
-}
-
-public fun <R> MathFunction<Double, RealField, R>.invoke(arg: Double): R = RealField.invoke(arg)
-
-/**
- * A suspendable function defined in algebraic context
- */
-// TODO make fun interface, when the new JVM IR is enabled
-public interface SuspendableMathFunction<T, C : Algebra<T>, R> {
-    public suspend operator fun C.invoke(arg: T): R
-}
-
-public suspend fun <R> SuspendableMathFunction<Double, RealField, R>.invoke(arg: Double): R = RealField.invoke(arg)
-
-/**
- * A parametric function with parameter
- */
-public fun interface ParametricFunction<T, P, C : Algebra<T>> {
-    public operator fun C.invoke(arg: T, parameter: P): T
-}

kmath-kotlingrad/build.gradle.kts

@@ -0,0 +1,9 @@
+plugins {
+    id("ru.mipt.npm.jvm")
+}
+
+dependencies {
+    implementation("com.github.breandan:kaliningraph:0.1.2")
+    implementation("com.github.breandan:kotlingrad:0.3.7")
+    api(project(":kmath-ast"))
+}

kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/DifferentiableMstExpression.kt

@@ -0,0 +1,53 @@
+package kscience.kmath.kotlingrad
+
+import edu.umontreal.kotlingrad.experimental.SFun
+import kscience.kmath.ast.MST
+import kscience.kmath.ast.MstAlgebra
+import kscience.kmath.ast.MstExpression
+import kscience.kmath.expressions.DifferentiableExpression
+import kscience.kmath.expressions.Symbol
+import kscience.kmath.operations.NumericAlgebra
+
+/**
+ * Represents wrapper of [MstExpression] implementing [DifferentiableExpression].
+ *
+ * The principle of this API is converting the [mst] to an [SFun], differentiating it with Kotlin∇, then converting
+ * [SFun] back to [MST].
+ *
+ * @param T the type of number.
+ * @param A the [NumericAlgebra] of [T].
+ * @property expr the underlying [MstExpression].
+ */
+public inline class DifferentiableMstExpression<T, A>(public val expr: MstExpression<T, A>) :
+    DifferentiableExpression<T, MstExpression<T, A>> where A : NumericAlgebra<T>, T : Number {
+    public constructor(algebra: A, mst: MST) : this(MstExpression(algebra, mst))
+
+    /**
+     * The [MstExpression.algebra] of [expr].
+     */
+    public val algebra: A
+        get() = expr.algebra
+
+    /**
+     * The [MstExpression.mst] of [expr].
+     */
+    public val mst: MST
+        get() = expr.mst
+
+    public override fun invoke(arguments: Map<Symbol, T>): T = expr(arguments)
+
+    public override fun derivativeOrNull(symbols: List<Symbol>): MstExpression<T, A> = MstExpression(
+        algebra,
+        symbols.map(Symbol::identity)
+            .map(MstAlgebra::symbol)
+            .map { it.toSVar<KMathNumber<T, A>>() }
+            .fold(mst.toSFun(), SFun<KMathNumber<T, A>>::d)
+            .toMst(),
+    )
+}
+
+/**
+ * Wraps this [MstExpression] into [DifferentiableMstExpression].
+ */
+public fun <T : Number, A : NumericAlgebra<T>> MstExpression<T, A>.differentiable(): DifferentiableMstExpression<T, A> =
+    DifferentiableMstExpression(this)

kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/KMathNumber.kt

@@ -0,0 +1,18 @@
+package kscience.kmath.kotlingrad
+
+import edu.umontreal.kotlingrad.experimental.RealNumber
+import edu.umontreal.kotlingrad.experimental.SConst
+import kscience.kmath.operations.NumericAlgebra
+
+/**
+ * Implements [RealNumber] by delegating its functionality to [NumericAlgebra].
+ *
+ * @param T the type of number.
+ * @param A the [NumericAlgebra] of [T].
+ * @property algebra the algebra.
+ * @param value the value of this number.
+ */
+public class KMathNumber<T, A>(public val algebra: A, value: T) :
+    RealNumber<KMathNumber<T, A>, T>(value) where T : Number, A : NumericAlgebra<T> {
+    public override fun wrap(number: Number): SConst<KMathNumber<T, A>> = SConst(algebra.number(number))
+}

kmath-kotlingrad/src/main/kotlin/kscience/kmath/kotlingrad/ScalarsAdapters.kt

@@ -0,0 +1,124 @@
+package kscience.kmath.kotlingrad
+
+import edu.umontreal.kotlingrad.experimental.*
+import kscience.kmath.ast.MST
+import kscience.kmath.ast.MstAlgebra
+import kscience.kmath.ast.MstExtendedField
+import kscience.kmath.ast.MstExtendedField.unaryMinus
+import kscience.kmath.operations.*
+
+/**
+ * Maps [SVar] to [MST.Symbolic] directly.
+ *
+ * @receiver the variable.
+ * @return a node.
+ */
+public fun <X : SFun<X>> SVar<X>.toMst(): MST.Symbolic = MstAlgebra.symbol(name)
+
+/**
+ * Maps [SVar] to [MST.Numeric] directly.
+ *
+ * @receiver the constant.
+ * @return a node.
+ */
+public fun <X : SFun<X>> SConst<X>.toMst(): MST.Numeric = MstAlgebra.number(doubleValue)
+
+/**
+ * Maps [SFun] objects to [MST]. Some unsupported operations like [Derivative] are bound and converted then.
+ * [Power] operation is limited to constant right-hand side arguments.
+ *
+ * Detailed mapping is:
+ *
+ * - [SVar] -> [MstExtendedField.symbol];
+ * - [SConst] -> [MstExtendedField.number];
+ * - [Sum] -> [MstExtendedField.add];
+ * - [Prod] -> [MstExtendedField.multiply];
+ * - [Power] -> [MstExtendedField.power] (limited to constant exponents only);
+ * - [Negative] -> [MstExtendedField.unaryMinus];
+ * - [Log] -> [MstExtendedField.ln] (left) / [MstExtendedField.ln] (right);
+ * - [Sine] -> [MstExtendedField.sin];
+ * - [Cosine] -> [MstExtendedField.cos];
+ * - [Tangent] -> [MstExtendedField.tan];
+ * - [DProd] is vector operation, and it is requested to be evaluated;
+ * - [SComposition] is also requested to be evaluated eagerly;
+ * - [VSumAll] is requested to be evaluated;
+ * - [Derivative] is requested to be evaluated.
+ *
+ * @receiver the scalar function.
+ * @return a node.
+ */
+public fun <X : SFun<X>> SFun<X>.toMst(): MST = MstExtendedField {
+    when (this@toMst) {
+        is SVar -> toMst()
+        is SConst -> toMst()
+        is Sum -> left.toMst() + right.toMst()
+        is Prod -> left.toMst() * right.toMst()
+        is Power -> left.toMst() pow ((right as? SConst<*>)?.doubleValue ?: (right() as SConst<*>).doubleValue)
+        is Negative -> -input.toMst()
+        is Log -> ln(left.toMst()) / ln(right.toMst())
+        is Sine -> sin(input.toMst())
+        is Cosine -> cos(input.toMst())
+        is Tangent -> tan(input.toMst())
+        is DProd -> this@toMst().toMst()
+        is SComposition -> this@toMst().toMst()
+        is VSumAll<X, *> -> this@toMst().toMst()
+        is Derivative -> this@toMst().toMst()
+    }
+}
+
+/**
+ * Maps [MST.Numeric] to [SConst] directly.
+ *
+ * @receiver the node.
+ * @return a new constant.
+ */
+public fun <X : SFun<X>> MST.Numeric.toSConst(): SConst<X> = SConst(value)
+
+/**
+ * Maps [MST.Symbolic] to [SVar] directly.
+ *
+ * @receiver the node.
+ * @param proto the prototype instance.
+ * @return a new variable.
+ */
+internal fun <X : SFun<X>> MST.Symbolic.toSVar(): SVar<X> = SVar(value)
+
+/**
+ * Maps [MST] objects to [SFun]. Unsupported operations throw [IllegalStateException].
+ *
+ * Detailed mapping is:
+ *
+ * - [MST.Numeric] -> [SConst];
+ * - [MST.Symbolic] -> [SVar];
+ * - [MST.Unary] -> [Negative], [Sine], [Cosine], [Tangent], [Power], [Log];
+ * - [MST.Binary] -> [Sum], [Prod], [Power].
+ *
+ * @receiver the node.
+ * @param proto the prototype instance.
+ * @return a scalar function.
+ */
+public fun <X : SFun<X>> MST.toSFun(): SFun<X> = when (this) {
+    is MST.Numeric -> toSConst()
+    is MST.Symbolic -> toSVar()
+
+    is MST.Unary -> when (operation) {
+        SpaceOperations.PLUS_OPERATION -> +value.toSFun<X>()
+        SpaceOperations.MINUS_OPERATION -> -value.toSFun<X>()
+        TrigonometricOperations.SIN_OPERATION -> sin(value.toSFun())
+        TrigonometricOperations.COS_OPERATION -> cos(value.toSFun())
+        TrigonometricOperations.TAN_OPERATION -> tan(value.toSFun())
+        PowerOperations.SQRT_OPERATION -> sqrt(value.toSFun())
+        ExponentialOperations.EXP_OPERATION -> exp(value.toSFun())
+        ExponentialOperations.LN_OPERATION -> value.toSFun<X>().ln()
+        else -> error("Unary operation $operation not defined in $this")
+    }
+
+    is MST.Binary -> when (operation) {
+        SpaceOperations.PLUS_OPERATION -> left.toSFun<X>() + right.toSFun()
+        SpaceOperations.MINUS_OPERATION -> left.toSFun<X>() - right.toSFun()
+        RingOperations.TIMES_OPERATION -> left.toSFun<X>() * right.toSFun()
+        FieldOperations.DIV_OPERATION -> left.toSFun<X>() / right.toSFun()
+        PowerOperations.POW_OPERATION -> left.toSFun<X>() pow (right as MST.Numeric).toSConst()
+        else -> error("Binary operation $operation not defined in $this")
+    }
+}

kmath-kotlingrad/src/test/kotlin/kscience/kmath/kotlingrad/AdaptingTests.kt

@@ -0,0 +1,64 @@
+package kscience.kmath.kotlingrad
+
+import edu.umontreal.kotlingrad.experimental.*
+import kscience.kmath.asm.compile
+import kscience.kmath.ast.MstAlgebra
+import kscience.kmath.ast.MstExpression
+import kscience.kmath.ast.parseMath
+import kscience.kmath.expressions.invoke
+import kscience.kmath.operations.RealField
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertTrue
+import kotlin.test.fail
+
+internal class AdaptingTests {
+    @Test
+    fun symbol() {
+        val c1 = MstAlgebra.symbol("x")
+        assertTrue(c1.toSVar<KMathNumber<Double, RealField>>().name == "x")
+        val c2 = "kitten".parseMath().toSFun<KMathNumber<Double, RealField>>()
+        if (c2 is SVar) assertTrue(c2.name == "kitten") else fail()
+    }
+
+    @Test
+    fun number() {
+        val c1 = MstAlgebra.number(12354324)
+        assertTrue(c1.toSConst<DReal>().doubleValue == 12354324.0)
+        val c2 = "0.234".parseMath().toSFun<KMathNumber<Double, RealField>>()
+        if (c2 is SConst) assertTrue(c2.doubleValue == 0.234) else fail()
+        val c3 = "1e-3".parseMath().toSFun<KMathNumber<Double, RealField>>()
+        if (c3 is SConst) assertEquals(0.001, c3.value) else fail()
+    }
+
+    @Test
+    fun simpleFunctionShape() {
+        val linear = "2*x+16".parseMath().toSFun<KMathNumber<Double, RealField>>()
+        if (linear !is Sum) fail()
+        if (linear.left !is Prod) fail()
+        if (linear.right !is SConst) fail()
+    }
+
+    @Test
+    fun simpleFunctionDerivative() {
+        val x = MstAlgebra.symbol("x").toSVar<KMathNumber<Double, RealField>>()
+        val quadratic = "x^2-4*x-44".parseMath().toSFun<KMathNumber<Double, RealField>>()
+        val actualDerivative = MstExpression(RealField, quadratic.d(x).toMst()).compile()
+        val expectedDerivative = MstExpression(RealField, "2*x-4".parseMath()).compile()
+        assertEquals(actualDerivative("x" to 123.0), expectedDerivative("x" to 123.0))
+    }
+
+    @Test
+    fun moreComplexDerivative() {
+        val x = MstAlgebra.symbol("x").toSVar<KMathNumber<Double, RealField>>()
+        val composition = "-sqrt(sin(x^2)-cos(x)^2-16*x)".parseMath().toSFun<KMathNumber<Double, RealField>>()
+        val actualDerivative = MstExpression(RealField, composition.d(x).toMst()).compile()
+
+        val expectedDerivative = MstExpression(
+            RealField,
+            "-(2*x*cos(x^2)+2*sin(x)*cos(x)-16)/(2*sqrt(sin(x^2)-16*x-cos(x)^2))".parseMath()
+        ).compile()
+
+        assertEquals(actualDerivative("x" to 0.1), expectedDerivative("x" to 0.1))
+    }
+}

kmath-nd4j/README.md

@@ -0,0 +1,82 @@
+# ND4J NDStructure implementation (`kmath-nd4j`)
+
+This subproject implements the following features:
+
+ - [nd4jarraystrucure](src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt) : NDStructure wrapper for INDArray
+ - [nd4jarrayrings](src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt) : Rings over Nd4jArrayStructure of Int and Long
+ - [nd4jarrayfields](src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt) : Fields over Nd4jArrayStructure of Float and Double
+
+
+> #### Artifact:
+>
+> This module artifact: `kscience.kmath:kmath-nd4j:0.2.0-dev-3`.
+>
+> Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-nd4j/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-nd4j/_latestVersion)
+>
+> Bintray development version:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-nd4j/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-nd4j/_latestVersion)
+>
+> **Gradle:**
+>
+> ```gradle
+> repositories {
+>     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" }
+>     maven { url 'https://dl.bintray.com/mipt-npm/kscience' }
+>     maven { url 'https://dl.bintray.com/mipt-npm/dev' }
+>     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+
+> }
+> 
+> dependencies {
+>     implementation 'kscience.kmath:kmath-nd4j:0.2.0-dev-3'
+> }
+> ```
+> **Gradle Kotlin DSL:**
+>
+> ```kotlin
+> repositories {
+>     maven("https://dl.bintray.com/kotlin/kotlin-eap")
+>     maven("https://dl.bintray.com/mipt-npm/kscience")
+>     maven("https://dl.bintray.com/mipt-npm/dev")
+>     maven("https://dl.bintray.com/hotkeytlt/maven")
+> }
+> 
+> dependencies {
+>     implementation("kscience.kmath:kmath-nd4j:0.2.0-dev-3")
+> }
+> ```
+
+## Examples
+
+NDStructure wrapper for INDArray:
+
+```kotlin
+import org.nd4j.linalg.factory.*
+import scientifik.kmath.nd4j.*
+import scientifik.kmath.structures.*
+
+val array = Nd4j.ones(2, 2).asRealStructure()
+println(array[0, 0]) // 1.0
+array[intArrayOf(0, 0)] = 24.0
+println(array[0, 0]) // 24.0
+```
+
+Fast element-wise and in-place arithmetics for INDArray:
+
+```kotlin
+import org.nd4j.linalg.factory.*
+import scientifik.kmath.nd4j.*
+import scientifik.kmath.operations.*
+
+val field = RealNd4jArrayField(intArrayOf(2, 2))
+val array = Nd4j.rand(2, 2).asRealStructure()
+
+val res = field {
+    (25.0 / array + 20) * 4
+}
+
+println(res.ndArray)
+// [[  250.6449,  428.5840], 
+//  [  269.7913,  202.2077]]
+```
+
+Contributed by [Iaroslav Postovalov](https://github.com/CommanderTvis).

kmath-nd4j/build.gradle.kts

@@ -0,0 +1,37 @@
+import ru.mipt.npm.gradle.Maturity
+
+plugins {
+    id("ru.mipt.npm.jvm")
+}
+
+dependencies {
+    api(project(":kmath-core"))
+    api("org.nd4j:nd4j-api:1.0.0-beta7")
+    testImplementation("org.deeplearning4j:deeplearning4j-core:1.0.0-beta7")
+    testImplementation("org.nd4j:nd4j-native-platform:1.0.0-beta7")
+    testImplementation("org.slf4j:slf4j-simple:1.7.30")
+}
+
+readme {
+    description = "ND4J NDStructure implementation and according NDAlgebra classes"
+    maturity = Maturity.EXPERIMENTAL
+    propertyByTemplate("artifact", rootProject.file("docs/templates/ARTIFACT-TEMPLATE.md"))
+
+    feature(
+        id = "nd4jarraystructure",
+        description = "NDStructure wrapper for INDArray",
+        ref = "src/commonMain/kotlin/kscience/kmath/operations/Algebra.kt"
+    )
+
+    feature(
+        id = "nd4jarrayrings",
+        description = "Rings over Nd4jArrayStructure of Int and Long",
+        ref = "src/commonMain/kotlin/kscience/kmath/structures/NDStructure.kt"
+    )
+
+    feature(
+        id = "nd4jarrayfields",
+        description = "Fields over Nd4jArrayStructure of Float and Double",
+        ref = "src/commonMain/kotlin/kscience/kmath/structures/Buffers.kt"
+    )
+}

kmath-nd4j/docs/README-TEMPLATE.md

@@ -0,0 +1,43 @@
+# ND4J NDStructure implementation (`kmath-nd4j`)
+
+This subproject implements the following features:
+
+${features}
+
+${artifact}
+
+## Examples
+
+NDStructure wrapper for INDArray:
+
+```kotlin
+import org.nd4j.linalg.factory.*
+import scientifik.kmath.nd4j.*
+import scientifik.kmath.structures.*
+
+val array = Nd4j.ones(2, 2).asRealStructure()
+println(array[0, 0]) // 1.0
+array[intArrayOf(0, 0)] = 24.0
+println(array[0, 0]) // 24.0
+```
+
+Fast element-wise and in-place arithmetics for INDArray:
+
+```kotlin
+import org.nd4j.linalg.factory.*
+import scientifik.kmath.nd4j.*
+import scientifik.kmath.operations.*
+
+val field = RealNd4jArrayField(intArrayOf(2, 2))
+val array = Nd4j.rand(2, 2).asRealStructure()
+
+val res = field {
+    (25.0 / array + 20) * 4
+}
+
+println(res.ndArray)
+// [[  250.6449,  428.5840], 
+//  [  269.7913,  202.2077]]
+```
+
+Contributed by [Iaroslav Postovalov](https://github.com/CommanderTvis).

kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayAlgebra.kt

@@ -0,0 +1,349 @@
+package kscience.kmath.nd4j
+
+import kscience.kmath.operations.*
+import kscience.kmath.structures.NDAlgebra
+import kscience.kmath.structures.NDField
+import kscience.kmath.structures.NDRing
+import kscience.kmath.structures.NDSpace
+import org.nd4j.linalg.api.ndarray.INDArray
+import org.nd4j.linalg.factory.Nd4j
+
+/**
+ * Represents [NDAlgebra] over [Nd4jArrayAlgebra].
+ *
+ * @param T the type of ND-structure element.
+ * @param C the type of the element context.
+ */
+public interface Nd4jArrayAlgebra<T, C> : NDAlgebra<T, C, Nd4jArrayStructure<T>> {
+    /**
+     * Wraps [INDArray] to [N].
+     */
+    public fun INDArray.wrap(): Nd4jArrayStructure<T>
+
+    public override fun produce(initializer: C.(IntArray) -> T): Nd4jArrayStructure<T> {
+        val struct = Nd4j.create(*shape)!!.wrap()
+        struct.indicesIterator().forEach { struct[it] = elementContext.initializer(it) }
+        return struct
+    }
+
+    public override fun map(arg: Nd4jArrayStructure<T>, transform: C.(T) -> T): Nd4jArrayStructure<T> {
+        check(arg)
+        val newStruct = arg.ndArray.dup().wrap()
+        newStruct.elements().forEach { (idx, value) -> newStruct[idx] = elementContext.transform(value) }
+        return newStruct
+    }
+
+    public override fun mapIndexed(
+        arg: Nd4jArrayStructure<T>,
+        transform: C.(index: IntArray, T) -> T
+    ): Nd4jArrayStructure<T> {
+        check(arg)
+        val new = Nd4j.create(*shape).wrap()
+        new.indicesIterator().forEach { idx -> new[idx] = elementContext.transform(idx, arg[idx]) }
+        return new
+    }
+
+    public override fun combine(
+        a: Nd4jArrayStructure<T>,
+        b: Nd4jArrayStructure<T>,
+        transform: C.(T, T) -> T
+    ): Nd4jArrayStructure<T> {
+        check(a, b)
+        val new = Nd4j.create(*shape).wrap()
+        new.indicesIterator().forEach { idx -> new[idx] = elementContext.transform(a[idx], b[idx]) }
+        return new
+    }
+}
+
+/**
+ * Represents [NDSpace] over [Nd4jArrayStructure].
+ *
+ * @param T the type of the element contained in ND structure.
+ * @param S the type of space of structure elements.
+ */
+public interface Nd4jArraySpace<T, S> : NDSpace<T, S, Nd4jArrayStructure<T>>,
+    Nd4jArrayAlgebra<T, S> where S : Space<T> {
+    public override val zero: Nd4jArrayStructure<T>
+        get() = Nd4j.zeros(*shape).wrap()
+
+    public override fun add(a: Nd4jArrayStructure<T>, b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
+        check(a, b)
+        return a.ndArray.add(b.ndArray).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<T>.minus(b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
+        check(this, b)
+        return ndArray.sub(b.ndArray).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<T>.unaryMinus(): Nd4jArrayStructure<T> {
+        check(this)
+        return ndArray.neg().wrap()
+    }
+
+    public override fun multiply(a: Nd4jArrayStructure<T>, k: Number): Nd4jArrayStructure<T> {
+        check(a)
+        return a.ndArray.mul(k).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<T>.div(k: Number): Nd4jArrayStructure<T> {
+        check(this)
+        return ndArray.div(k).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<T>.times(k: Number): Nd4jArrayStructure<T> {
+        check(this)
+        return ndArray.mul(k).wrap()
+    }
+}
+
+/**
+ * Represents [NDRing] over [Nd4jArrayStructure].
+ *
+ * @param T the type of the element contained in ND structure.
+ * @param R the type of ring of structure elements.
+ */
+public interface Nd4jArrayRing<T, R> : NDRing<T, R, Nd4jArrayStructure<T>>, Nd4jArraySpace<T, R> where R : Ring<T> {
+    public override val one: Nd4jArrayStructure<T>
+        get() = Nd4j.ones(*shape).wrap()
+
+    public override fun multiply(a: Nd4jArrayStructure<T>, b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
+        check(a, b)
+        return a.ndArray.mul(b.ndArray).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<T>.minus(b: Number): Nd4jArrayStructure<T> {
+        check(this)
+        return ndArray.sub(b).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<T>.plus(b: Number): Nd4jArrayStructure<T> {
+        check(this)
+        return ndArray.add(b).wrap()
+    }
+
+    public override operator fun Number.minus(b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
+        check(b)
+        return b.ndArray.rsub(this).wrap()
+    }
+
+    public companion object {
+        private val intNd4jArrayRingCache: ThreadLocal<MutableMap<IntArray, IntNd4jArrayRing>> =
+            ThreadLocal.withInitial { hashMapOf() }
+
+        private val longNd4jArrayRingCache: ThreadLocal<MutableMap<IntArray, LongNd4jArrayRing>> =
+            ThreadLocal.withInitial { hashMapOf() }
+
+        /**
+         * Creates an [NDRing] for [Int] values or pull it from cache if it was created previously.
+         */
+        public fun int(vararg shape: Int): Nd4jArrayRing<Int, IntRing> =
+            intNd4jArrayRingCache.get().getOrPut(shape) { IntNd4jArrayRing(shape) }
+
+        /**
+         * Creates an [NDRing] for [Long] values or pull it from cache if it was created previously.
+         */
+        public fun long(vararg shape: Int): Nd4jArrayRing<Long, LongRing> =
+            longNd4jArrayRingCache.get().getOrPut(shape) { LongNd4jArrayRing(shape) }
+
+        /**
+         * Creates a most suitable implementation of [NDRing] using reified class.
+         */
+        @Suppress("UNCHECKED_CAST")
+        public inline fun <reified T : Any> auto(vararg shape: Int): Nd4jArrayRing<T, out Ring<T>> = when {
+            T::class == Int::class -> int(*shape) as Nd4jArrayRing<T, out Ring<T>>
+            T::class == Long::class -> long(*shape) as Nd4jArrayRing<T, out Ring<T>>
+            else -> throw UnsupportedOperationException("This factory method only supports Int and Long types.")
+        }
+    }
+}
+
+/**
+ * Represents [NDField] over [Nd4jArrayStructure].
+ *
+ * @param T the type of the element contained in ND structure.
+ * @param N the type of ND structure.
+ * @param F the type field of structure elements.
+ */
+public interface Nd4jArrayField<T, F> : NDField<T, F, Nd4jArrayStructure<T>>, Nd4jArrayRing<T, F> where F : Field<T> {
+    public override fun divide(a: Nd4jArrayStructure<T>, b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
+        check(a, b)
+        return a.ndArray.div(b.ndArray).wrap()
+    }
+
+    public override operator fun Number.div(b: Nd4jArrayStructure<T>): Nd4jArrayStructure<T> {
+        check(b)
+        return b.ndArray.rdiv(this).wrap()
+    }
+
+
+    public companion object {
+        private val floatNd4jArrayFieldCache: ThreadLocal<MutableMap<IntArray, FloatNd4jArrayField>> =
+            ThreadLocal.withInitial { hashMapOf() }
+
+        private val realNd4jArrayFieldCache: ThreadLocal<MutableMap<IntArray, RealNd4jArrayField>> =
+            ThreadLocal.withInitial { hashMapOf() }
+
+        /**
+         * Creates an [NDField] for [Float] values or pull it from cache if it was created previously.
+         */
+        public fun float(vararg shape: Int): Nd4jArrayRing<Float, FloatField> =
+            floatNd4jArrayFieldCache.get().getOrPut(shape) { FloatNd4jArrayField(shape) }
+
+        /**
+         * Creates an [NDField] for [Double] values or pull it from cache if it was created previously.
+         */
+        public fun real(vararg shape: Int): Nd4jArrayRing<Double, RealField> =
+            realNd4jArrayFieldCache.get().getOrPut(shape) { RealNd4jArrayField(shape) }
+
+        /**
+         * Creates a most suitable implementation of [NDRing] using reified class.
+         */
+        @Suppress("UNCHECKED_CAST")
+        public inline fun <reified T : Any> auto(vararg shape: Int): Nd4jArrayField<T, out Field<T>> = when {
+            T::class == Float::class -> float(*shape) as Nd4jArrayField<T, out Field<T>>
+            T::class == Double::class -> real(*shape) as Nd4jArrayField<T, out Field<T>>
+            else -> throw UnsupportedOperationException("This factory method only supports Float and Double types.")
+        }
+    }
+}
+
+/**
+ * Represents [NDField] over [Nd4jArrayRealStructure].
+ */
+public class RealNd4jArrayField(public override val shape: IntArray) : Nd4jArrayField<Double, RealField> {
+    public override val elementContext: RealField
+        get() = RealField
+
+    public override fun INDArray.wrap(): Nd4jArrayStructure<Double> = check(asRealStructure())
+
+    public override operator fun Nd4jArrayStructure<Double>.div(arg: Double): Nd4jArrayStructure<Double> {
+        check(this)
+        return ndArray.div(arg).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<Double>.plus(arg: Double): Nd4jArrayStructure<Double> {
+        check(this)
+        return ndArray.add(arg).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<Double>.minus(arg: Double): Nd4jArrayStructure<Double> {
+        check(this)
+        return ndArray.sub(arg).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<Double>.times(arg: Double): Nd4jArrayStructure<Double> {
+        check(this)
+        return ndArray.mul(arg).wrap()
+    }
+
+    public override operator fun Double.div(arg: Nd4jArrayStructure<Double>): Nd4jArrayStructure<Double> {
+        check(arg)
+        return arg.ndArray.rdiv(this).wrap()
+    }
+
+    public override operator fun Double.minus(arg: Nd4jArrayStructure<Double>): Nd4jArrayStructure<Double> {
+        check(arg)
+        return arg.ndArray.rsub(this).wrap()
+    }
+}
+
+/**
+ * Represents [NDField] over [Nd4jArrayStructure] of [Float].
+ */
+public class FloatNd4jArrayField(public override val shape: IntArray) : Nd4jArrayField<Float, FloatField> {
+    public override val elementContext: FloatField
+        get() = FloatField
+
+    public override fun INDArray.wrap(): Nd4jArrayStructure<Float> = check(asFloatStructure())
+
+    public override operator fun Nd4jArrayStructure<Float>.div(arg: Float): Nd4jArrayStructure<Float> {
+        check(this)
+        return ndArray.div(arg).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<Float>.plus(arg: Float): Nd4jArrayStructure<Float> {
+        check(this)
+        return ndArray.add(arg).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<Float>.minus(arg: Float): Nd4jArrayStructure<Float> {
+        check(this)
+        return ndArray.sub(arg).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<Float>.times(arg: Float): Nd4jArrayStructure<Float> {
+        check(this)
+        return ndArray.mul(arg).wrap()
+    }
+
+    public override operator fun Float.div(arg: Nd4jArrayStructure<Float>): Nd4jArrayStructure<Float> {
+        check(arg)
+        return arg.ndArray.rdiv(this).wrap()
+    }
+
+    public override operator fun Float.minus(arg: Nd4jArrayStructure<Float>): Nd4jArrayStructure<Float> {
+        check(arg)
+        return arg.ndArray.rsub(this).wrap()
+    }
+}
+
+/**
+ * Represents [NDRing] over [Nd4jArrayIntStructure].
+ */
+public class IntNd4jArrayRing(public override val shape: IntArray) : Nd4jArrayRing<Int, IntRing> {
+    public override val elementContext: IntRing
+        get() = IntRing
+
+    public override fun INDArray.wrap(): Nd4jArrayStructure<Int> = check(asIntStructure())
+
+    public override operator fun Nd4jArrayStructure<Int>.plus(arg: Int): Nd4jArrayStructure<Int> {
+        check(this)
+        return ndArray.add(arg).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<Int>.minus(arg: Int): Nd4jArrayStructure<Int> {
+        check(this)
+        return ndArray.sub(arg).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<Int>.times(arg: Int): Nd4jArrayStructure<Int> {
+        check(this)
+        return ndArray.mul(arg).wrap()
+    }
+
+    public override operator fun Int.minus(arg: Nd4jArrayStructure<Int>): Nd4jArrayStructure<Int> {
+        check(arg)
+        return arg.ndArray.rsub(this).wrap()
+    }
+}
+
+/**
+ * Represents [NDRing] over [Nd4jArrayStructure] of [Long].
+ */
+public class LongNd4jArrayRing(public override val shape: IntArray) : Nd4jArrayRing<Long, LongRing> {
+    public override val elementContext: LongRing
+        get() = LongRing
+
+    public override fun INDArray.wrap(): Nd4jArrayStructure<Long> = check(asLongStructure())
+
+    public override operator fun Nd4jArrayStructure<Long>.plus(arg: Long): Nd4jArrayStructure<Long> {
+        check(this)
+        return ndArray.add(arg).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<Long>.minus(arg: Long): Nd4jArrayStructure<Long> {
+        check(this)
+        return ndArray.sub(arg).wrap()
+    }
+
+    public override operator fun Nd4jArrayStructure<Long>.times(arg: Long): Nd4jArrayStructure<Long> {
+        check(this)
+        return ndArray.mul(arg).wrap()
+    }
+
+    public override operator fun Long.minus(arg: Nd4jArrayStructure<Long>): Nd4jArrayStructure<Long> {
+        check(arg)
+        return arg.ndArray.rsub(this).wrap()
+    }
+}

kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayIterator.kt

@@ -0,0 +1,62 @@
+package kscience.kmath.nd4j
+
+import org.nd4j.linalg.api.ndarray.INDArray
+import org.nd4j.linalg.api.shape.Shape
+
+private class Nd4jArrayIndicesIterator(private val iterateOver: INDArray) : Iterator<IntArray> {
+    private var i: Int = 0
+
+    override fun hasNext(): Boolean = i < iterateOver.length()
+
+    override fun next(): IntArray {
+        val la = if (iterateOver.ordering() == 'c')
+            Shape.ind2subC(iterateOver, i++.toLong())!!
+        else
+            Shape.ind2sub(iterateOver, i++.toLong())!!
+
+        return la.toIntArray()
+    }
+}
+
+internal fun INDArray.indicesIterator(): Iterator<IntArray> = Nd4jArrayIndicesIterator(this)
+
+private sealed class Nd4jArrayIteratorBase<T>(protected val iterateOver: INDArray) : Iterator<Pair<IntArray, T>> {
+    private var i: Int = 0
+
+    final override fun hasNext(): Boolean = i < iterateOver.length()
+
+    abstract fun getSingle(indices: LongArray): T
+
+    final override fun next(): Pair<IntArray, T> {
+        val la = if (iterateOver.ordering() == 'c')
+            Shape.ind2subC(iterateOver, i++.toLong())!!
+        else
+            Shape.ind2sub(iterateOver, i++.toLong())!!
+
+        return la.toIntArray() to getSingle(la)
+    }
+}
+
+private class Nd4jArrayRealIterator(iterateOver: INDArray) : Nd4jArrayIteratorBase<Double>(iterateOver) {
+    override fun getSingle(indices: LongArray): Double = iterateOver.getDouble(*indices)
+}
+
+internal fun INDArray.realIterator(): Iterator<Pair<IntArray, Double>> = Nd4jArrayRealIterator(this)
+
+private class Nd4jArrayLongIterator(iterateOver: INDArray) : Nd4jArrayIteratorBase<Long>(iterateOver) {
+    override fun getSingle(indices: LongArray) = iterateOver.getLong(*indices)
+}
+
+internal fun INDArray.longIterator(): Iterator<Pair<IntArray, Long>> = Nd4jArrayLongIterator(this)
+
+private class Nd4jArrayIntIterator(iterateOver: INDArray) : Nd4jArrayIteratorBase<Int>(iterateOver) {
+    override fun getSingle(indices: LongArray) = iterateOver.getInt(*indices.toIntArray())
+}
+
+internal fun INDArray.intIterator(): Iterator<Pair<IntArray, Int>> = Nd4jArrayIntIterator(this)
+
+private class Nd4jArrayFloatIterator(iterateOver: INDArray) : Nd4jArrayIteratorBase<Float>(iterateOver) {
+    override fun getSingle(indices: LongArray) = iterateOver.getFloat(*indices)
+}
+
+internal fun INDArray.floatIterator(): Iterator<Pair<IntArray, Float>> = Nd4jArrayFloatIterator(this)

kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/Nd4jArrayStructure.kt

@@ -0,0 +1,68 @@
+package kscience.kmath.nd4j
+
+import kscience.kmath.structures.MutableNDStructure
+import kscience.kmath.structures.NDStructure
+import org.nd4j.linalg.api.ndarray.INDArray
+
+/**
+ * Represents a [NDStructure] wrapping an [INDArray] object.
+ *
+ * @param T the type of items.
+ */
+public sealed class Nd4jArrayStructure<T> : MutableNDStructure<T> {
+    /**
+     * The wrapped [INDArray].
+     */
+    public abstract val ndArray: INDArray
+
+    public override val shape: IntArray
+        get() = ndArray.shape().toIntArray()
+
+    internal abstract fun elementsIterator(): Iterator<Pair<IntArray, T>>
+    internal fun indicesIterator(): Iterator<IntArray> = ndArray.indicesIterator()
+    public override fun elements(): Sequence<Pair<IntArray, T>> = Sequence(::elementsIterator)
+}
+
+private data class Nd4jArrayIntStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Int>() {
+    override fun elementsIterator(): Iterator<Pair<IntArray, Int>> = ndArray.intIterator()
+    override fun get(index: IntArray): Int = ndArray.getInt(*index)
+    override fun set(index: IntArray, value: Int): Unit = run { ndArray.putScalar(index, value) }
+}
+
+/**
+ * Wraps this [INDArray] to [Nd4jArrayStructure].
+ */
+public fun INDArray.asIntStructure(): Nd4jArrayStructure<Int> = Nd4jArrayIntStructure(this)
+
+private data class Nd4jArrayLongStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Long>() {
+    override fun elementsIterator(): Iterator<Pair<IntArray, Long>> = ndArray.longIterator()
+    override fun get(index: IntArray): Long = ndArray.getLong(*index.toLongArray())
+    override fun set(index: IntArray, value: Long): Unit = run { ndArray.putScalar(index, value.toDouble()) }
+}
+
+/**
+ * Wraps this [INDArray] to [Nd4jArrayStructure].
+ */
+public fun INDArray.asLongStructure(): Nd4jArrayStructure<Long> = Nd4jArrayLongStructure(this)
+
+private data class Nd4jArrayRealStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Double>() {
+    override fun elementsIterator(): Iterator<Pair<IntArray, Double>> = ndArray.realIterator()
+    override fun get(index: IntArray): Double = ndArray.getDouble(*index)
+    override fun set(index: IntArray, value: Double): Unit = run { ndArray.putScalar(index, value) }
+}
+
+/**
+ * Wraps this [INDArray] to [Nd4jArrayStructure].
+ */
+public fun INDArray.asRealStructure(): Nd4jArrayStructure<Double> = Nd4jArrayRealStructure(this)
+
+private data class Nd4jArrayFloatStructure(override val ndArray: INDArray) : Nd4jArrayStructure<Float>() {
+    override fun elementsIterator(): Iterator<Pair<IntArray, Float>> = ndArray.floatIterator()
+    override fun get(index: IntArray): Float = ndArray.getFloat(*index)
+    override fun set(index: IntArray, value: Float): Unit = run { ndArray.putScalar(index, value) }
+}
+
+/**
+ * Wraps this [INDArray] to [Nd4jArrayStructure].
+ */
+public fun INDArray.asFloatStructure(): Nd4jArrayStructure<Float> = Nd4jArrayFloatStructure(this)

kmath-nd4j/src/main/kotlin/kscience.kmath.nd4j/arrays.kt

@@ -0,0 +1,4 @@
+package kscience.kmath.nd4j
+
+internal fun IntArray.toLongArray(): LongArray = LongArray(size) { this[it].toLong() }
+internal fun LongArray.toIntArray(): IntArray = IntArray(size) { this[it].toInt() }

kmath-nd4j/src/test/kotlin/kscience/kmath/nd4j/Nd4jArrayAlgebraTest.kt

@@ -0,0 +1,42 @@
+package kscience.kmath.nd4j
+
+import org.nd4j.linalg.factory.Nd4j
+import kscience.kmath.operations.invoke
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.fail
+
+internal class Nd4jArrayAlgebraTest {
+    @Test
+    fun testProduce() {
+        val res = (RealNd4jArrayField(intArrayOf(2, 2))) { produce { it.sum().toDouble() } }
+        val expected = (Nd4j.create(2, 2) ?: fail()).asRealStructure()
+        expected[intArrayOf(0, 0)] = 0.0
+        expected[intArrayOf(0, 1)] = 1.0
+        expected[intArrayOf(1, 0)] = 1.0
+        expected[intArrayOf(1, 1)] = 2.0
+        assertEquals(expected, res)
+    }
+
+    @Test
+    fun testMap() {
+        val res = (IntNd4jArrayRing(intArrayOf(2, 2))) { map(one) { it + it * 2 } }
+        val expected = (Nd4j.create(2, 2) ?: fail()).asIntStructure()
+        expected[intArrayOf(0, 0)] = 3
+        expected[intArrayOf(0, 1)] = 3
+        expected[intArrayOf(1, 0)] = 3
+        expected[intArrayOf(1, 1)] = 3
+        assertEquals(expected, res)
+    }
+
+    @Test
+    fun testAdd() {
+        val res = (IntNd4jArrayRing(intArrayOf(2, 2))) { one + 25 }
+        val expected = (Nd4j.create(2, 2) ?: fail()).asIntStructure()
+        expected[intArrayOf(0, 0)] = 26
+        expected[intArrayOf(0, 1)] = 26
+        expected[intArrayOf(1, 0)] = 26
+        expected[intArrayOf(1, 1)] = 26
+        assertEquals(expected, res)
+    }
+}

kmath-nd4j/src/test/kotlin/kscience/kmath/nd4j/Nd4jArrayStructureTest.kt

@@ -0,0 +1,72 @@
+package kscience.kmath.nd4j
+
+import kscience.kmath.structures.get
+import org.nd4j.linalg.factory.Nd4j
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertNotEquals
+import kotlin.test.fail
+
+internal class Nd4jArrayStructureTest {
+    @Test
+    fun testElements() {
+        val nd = Nd4j.create(doubleArrayOf(1.0, 2.0, 3.0))!!
+        val struct = nd.asRealStructure()
+        val res = struct.elements().map(Pair<IntArray, Double>::second).toList()
+        assertEquals(listOf(1.0, 2.0, 3.0), res)
+    }
+
+    @Test
+    fun testShape() {
+        val nd = Nd4j.rand(10, 2, 3, 6) ?: fail()
+        val struct = nd.asRealStructure()
+        assertEquals(intArrayOf(10, 2, 3, 6).toList(), struct.shape.toList())
+    }
+
+    @Test
+    fun testEquals() {
+        val nd1 = Nd4j.create(doubleArrayOf(1.0, 2.0, 3.0)) ?: fail()
+        val struct1 = nd1.asRealStructure()
+        assertEquals(struct1, struct1)
+        assertNotEquals(struct1 as Any?, null)
+        val nd2 = Nd4j.create(doubleArrayOf(1.0, 2.0, 3.0)) ?: fail()
+        val struct2 = nd2.asRealStructure()
+        assertEquals(struct1, struct2)
+        assertEquals(struct2, struct1)
+        val nd3 = Nd4j.create(doubleArrayOf(1.0, 2.0, 3.0)) ?: fail()
+        val struct3 = nd3.asRealStructure()
+        assertEquals(struct2, struct3)
+        assertEquals(struct1, struct3)
+    }
+
+    @Test
+    fun testHashCode() {
+        val nd1 = Nd4j.create(doubleArrayOf(1.0, 2.0, 3.0))?:fail()
+        val struct1 = nd1.asRealStructure()
+        val nd2 = Nd4j.create(doubleArrayOf(1.0, 2.0, 3.0))?:fail()
+        val struct2 = nd2.asRealStructure()
+        assertEquals(struct1.hashCode(), struct2.hashCode())
+    }
+
+    @Test
+    fun testDimension() {
+        val nd = Nd4j.rand(8, 16, 3, 7, 1)!!
+        val struct = nd.asFloatStructure()
+        assertEquals(5, struct.dimension)
+    }
+
+    @Test
+    fun testGet() {
+        val nd = Nd4j.rand(10, 2, 3, 6)?:fail()
+        val struct = nd.asIntStructure()
+        assertEquals(nd.getInt(0, 0, 0, 0), struct[0, 0, 0, 0])
+    }
+
+    @Test
+    fun testSet() {
+        val nd = Nd4j.rand(17, 12, 4, 8)!!
+        val struct = nd.asLongStructure()
+        struct[intArrayOf(1, 2, 3, 4)] = 777
+        assertEquals(777, struct[1, 2, 3, 4])
+    }
+}

kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/Distribution.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import kotlinx.coroutines.flow.first
 import kscience.kmath.chains.Chain

kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/FactorizedDistribution.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import kscience.kmath.chains.Chain
 import kscience.kmath.chains.SimpleChain

kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/Fitting.kt

@@ -0,0 +1,63 @@
+package kscience.kmath.stat
+
+import kscience.kmath.expressions.*
+import kscience.kmath.operations.ExtendedField
+import kscience.kmath.structures.Buffer
+import kscience.kmath.structures.indices
+import kotlin.math.pow
+
+public object Fitting {
+
+    /**
+     * Generate a chi squared expression from given x-y-sigma data and inline model. Provides automatic differentiation
+     */
+    public fun <T : Any, I : Any, A> chiSquared(
+        autoDiff: AutoDiffProcessor<T, I, A, Expression<T>>,
+        x: Buffer<T>,
+        y: Buffer<T>,
+        yErr: Buffer<T>,
+        model: A.(I) -> I,
+    ): DifferentiableExpression<T, Expression<T>> where A : ExtendedField<I>, A : ExpressionAlgebra<T, I> {
+        require(x.size == y.size) { "X and y buffers should be of the same size" }
+        require(y.size == yErr.size) { "Y and yErr buffer should of the same size" }
+
+        return autoDiff.process {
+            var sum = zero
+
+            x.indices.forEach {
+                val xValue = const(x[it])
+                val yValue = const(y[it])
+                val yErrValue = const(yErr[it])
+                val modelValue = model(xValue)
+                sum += ((yValue - modelValue) / yErrValue).pow(2)
+            }
+
+            sum
+        }
+    }
+
+    /**
+     * Generate a chi squared expression from given x-y-sigma model represented by an expression. Does not provide derivatives
+     */
+    public fun chiSquared(
+        x: Buffer<Double>,
+        y: Buffer<Double>,
+        yErr: Buffer<Double>,
+        model: Expression<Double>,
+        xSymbol: Symbol = StringSymbol("x"),
+    ): Expression<Double> {
+        require(x.size == y.size) { "X and y buffers should be of the same size" }
+        require(y.size == yErr.size) { "Y and yErr buffer should of the same size" }
+
+        return Expression { arguments ->
+            x.indices.sumByDouble {
+                val xValue = x[it]
+                val yValue = y[it]
+                val yErrValue = yErr[it]
+                val modifiedArgs = arguments + (xSymbol to xValue)
+                val modelValue = model(modifiedArgs)
+                ((yValue - modelValue) / yErrValue).pow(2)
+            }
+        }
+    }
+}

kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/OptimizationProblem.kt

@@ -0,0 +1,88 @@
+package kscience.kmath.stat
+
+import kscience.kmath.expressions.DifferentiableExpression
+import kscience.kmath.expressions.Expression
+import kscience.kmath.expressions.Symbol
+
+public interface OptimizationFeature
+
+public class OptimizationResult<T>(
+    public val point: Map<Symbol, T>,
+    public val value: T,
+    public val features: Set<OptimizationFeature> = emptySet(),
+) {
+    override fun toString(): String {
+        return "OptimizationResult(point=$point, value=$value)"
+    }
+}
+
+public operator fun <T> OptimizationResult<T>.plus(
+    feature: OptimizationFeature,
+): OptimizationResult<T> = OptimizationResult(point, value, features + feature)
+
+/**
+ * A configuration builder for optimization problem
+ */
+public interface OptimizationProblem<T : Any> {
+    /**
+     * Define the initial guess for the optimization problem
+     */
+    public fun initialGuess(map: Map<Symbol, T>)
+
+    /**
+     * Set an objective function expression
+     */
+    public fun expression(expression: Expression<T>)
+
+    /**
+     * Set a differentiable expression as objective function as function and gradient provider
+     */
+    public fun diffExpression(expression: DifferentiableExpression<T, Expression<T>>)
+
+    /**
+     * Update the problem from previous optimization run
+     */
+    public fun update(result: OptimizationResult<T>)
+
+    /**
+     * Make an optimization run
+     */
+    public fun optimize(): OptimizationResult<T>
+}
+
+public fun interface OptimizationProblemFactory<T : Any, out P : OptimizationProblem<T>> {
+    public fun build(symbols: List<Symbol>): P
+}
+
+public operator fun <T : Any, P : OptimizationProblem<T>> OptimizationProblemFactory<T, P>.invoke(
+    symbols: List<Symbol>,
+    block: P.() -> Unit,
+): P = build(symbols).apply(block)
+
+/**
+ * Optimize expression without derivatives using specific [OptimizationProblemFactory]
+ */
+public fun <T : Any, F : OptimizationProblem<T>> Expression<T>.optimizeWith(
+    factory: OptimizationProblemFactory<T, F>,
+    vararg symbols: Symbol,
+    configuration: F.() -> Unit,
+): OptimizationResult<T> {
+    require(symbols.isNotEmpty()) { "Must provide a list of symbols for optimization" }
+    val problem = factory(symbols.toList(),configuration)
+    problem.expression(this)
+    return problem.optimize()
+}
+
+/**
+ * Optimize differentiable expression using specific [OptimizationProblemFactory]
+ */
+public fun <T : Any, F : OptimizationProblem<T>>  DifferentiableExpression<T, Expression<T>>.optimizeWith(
+    factory: OptimizationProblemFactory<T, F>,
+    vararg symbols: Symbol,
+    configuration: F.() -> Unit,
+): OptimizationResult<T> {
+    require(symbols.isNotEmpty()) { "Must provide a list of symbols for optimization" }
+    val problem = factory(symbols.toList(), configuration)
+    problem.diffExpression(this)
+    return problem.optimize()
+}

kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/RandomChain.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import kscience.kmath.chains.BlockingIntChain
 import kscience.kmath.chains.BlockingRealChain

kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/RandomGenerator.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import kotlin.random.Random
 

kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/SamplerAlgebra.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import kscience.kmath.chains.Chain
 import kscience.kmath.chains.ConstantChain

kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/Statistic.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import kotlinx.coroutines.CoroutineDispatcher
 import kotlinx.coroutines.Dispatchers

kmath-stat/src/commonMain/kotlin/kscience/kmath/stat/UniformDistribution.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import kscience.kmath.chains.Chain
 import kscience.kmath.chains.SimpleChain

kmath-stat/src/jvmMain/kotlin/kscience/kmath/stat/RandomSourceGenerator.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import org.apache.commons.rng.UniformRandomProvider
 import org.apache.commons.rng.simple.RandomSource

kmath-stat/src/jvmTest/kotlin/kscience/kmath/stat/CommonsDistributionsTest.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import kotlinx.coroutines.flow.take
 import kotlinx.coroutines.flow.toList

kmath-stat/src/jvmTest/kotlin/kscience/kmath/stat/SamplerTest.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import kotlinx.coroutines.runBlocking
 import kotlin.test.Test

kmath-stat/src/jvmTest/kotlin/kscience/kmath/stat/StatisticTest.kt

@@ -1,4 +1,4 @@
-package kscience.kmath.prob
+package kscience.kmath.stat
 
 import kotlinx.coroutines.flow.drop
 import kotlinx.coroutines.flow.first

kmath-viktor/build.gradle.kts

@@ -1,4 +1,6 @@
-plugins { id("ru.mipt.npm.jvm") }
+plugins {
+    id("ru.mipt.npm.jvm")
+}
 
 description = "Binding for https://github.com/JetBrains-Research/viktor"
 

settings.gradle.kts

@@ -1,17 +1,15 @@
 pluginManagement {
     repositories {
-        mavenLocal()
-        jcenter()
         gradlePluginPortal()
+        jcenter()
         maven("https://dl.bintray.com/kotlin/kotlin-eap")
         maven("https://dl.bintray.com/mipt-npm/kscience")
         maven("https://dl.bintray.com/mipt-npm/dev")
         maven("https://dl.bintray.com/kotlin/kotlinx")
-        maven("https://dl.bintray.com/kotlin/kotlin-dev/")
     }
 
-    val toolsVersion = "0.6.1-dev-1.4.20-M1"
+    val toolsVersion = "0.6.4-dev-1.4.20-M2"
-    val kotlinVersion = "1.4.20-M1"
+    val kotlinVersion = "1.4.20-M2"
 
     plugins {
         id("kotlinx.benchmark") version "0.2.0-dev-20"
@@ -30,16 +28,17 @@ include(
     ":kmath-memory",
     ":kmath-core",
     ":kmath-functions",
-//    ":kmath-io",
     ":kmath-coroutines",
     ":kmath-histograms",
     ":kmath-commons",
     ":kmath-viktor",
-    ":kmath-prob",
+    ":kmath-stat",
+    ":kmath-nd4j",
     ":kmath-dimensions",
     ":kmath-for-real",
     ":kmath-geometry",
     ":kmath-ast",
-    ":examples",
+    ":kmath-ejml",
-    ":kmath-ejml"
+    ":kmath-kotlingrad",
+    ":examples"
 )