CHANGELOG.md

@@ -13,6 +13,7 @@
 - Buffer factories for primitives moved to MutableBuffer.Companion
 - NDStructure and NDAlgebra to StructureND and AlgebraND respectively
 - Real -> Double
+- DataSets are moved from functions to core
 
 ### Deprecated
 

README.md

@@ -1,9 +1,8 @@
 [![JetBrains Research](https://jb.gg/badges/research.svg)](https://confluence.jetbrains.com/display/ALL/JetBrains+on+GitHub)
 [![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382)
-
 ![Gradle build](https://github.com/mipt-npm/kmath/workflows/Gradle%20build/badge.svg)
-
+[![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22)
-[![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22%20AND%20a:%22kmath-core%22)
+[![Space](https://img.shields.io/maven-metadata/v?label=Space&metadataUrl=https%3A%2F%2Fmaven.pkg.jetbrains.space%2Fmipt-npm%2Fp%2Fsci%2Fmaven%2Fkscience%2Fkmath%2Fkmath-core%2Fmaven-metadata.xml)](https://maven.pkg.jetbrains.space/mipt-npm/p/sci/maven/space/kscience/)
 
 # KMath
 
@@ -147,6 +146,12 @@ performance calculations to code generation.
 > 
 >
 > **Maturity**: PROTOTYPE
+>
+> **Features:**
+> - [ejml-vector](kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt) : The Point implementation using SimpleMatrix.
+> - [ejml-matrix](kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt) : The Matrix implementation using SimpleMatrix.
+> - [ejml-linear-space](kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt) : The LinearSpace implementation using SimpleMatrix.
+
 <hr/>
 
 * ### [kmath-for-real](kmath-for-real)

build.gradle.kts

@@ -1,4 +1,6 @@
+import org.jetbrains.dokka.gradle.DokkaTask
 import ru.mipt.npm.gradle.KSciencePublishingPlugin
+import java.net.URL
 
 plugins {
     id("ru.mipt.npm.gradle.project")
@@ -10,21 +12,38 @@ allprojects {
         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/mipt-npm/kscience")
         maven("https://jitpack.io")
         maven("http://logicrunch.research.it.uu.se/maven/")
         mavenCentral()
     }
 
     group = "space.kscience"
-    version = "0.3.0-dev-3"
+    version = "0.3.0-dev-4"
 }
 
 subprojects {
     if (name.startsWith("kmath")) apply<KSciencePublishingPlugin>()
+
+    afterEvaluate {
+        tasks.withType<DokkaTask> {
+            dokkaSourceSets.all {
+                val readmeFile = File(this@subprojects.projectDir, "./README.md")
+                if (readmeFile.exists())
+                    includes.setFrom(includes + readmeFile.absolutePath)
+
+                arrayOf(
+                    "http://ejml.org/javadoc/",
+                    "https://commons.apache.org/proper/commons-math/javadocs/api-3.6.1/",
+                    "https://deeplearning4j.org/api/latest/"
+                ).map { URL("${it}package-list") to URL(it) }.forEach { (a, b) ->
+                    externalDocumentationLink {
+                        packageListUrl.set(a)
+                        url.set(b)
+                    }
+                }
+            }
+        }
+    }
 }
 
 readme {

docs/templates/ARTIFACT-TEMPLATE.md

@@ -1,34 +1,28 @@
-> #### Artifact:
+## Artifact:
->
+
-> This module artifact: `${group}:${name}:${version}`.
+The Maven coordinates of this project are `${group}:${name}:${version}`.
->
+
-> Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/${name}/images/download.svg) ](https://bintray.com/mipt-npm/kscience/${name}/_latestVersion)
+**Gradle:**
->
+```gradle
-> Bintray development version:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/${name}/images/download.svg) ](https://bintray.com/mipt-npm/dev/${name}/_latestVersion)
+repositories {
->
+    maven { url 'https://repo.kotlin.link' }
-> **Gradle:**
+    maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
->
+    maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
-> ```gradle
+}
-> repositories {
+
->     maven { url 'https://repo.kotlin.link' }
+dependencies {
->     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+    implementation '${group}:${name}:${version}'
->     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
+}
-> }
+```
-> 
+**Gradle Kotlin DSL:**
-> dependencies {
+```kotlin
->     implementation '${group}:${name}:${version}'
+repositories {
-> }
+    maven("https://repo.kotlin.link")
-> ```
+    maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
-> **Gradle Kotlin DSL:**
+    maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
->
+}
-> ```kotlin
+
-> repositories {
+dependencies {
->     maven("https://repo.kotlin.link")
+    implementation("${group}:${name}:${version}")
->     maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
+}
->     maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
+```
-> }
-> 
-> dependencies {
->     implementation("${group}:${name}:${version}")
-> }
-> ```

docs/templates/README-TEMPLATE.md

@@ -1,9 +1,8 @@
 [![JetBrains Research](https://jb.gg/badges/research.svg)](https://confluence.jetbrains.com/display/ALL/JetBrains+on+GitHub)
 [![DOI](https://zenodo.org/badge/129486382.svg)](https://zenodo.org/badge/latestdoi/129486382)
-
 ![Gradle build](https://github.com/mipt-npm/kmath/workflows/Gradle%20build/badge.svg)
-
+[![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22)
-[![Maven Central](https://img.shields.io/maven-central/v/space.kscience/kmath-core.svg?label=Maven%20Central)](https://search.maven.org/search?q=g:%22space.kscience%22%20AND%20a:%22kmath-core%22)
+[![Space](https://img.shields.io/maven-metadata/v?label=Space&metadataUrl=https%3A%2F%2Fmaven.pkg.jetbrains.space%2Fmipt-npm%2Fp%2Fsci%2Fmaven%2Fkscience%2Fkmath%2Fkmath-core%2Fmaven-metadata.xml)](https://maven.pkg.jetbrains.space/mipt-npm/p/sci/maven/space/kscience/)
 
 # KMath
 

examples/build.gradle.kts

@@ -106,6 +106,7 @@ kotlin.sourceSets.all {
     with(languageSettings) {
         useExperimentalAnnotation("kotlin.contracts.ExperimentalContracts")
         useExperimentalAnnotation("kotlin.ExperimentalUnsignedTypes")
+        useExperimentalAnnotation("space.kscience.kmath.misc.UnstableKMathAPI")
     }
 }
 

examples/src/benchmarks/kotlin/space/kscience/kmath/benchmarks/ExpressionsInterpretersBenchmark.kt

@@ -9,7 +9,7 @@ import space.kscience.kmath.ast.mstInField
 import space.kscience.kmath.expressions.Expression
 import space.kscience.kmath.expressions.expressionInField
 import space.kscience.kmath.expressions.invoke
-import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.misc.symbol
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.bindSymbol
 import kotlin.random.Random

examples/src/main/kotlin/space/kscience/kmath/ast/kotlingradSupport.kt

@@ -3,8 +3,8 @@ package space.kscience.kmath.ast
 import space.kscience.kmath.asm.compile
 import space.kscience.kmath.expressions.derivative
 import space.kscience.kmath.expressions.invoke
-import space.kscience.kmath.expressions.symbol
 import space.kscience.kmath.kotlingrad.differentiable
+import space.kscience.kmath.misc.symbol
 import space.kscience.kmath.operations.DoubleField
 
 /**

examples/src/main/kotlin/space/kscience/kmath/commons/fit/fitWithAutoDiff.kt

@@ -7,11 +7,14 @@ import kscience.plotly.models.ScatterMode
 import kscience.plotly.models.TraceValues
 import space.kscience.kmath.commons.optimization.chiSquared
 import space.kscience.kmath.commons.optimization.minimize
-import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.misc.symbol
+import space.kscience.kmath.optimization.FunctionOptimization
+import space.kscience.kmath.optimization.OptimizationResult
 import space.kscience.kmath.real.DoubleVector
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
-import space.kscience.kmath.stat.*
+import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.stat.distributions.NormalDistribution
 import space.kscience.kmath.structures.asIterable
 import space.kscience.kmath.structures.toList
 import kotlin.math.pow
@@ -33,10 +36,9 @@ operator fun TraceValues.invoke(vector: DoubleVector) {
 /**
  * Least squares fie with auto-differentiation. Uses `kmath-commons` and `kmath-for-real` modules.
  */
-fun main() {
+suspend fun main() {
-
     //A generator for a normally distributed values
-    val generator = Distribution.normal()
+    val generator = NormalDistribution(2.0, 7.0)
 
     //A chain/flow of random values with the given seed
     val chain = generator.sample(RandomGenerator.default(112667))
@@ -49,7 +51,7 @@ fun main() {
     //Perform an operation on each x value (much more effective, than numpy)
     val y = x.map {
         val value = it.pow(2) + it + 1
-        value + chain.nextDouble() * sqrt(value)
+        value + chain.next() * sqrt(value)
     }
     // this will also work, but less effective:
     // val y = x.pow(2)+ x + 1 + chain.nextDouble()
@@ -58,7 +60,7 @@ fun main() {
     val yErr = y.map { sqrt(it) }//RealVector.same(x.size, sigma)
 
     // compute differentiable chi^2 sum for given model ax^2 + bx + c
-    val chi2 = Fitting.chiSquared(x, y, yErr) { x1 ->
+    val chi2 = FunctionOptimization.chiSquared(x, y, yErr) { x1 ->
         //bind variables to autodiff context
         val a = bind(a)
         val b = bind(b)

examples/src/main/kotlin/space/kscience/kmath/stat/DistributionBenchmark.kt

@@ -1,23 +1,24 @@
-package kscience.kmath.commons.prob
+package space.kscience.kmath.stat
 
 import kotlinx.coroutines.Dispatchers
 import kotlinx.coroutines.async
 import kotlinx.coroutines.runBlocking
-import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler
+import space.kscience.kmath.stat.samplers.GaussianSampler
 import org.apache.commons.rng.simple.RandomSource
-import space.kscience.kmath.stat.*
 import java.time.Duration
 import java.time.Instant
+import org.apache.commons.rng.sampling.distribution.GaussianSampler as CMGaussianSampler
+import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler as CMZigguratNormalizedGaussianSampler
 
-private fun runChain(): Duration {
+private suspend fun runKMathChained(): Duration {
     val generator = RandomGenerator.fromSource(RandomSource.MT, 123L)
-    val normal = Distribution.normal(NormalSamplerMethod.Ziggurat)
+    val normal = GaussianSampler.of(7.0, 2.0)
-    val chain = normal.sample(generator)
+    val chain = normal.sample(generator).blocking()
     val startTime = Instant.now()
     var sum = 0.0
 
     repeat(10000001) { counter ->
-        sum += chain.nextDouble()
+        sum += chain.next()
 
         if (counter % 100000 == 0) {
             val duration = Duration.between(startTime, Instant.now())
@@ -29,9 +30,15 @@ private fun runChain(): Duration {
     return Duration.between(startTime, Instant.now())
 }
 
-private fun runDirect(): Duration {
+private fun runApacheDirect(): Duration {
-    val provider = RandomSource.create(RandomSource.MT, 123L)
+    val rng = RandomSource.create(RandomSource.MT, 123L)
-    val sampler = ZigguratNormalizedGaussianSampler(provider)
+
+    val sampler = CMGaussianSampler.of(
+        CMZigguratNormalizedGaussianSampler.of(rng),
+        7.0,
+        2.0
+    )
+
     val startTime = Instant.now()
     var sum = 0.0
 
@@ -51,11 +58,9 @@ private fun runDirect(): Duration {
 /**
  * Comparing chain sampling performance with direct sampling performance
  */
-fun main() {
+fun main(): Unit = runBlocking(Dispatchers.Default) {
-    runBlocking(Dispatchers.Default) {
+    val chainJob = async { runKMathChained() }
-        val chainJob = async { runChain() }
+    val directJob = async { runApacheDirect() }
-        val directJob = async { runDirect() }
+    println("KMath Chained: ${chainJob.await()}")
-        println("Chain: ${chainJob.await()}")
+    println("Apache Direct: ${directJob.await()}")
-        println("Direct: ${directJob.await()}")
-    }
 }

examples/src/main/kotlin/space/kscience/kmath/stat/DistributionDemo.kt

@@ -3,14 +3,15 @@ package space.kscience.kmath.stat
 import kotlinx.coroutines.runBlocking
 import space.kscience.kmath.chains.Chain
 import space.kscience.kmath.chains.collectWithState
+import space.kscience.kmath.stat.distributions.NormalDistribution
 
 /**
- * The state of distribution averager
+ * The state of distribution averager.
  */
 private data class AveragingChainState(var num: Int = 0, var value: Double = 0.0)
 
 /**
- * Averaging
+ * Averaging.
  */
 private fun Chain<Double>.mean(): Chain<Double> = collectWithState(AveragingChainState(), { it.copy() }) { chain ->
     val next = chain.next()
@@ -21,7 +22,7 @@ private fun Chain<Double>.mean(): Chain<Double> = collectWithState(AveragingChai
 
 
 fun main() {
-    val normal = Distribution.normal()
+    val normal = NormalDistribution(0.0, 2.0)
     val chain = normal.sample(RandomGenerator.default).mean()
 
     runBlocking {

gradle.properties

@@ -4,5 +4,5 @@ kotlin.mpp.stability.nowarn=true
 kotlin.native.enableDependencyPropagation=false
 kotlin.parallel.tasks.in.project=true
 org.gradle.configureondemand=true
-org.gradle.jvmargs=-XX:MaxMetaspaceSize=2G
+org.gradle.jvmargs=-XX:MaxMetaspaceSize=9G
 org.gradle.parallel=true

kmath-ast/README.md

@@ -1,6 +1,6 @@
-# Abstract Syntax Tree Expression Representation and Operations (`kmath-ast`)
+# Module kmath-ast
 
-This subproject implements the following features:
+Abstract syntax tree expression representation and related optimizations.
 
  - [expression-language](src/jvmMain/kotlin/space/kscience/kmath/ast/parser.kt) : Expression language and its parser
  - [mst](src/commonMain/kotlin/space/kscience/kmath/ast/MST.kt) : MST (Mathematical Syntax Tree) as expression language's syntax intermediate representation
@@ -10,40 +10,34 @@ This subproject implements the following features:
  - [mst-js-codegen](src/jsMain/kotlin/space/kscience/kmath/estree/estree.kt) : Dynamic MST to JS compiler
 
 
-> #### Artifact:
+## Artifact:
->
+
-> This module artifact: `space.kscience:kmath-ast:0.3.0-dev-3`.
+The Maven coordinates of this project are `space.kscience:kmath-ast:0.3.0-dev-3`.
->
+
-> Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-ast/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-ast/_latestVersion)
+**Gradle:**
->
+```gradle
-> Bintray development version:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-ast/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-ast/_latestVersion)
+repositories {
->
+    maven { url 'https://repo.kotlin.link' }
-> **Gradle:**
+    maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
->
+    maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
-> ```gradle
+}
-> repositories {
+
->     maven { url 'https://repo.kotlin.link' }
+dependencies {
->     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+    implementation 'space.kscience:kmath-ast:0.3.0-dev-3'
->     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
+}
-> }
+```
-> 
+**Gradle Kotlin DSL:**
-> dependencies {
+```kotlin
->     implementation 'space.kscience:kmath-ast:0.3.0-dev-3'
+repositories {
-> }
+    maven("https://repo.kotlin.link")
-> ```
+    maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
-> **Gradle Kotlin DSL:**
+    maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
->
+}
-> ```kotlin
+
-> repositories {
+dependencies {
->     maven("https://repo.kotlin.link")
+    implementation("space.kscience:kmath-ast:0.3.0-dev-3")
->     maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
+}
->     maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
+```
-> }
-> 
-> dependencies {
->     implementation("space.kscience:kmath-ast:0.3.0-dev-3")
-> }
-> ```
 
 ## Dynamic expression code generation
 

kmath-ast/docs/README-TEMPLATE.md

@@ -1,6 +1,6 @@
-# Abstract Syntax Tree Expression Representation and Operations (`kmath-ast`)
+# Module kmath-ast
 
-This subproject implements the following features:
+Abstract syntax tree expression representation and related optimizations.
 
 ${features}
 

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

@@ -1,6 +1,8 @@
 package space.kscience.kmath.ast
 
 import space.kscience.kmath.expressions.*
+import space.kscience.kmath.misc.StringSymbol
+import space.kscience.kmath.misc.Symbol
 import space.kscience.kmath.operations.*
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract

kmath-ast/src/jsMain/kotlin/space/kscience/kmath/estree/internal/ESTreeBuilder.kt

@@ -3,7 +3,7 @@ package space.kscience.kmath.estree.internal
 import space.kscience.kmath.estree.internal.astring.generate
 import space.kscience.kmath.estree.internal.estree.*
 import space.kscience.kmath.expressions.Expression
-import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.misc.Symbol
 
 internal class ESTreeBuilder<T>(val bodyCallback: ESTreeBuilder<T>.() -> BaseExpression) {
     private class GeneratedExpression<T>(val executable: dynamic, val constants: Array<dynamic>) : Expression<T> {

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

@@ -2,8 +2,8 @@
 
 package space.kscience.kmath.asm.internal
 
-import space.kscience.kmath.expressions.StringSymbol
+import space.kscience.kmath.misc.StringSymbol
-import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.misc.Symbol
 
 /**
  * Gets value with given [key] or throws [NoSuchElementException] whenever it is not present.

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

@@ -2,6 +2,8 @@ package space.kscience.kmath.commons.expressions
 
 import org.apache.commons.math3.analysis.differentiation.DerivativeStructure
 import space.kscience.kmath.expressions.*
+import space.kscience.kmath.misc.StringSymbol
+import space.kscience.kmath.misc.Symbol
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.ExtendedField
 import space.kscience.kmath.operations.NumbersAddOperations

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/linear/CMMatrix.kt

@@ -3,6 +3,7 @@ package space.kscience.kmath.commons.linear
 import org.apache.commons.math3.linear.*
 import space.kscience.kmath.linear.*
 import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.nd.StructureFeature
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.structures.DoubleBuffer
 import kotlin.reflect.KClass
@@ -89,7 +90,7 @@ public object CMLinearSpace : LinearSpace<Double, DoubleField> {
         v * this
 
     @UnstableKMathAPI
-    override fun <F : Any> getFeature(structure: Matrix<Double>, type: KClass<F>): F? {
+    override fun <F : StructureFeature> getFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
         //Return the feature if it is intrinsic to the structure
         structure.getFeature(type)?.let { return it }
 

kmath-commons/src/main/kotlin/space/kscience/kmath/commons/optimization/CMOptimization.kt

@@ -9,22 +9,36 @@ import org.apache.commons.math3.optim.nonlinear.scalar.gradient.NonLinearConjuga
 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 space.kscience.kmath.expressions.*
+import space.kscience.kmath.expressions.DifferentiableExpression
-import space.kscience.kmath.stat.OptimizationFeature
+import space.kscience.kmath.expressions.Expression
-import space.kscience.kmath.stat.OptimizationProblem
+import space.kscience.kmath.expressions.SymbolIndexer
-import space.kscience.kmath.stat.OptimizationProblemFactory
+import space.kscience.kmath.expressions.derivative
-import space.kscience.kmath.stat.OptimizationResult
+import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.optimization.*
 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>) :
+@OptIn(UnstableKMathAPI::class)
-    OptimizationProblem<Double>, SymbolIndexer, OptimizationFeature {
+public class CMOptimization(
+    override val symbols: List<Symbol>,
+) : FunctionOptimization<Double>, NoDerivFunctionOptimization<Double>, SymbolIndexer, OptimizationFeature {
+
     private val optimizationData: HashMap<KClass<out OptimizationData>, OptimizationData> = HashMap()
-    private var optimizatorBuilder: (() -> MultivariateOptimizer)? = null
+    private var optimizerBuilder: (() -> MultivariateOptimizer)? = null
-    public var convergenceChecker: ConvergenceChecker<PointValuePair> = SimpleValueChecker(DEFAULT_RELATIVE_TOLERANCE,
+    public var convergenceChecker: ConvergenceChecker<PointValuePair> = SimpleValueChecker(
-        DEFAULT_ABSOLUTE_TOLERANCE, DEFAULT_MAX_ITER)
+        DEFAULT_RELATIVE_TOLERANCE,
+        DEFAULT_ABSOLUTE_TOLERANCE,
+        DEFAULT_MAX_ITER
+    )
+
+    override var maximize: Boolean
+        get() = optimizationData[GoalType::class] == GoalType.MAXIMIZE
+        set(value) {
+            optimizationData[GoalType::class] = if (value) GoalType.MAXIMIZE else GoalType.MINIMIZE
+        }
 
     public fun addOptimizationData(data: OptimizationData) {
         optimizationData[data::class] = data
@@ -40,7 +54,7 @@ public class CMOptimizationProblem(override val symbols: List<Symbol>) :
         addOptimizationData(InitialGuess(map.toDoubleArray()))
     }
 
-    public override fun expression(expression: Expression<Double>): Unit {
+    public override fun function(expression: Expression<Double>): Unit {
         val objectiveFunction = ObjectiveFunction {
             val args = it.toMap()
             expression(args)
@@ -48,8 +62,8 @@ public class CMOptimizationProblem(override val symbols: List<Symbol>) :
         addOptimizationData(objectiveFunction)
     }
 
-    public override fun diffExpression(expression: DifferentiableExpression<Double, Expression<Double>>) {
+    public override fun diffFunction(expression: DifferentiableExpression<Double, Expression<Double>>) {
-        expression(expression)
+        function(expression)
         val gradientFunction = ObjectiveFunctionGradient {
             val args = it.toMap()
             DoubleArray(symbols.size) { index ->
@@ -57,8 +71,8 @@ public class CMOptimizationProblem(override val symbols: List<Symbol>) :
             }
         }
         addOptimizationData(gradientFunction)
-        if (optimizatorBuilder == null) {
+        if (optimizerBuilder == null) {
-            optimizatorBuilder = {
+            optimizerBuilder = {
                 NonLinearConjugateGradientOptimizer(
                     NonLinearConjugateGradientOptimizer.Formula.FLETCHER_REEVES,
                     convergenceChecker
@@ -70,8 +84,8 @@ public class CMOptimizationProblem(override val symbols: List<Symbol>) :
     public fun simplex(simplex: AbstractSimplex) {
         addOptimizationData(simplex)
         //Set optimization builder to simplex if it is not present
-        if (optimizatorBuilder == null) {
+        if (optimizerBuilder == null) {
-            optimizatorBuilder = { SimplexOptimizer(convergenceChecker) }
+            optimizerBuilder = { SimplexOptimizer(convergenceChecker) }
         }
     }
 
@@ -84,7 +98,7 @@ public class CMOptimizationProblem(override val symbols: List<Symbol>) :
     }
 
     public fun optimizer(block: () -> MultivariateOptimizer) {
-        optimizatorBuilder = block
+        optimizerBuilder = block
     }
 
     override fun update(result: OptimizationResult<Double>) {
@@ -92,19 +106,19 @@ public class CMOptimizationProblem(override val symbols: List<Symbol>) :
     }
 
     override fun optimize(): OptimizationResult<Double> {
-        val optimizer = optimizatorBuilder?.invoke() ?: error("Optimizer not defined")
+        val optimizer = optimizerBuilder?.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 companion object : OptimizationProblemFactory<Double, CMOptimization> {
         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)
+        override fun build(symbols: List<Symbol>): CMOptimization = CMOptimization(symbols)
     }
 }
 
-public fun CMOptimizationProblem.initialGuess(vararg pairs: Pair<Symbol, Double>): Unit = initialGuess(pairs.toMap())
+public fun CMOptimization.initialGuess(vararg pairs: Pair<Symbol, Double>): Unit = initialGuess(pairs.toMap())
-public fun CMOptimizationProblem.simplexSteps(vararg pairs: Pair<Symbol, Double>): Unit = simplexSteps(pairs.toMap())
+public fun CMOptimization.simplexSteps(vararg pairs: Pair<Symbol, Double>): Unit = simplexSteps(pairs.toMap())

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

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

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

@@ -1,6 +1,10 @@
 package space.kscience.kmath.commons.expressions
 
-import space.kscience.kmath.expressions.*
+import space.kscience.kmath.expressions.binding
+import space.kscience.kmath.expressions.derivative
+import space.kscience.kmath.expressions.invoke
+import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.misc.symbol
 import kotlin.contracts.InvocationKind
 import kotlin.contracts.contract
 import kotlin.test.Test

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

@@ -1,13 +1,13 @@
 package space.kscience.kmath.commons.optimization
 
-import org.junit.jupiter.api.Test
+import kotlinx.coroutines.runBlocking
 import space.kscience.kmath.commons.expressions.DerivativeStructureExpression
-import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.misc.symbol
-import space.kscience.kmath.stat.Distribution
+import space.kscience.kmath.optimization.FunctionOptimization
-import space.kscience.kmath.stat.Fitting
 import space.kscience.kmath.stat.RandomGenerator
-import space.kscience.kmath.stat.normal
+import space.kscience.kmath.stat.distributions.NormalDistribution
 import kotlin.math.pow
+import kotlin.test.Test
 
 internal class OptimizeTest {
     val x by symbol
@@ -34,28 +34,29 @@ internal class OptimizeTest {
             simplexSteps(x to 2.0, y to 0.5)
             //this sets simplex optimizer
         }
+
         println(result.point)
         println(result.value)
     }
 
     @Test
-    fun testCmFit() {
+    fun testCmFit() = runBlocking {
         val a by symbol
         val b by symbol
         val c by symbol
 
         val sigma = 1.0
-        val generator = Distribution.normal(0.0, sigma)
+        val generator = NormalDistribution(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()
+            it.pow(2) + it + 1 + chain.next()
         }
 
         val yErr = List(x.size) { sigma }
 
-        val chi2 = Fitting.chiSquared(x, y, yErr) { x1 ->
+        val chi2 = FunctionOptimization.chiSquared(x, y, yErr) { x1 ->
             val cWithDefault = bindSymbolOrNull(c) ?: one
             bind(a) * x1.pow(2) + bind(b) * x1 + cWithDefault
         }
@@ -64,5 +65,4 @@ internal class OptimizeTest {
         println(result)
         println("Chi2/dof = ${result.value / (x.size - 3)}")
     }
-
 }

kmath-complex/README.md

@@ -1,42 +1,36 @@
-# The Core Module (`kmath-core`)
+# Module kmath-complex
 
-Complex and hypercomplex number systems in KMath:
+Complex and hypercomplex number systems in KMath.
 
  - [complex](src/commonMain/kotlin/space/kscience/kmath/complex/Complex.kt) : Complex Numbers
  - [quaternion](src/commonMain/kotlin/space/kscience/kmath/complex/Quaternion.kt) : Quaternions
 
 
-> #### Artifact:
+## Artifact:
->
+
-> This module artifact: `space.kscience:kmath-complex:0.3.0-dev-3`.
+The Maven coordinates of this project are `space.kscience:kmath-complex:0.3.0-dev-3`.
->
+
-> Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-complex/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-complex/_latestVersion)
+**Gradle:**
->
+```gradle
-> Bintray development version:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-complex/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-complex/_latestVersion)
+repositories {
->
+    maven { url 'https://repo.kotlin.link' }
-> **Gradle:**
+    maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
->
+    maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
-> ```gradle
+}
-> repositories {
+
->     maven { url 'https://repo.kotlin.link' }
+dependencies {
->     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+    implementation 'space.kscience:kmath-complex:0.3.0-dev-3'
->     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
+}
-> }
+```
-> 
+**Gradle Kotlin DSL:**
-> dependencies {
+```kotlin
->     implementation 'space.kscience:kmath-complex:0.3.0-dev-3'
+repositories {
-> }
+    maven("https://repo.kotlin.link")
-> ```
+    maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
-> **Gradle Kotlin DSL:**
+    maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
->
+}
-> ```kotlin
+
-> repositories {
+dependencies {
->     maven("https://repo.kotlin.link")
+    implementation("space.kscience:kmath-complex:0.3.0-dev-3")
->     maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
+}
->     maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
+```
-> }
-> 
-> dependencies {
->     implementation("space.kscience:kmath-complex:0.3.0-dev-3")
-> }
-> ```

kmath-complex/docs/README-TEMPLATE.md

@@ -1,6 +1,6 @@
-# The Core Module (`kmath-core`)
+# Module kmath-complex
 
-Complex and hypercomplex number systems in KMath:
+Complex and hypercomplex number systems in KMath.
 
 ${features}
 

kmath-complex/src/commonTest/kotlin/space/kscience/kmath/complex/ExpressionFieldForComplexTest.kt

@@ -3,7 +3,7 @@ package space.kscience.kmath.complex
 import space.kscience.kmath.expressions.FunctionalExpressionField
 import space.kscience.kmath.expressions.bindSymbol
 import space.kscience.kmath.expressions.invoke
-import space.kscience.kmath.expressions.symbol
+import space.kscience.kmath.misc.symbol
 import kotlin.test.Test
 import kotlin.test.assertEquals
 

kmath-core/README.md

@@ -1,6 +1,6 @@
-# The Core Module (`kmath-core`)
+# Module kmath-core
 
-The core features of KMath:
+The core interfaces of KMath.
 
  - [algebras](src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt) : Algebraic structures like rings, spaces and fields.
  - [nd](src/commonMain/kotlin/space/kscience/kmath/structures/StructureND.kt) : Many-dimensional structures and operations on them.
@@ -13,37 +13,31 @@ performance calculations to code generation.
  - [autodif](src/commonMain/kotlin/space/kscience/kmath/expressions/SimpleAutoDiff.kt) : Automatic differentiation
 
 
-> #### Artifact:
+## Artifact:
->
+
-> This module artifact: `space.kscience:kmath-core:0.3.0-dev-3`.
+The Maven coordinates of this project are `space.kscience:kmath-core:0.3.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)
+**Gradle:**
->
+```gradle
-> Bintray development version:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-core/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-core/_latestVersion)
+repositories {
->
+    maven { url 'https://repo.kotlin.link' }
-> **Gradle:**
+    maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
->
+    maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
-> ```gradle
+}
-> repositories {
+
->     maven { url 'https://repo.kotlin.link' }
+dependencies {
->     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+    implementation 'space.kscience:kmath-core:0.3.0-dev-3'
->     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
+}
-> }
+```
-> 
+**Gradle Kotlin DSL:**
-> dependencies {
+```kotlin
->     implementation 'space.kscience:kmath-core:0.3.0-dev-3'
+repositories {
-> }
+    maven("https://repo.kotlin.link")
-> ```
+    maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
-> **Gradle Kotlin DSL:**
+    maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
->
+}
-> ```kotlin
+
-> repositories {
+dependencies {
->     maven("https://repo.kotlin.link")
+    implementation("space.kscience:kmath-core:0.3.0-dev-3")
->     maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
+}
->     maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
+```
-> }
-> 
-> dependencies {
->     implementation("space.kscience:kmath-core:0.3.0-dev-3")
-> }
-> ```

kmath-core/api/kmath-core.api

@@ -1,3 +1,18 @@
+public final class space/kscience/kmath/data/ColumnarDataKt {
+}
+
+public final class space/kscience/kmath/data/XYColumnarData$DefaultImpls {
+	public static fun get (Lspace/kscience/kmath/data/XYColumnarData;Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/structures/Buffer;
+}
+
+public final class space/kscience/kmath/data/XYColumnarDataKt {
+	public static synthetic fun asXYData$default (Lspace/kscience/kmath/nd/Structure2D;IIILjava/lang/Object;)Lspace/kscience/kmath/data/XYColumnarData;
+}
+
+public final class space/kscience/kmath/data/XYZColumnarData$DefaultImpls {
+	public static fun get (Lspace/kscience/kmath/data/XYZColumnarData;Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/structures/Buffer;
+}
+
 public abstract interface class space/kscience/kmath/domains/Domain {
 	public abstract fun contains (Lspace/kscience/kmath/structures/Buffer;)Z
 	public abstract fun getDimension ()I
@@ -14,7 +29,7 @@ public class space/kscience/kmath/expressions/AutoDiffValue {
 
 public final class space/kscience/kmath/expressions/DerivationResult {
 	public fun <init> (Ljava/lang/Object;Ljava/util/Map;Lspace/kscience/kmath/operations/Field;)V
-	public final fun derivative (Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
+	public final fun derivative (Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object;
 	public final fun div ()Ljava/lang/Object;
 	public final fun getContext ()Lspace/kscience/kmath/operations/Field;
 	public final fun getValue ()Ljava/lang/Object;
@@ -27,7 +42,7 @@ public abstract interface class space/kscience/kmath/expressions/DifferentiableE
 public final class space/kscience/kmath/expressions/DifferentiableExpressionKt {
 	public static final fun derivative (Lspace/kscience/kmath/expressions/DifferentiableExpression;Ljava/lang/String;)Lspace/kscience/kmath/expressions/Expression;
 	public static final fun derivative (Lspace/kscience/kmath/expressions/DifferentiableExpression;Ljava/util/List;)Lspace/kscience/kmath/expressions/Expression;
-	public static final fun derivative (Lspace/kscience/kmath/expressions/DifferentiableExpression;[Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/expressions/Expression;
+	public static final fun derivative (Lspace/kscience/kmath/expressions/DifferentiableExpression;[Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/expressions/Expression;
 }
 
 public abstract interface class space/kscience/kmath/expressions/Expression {
@@ -36,7 +51,7 @@ public abstract interface class space/kscience/kmath/expressions/Expression {
 
 public abstract interface class space/kscience/kmath/expressions/ExpressionAlgebra : space/kscience/kmath/operations/Algebra {
 	public abstract fun bindSymbol (Ljava/lang/String;)Ljava/lang/Object;
-	public abstract fun bindSymbolOrNull (Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
+	public abstract fun bindSymbolOrNull (Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object;
 	public abstract fun const (Ljava/lang/Object;)Ljava/lang/Object;
 }
 
@@ -56,18 +71,17 @@ public final class space/kscience/kmath/expressions/ExpressionBuildersKt {
 }
 
 public final class space/kscience/kmath/expressions/ExpressionKt {
-	public static final fun bindSymbol (Lspace/kscience/kmath/expressions/ExpressionAlgebra;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
+	public static final fun bindSymbol (Lspace/kscience/kmath/expressions/ExpressionAlgebra;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object;
 	public static final fun binding (Lspace/kscience/kmath/expressions/ExpressionAlgebra;)Lkotlin/properties/ReadOnlyProperty;
 	public static final fun callByString (Lspace/kscience/kmath/expressions/Expression;[Lkotlin/Pair;)Ljava/lang/Object;
 	public static final fun callBySymbol (Lspace/kscience/kmath/expressions/Expression;[Lkotlin/Pair;)Ljava/lang/Object;
-	public static final fun getSymbol ()Lkotlin/properties/ReadOnlyProperty;
 	public static final fun invoke (Lspace/kscience/kmath/expressions/Expression;)Ljava/lang/Object;
 }
 
 public abstract class space/kscience/kmath/expressions/FirstDerivativeExpression : space/kscience/kmath/expressions/DifferentiableExpression {
 	public fun <init> ()V
 	public final fun derivativeOrNull (Ljava/util/List;)Lspace/kscience/kmath/expressions/Expression;
-	public abstract fun derivativeOrNull (Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/expressions/Expression;
+	public abstract fun derivativeOrNull (Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/expressions/Expression;
 }
 
 public abstract class space/kscience/kmath/expressions/FunctionalExpressionAlgebra : space/kscience/kmath/expressions/ExpressionAlgebra {
@@ -77,8 +91,8 @@ public abstract class space/kscience/kmath/expressions/FunctionalExpressionAlgeb
 	public fun binaryOperationFunction (Ljava/lang/String;)Lkotlin/jvm/functions/Function2;
 	public synthetic fun bindSymbol (Ljava/lang/String;)Ljava/lang/Object;
 	public fun bindSymbol (Ljava/lang/String;)Lspace/kscience/kmath/expressions/Expression;
-	public synthetic fun bindSymbolOrNull (Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
+	public synthetic fun bindSymbolOrNull (Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object;
-	public fun bindSymbolOrNull (Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/expressions/Expression;
+	public fun bindSymbolOrNull (Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/expressions/Expression;
 	public synthetic fun const (Ljava/lang/Object;)Ljava/lang/Object;
 	public fun const (Ljava/lang/Object;)Lspace/kscience/kmath/expressions/Expression;
 	public final fun getAlgebra ()Lspace/kscience/kmath/operations/Algebra;
@@ -203,7 +217,7 @@ public class space/kscience/kmath/expressions/FunctionalExpressionRing : space/k
 
 public final class space/kscience/kmath/expressions/SimpleAutoDiffExpression : space/kscience/kmath/expressions/FirstDerivativeExpression {
 	public fun <init> (Lspace/kscience/kmath/operations/Field;Lkotlin/jvm/functions/Function1;)V
-	public fun derivativeOrNull (Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/expressions/Expression;
+	public fun derivativeOrNull (Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/expressions/Expression;
 	public final fun getField ()Lspace/kscience/kmath/operations/Field;
 	public final fun getFunction ()Lkotlin/jvm/functions/Function1;
 	public fun invoke (Ljava/util/Map;)Ljava/lang/Object;
@@ -264,8 +278,8 @@ public class space/kscience/kmath/expressions/SimpleAutoDiffField : space/kscien
 	public fun binaryOperationFunction (Ljava/lang/String;)Lkotlin/jvm/functions/Function2;
 	public synthetic fun bindSymbol (Ljava/lang/String;)Ljava/lang/Object;
 	public fun bindSymbol (Ljava/lang/String;)Lspace/kscience/kmath/expressions/AutoDiffValue;
-	public synthetic fun bindSymbolOrNull (Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
+	public synthetic fun bindSymbolOrNull (Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object;
-	public fun bindSymbolOrNull (Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/expressions/AutoDiffValue;
+	public fun bindSymbolOrNull (Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/expressions/AutoDiffValue;
 	public synthetic fun const (Ljava/lang/Object;)Ljava/lang/Object;
 	public fun const (Ljava/lang/Object;)Lspace/kscience/kmath/expressions/AutoDiffValue;
 	public final fun const (Lkotlin/jvm/functions/Function1;)Lspace/kscience/kmath/expressions/AutoDiffValue;
@@ -332,7 +346,7 @@ public final class space/kscience/kmath/expressions/SimpleAutoDiffKt {
 	public static final fun cos (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;)Lspace/kscience/kmath/expressions/AutoDiffValue;
 	public static final fun cosh (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;)Lspace/kscience/kmath/expressions/AutoDiffValue;
 	public static final fun exp (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;)Lspace/kscience/kmath/expressions/AutoDiffValue;
-	public static final fun grad (Lspace/kscience/kmath/expressions/DerivationResult;[Lspace/kscience/kmath/expressions/Symbol;)Lspace/kscience/kmath/structures/Buffer;
+	public static final fun grad (Lspace/kscience/kmath/expressions/DerivationResult;[Lspace/kscience/kmath/misc/Symbol;)Lspace/kscience/kmath/structures/Buffer;
 	public static final fun ln (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;)Lspace/kscience/kmath/expressions/AutoDiffValue;
 	public static final fun pow (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;D)Lspace/kscience/kmath/expressions/AutoDiffValue;
 	public static final fun pow (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;I)Lspace/kscience/kmath/expressions/AutoDiffValue;
@@ -348,79 +362,13 @@ public final class space/kscience/kmath/expressions/SimpleAutoDiffKt {
 	public static final fun tanh (Lspace/kscience/kmath/expressions/SimpleAutoDiffField;Lspace/kscience/kmath/expressions/AutoDiffValue;)Lspace/kscience/kmath/expressions/AutoDiffValue;
 }
 
-public final class space/kscience/kmath/expressions/SimpleSymbolIndexer : space/kscience/kmath/expressions/SymbolIndexer {
-	public static final synthetic fun box-impl (Ljava/util/List;)Lspace/kscience/kmath/expressions/SimpleSymbolIndexer;
-	public static fun constructor-impl (Ljava/util/List;)Ljava/util/List;
-	public fun equals (Ljava/lang/Object;)Z
-	public static fun equals-impl (Ljava/util/List;Ljava/lang/Object;)Z
-	public static final fun equals-impl0 (Ljava/util/List;Ljava/util/List;)Z
-	public fun get (Ljava/util/List;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public fun get (Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/expressions/Symbol;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public fun get (Lspace/kscience/kmath/structures/Buffer;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public fun get ([DLspace/kscience/kmath/expressions/Symbol;)D
-	public fun get ([Ljava/lang/Object;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public static fun get-impl (Ljava/util/List;Ljava/util/List;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public static fun get-impl (Ljava/util/List;Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/expressions/Symbol;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public static fun get-impl (Ljava/util/List;Lspace/kscience/kmath/structures/Buffer;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public static fun get-impl (Ljava/util/List;[DLspace/kscience/kmath/expressions/Symbol;)D
-	public static fun get-impl (Ljava/util/List;[Ljava/lang/Object;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public fun getSymbols ()Ljava/util/List;
-	public fun hashCode ()I
-	public static fun hashCode-impl (Ljava/util/List;)I
-	public fun indexOf (Lspace/kscience/kmath/expressions/Symbol;)I
-	public static fun indexOf-impl (Ljava/util/List;Lspace/kscience/kmath/expressions/Symbol;)I
-	public fun toDoubleArray (Ljava/util/Map;)[D
-	public static fun toDoubleArray-impl (Ljava/util/List;Ljava/util/Map;)[D
-	public fun toList (Ljava/util/Map;)Ljava/util/List;
-	public static fun toList-impl (Ljava/util/List;Ljava/util/Map;)Ljava/util/List;
-	public fun toMap ([D)Ljava/util/Map;
-	public static fun toMap-impl (Ljava/util/List;[D)Ljava/util/Map;
-	public fun toPoint (Ljava/util/Map;Lkotlin/jvm/functions/Function2;)Lspace/kscience/kmath/structures/Buffer;
-	public static fun toPoint-impl (Ljava/util/List;Ljava/util/Map;Lkotlin/jvm/functions/Function2;)Lspace/kscience/kmath/structures/Buffer;
-	public fun toString ()Ljava/lang/String;
-	public static fun toString-impl (Ljava/util/List;)Ljava/lang/String;
-	public final synthetic fun unbox-impl ()Ljava/util/List;
-}
-
-public final class space/kscience/kmath/expressions/StringSymbol : space/kscience/kmath/expressions/Symbol {
-	public static final synthetic fun box-impl (Ljava/lang/String;)Lspace/kscience/kmath/expressions/StringSymbol;
-	public static fun constructor-impl (Ljava/lang/String;)Ljava/lang/String;
-	public fun equals (Ljava/lang/Object;)Z
-	public static fun equals-impl (Ljava/lang/String;Ljava/lang/Object;)Z
-	public static final fun equals-impl0 (Ljava/lang/String;Ljava/lang/String;)Z
-	public fun getIdentity ()Ljava/lang/String;
-	public fun hashCode ()I
-	public static fun hashCode-impl (Ljava/lang/String;)I
-	public fun toString ()Ljava/lang/String;
-	public static fun toString-impl (Ljava/lang/String;)Ljava/lang/String;
-	public final synthetic fun unbox-impl ()Ljava/lang/String;
-}
-
-public abstract interface class space/kscience/kmath/expressions/Symbol {
-	public abstract fun getIdentity ()Ljava/lang/String;
-}
-
-public abstract interface class space/kscience/kmath/expressions/SymbolIndexer {
-	public abstract fun get (Ljava/util/List;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public abstract fun get (Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/expressions/Symbol;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public abstract fun get (Lspace/kscience/kmath/structures/Buffer;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public abstract fun get ([DLspace/kscience/kmath/expressions/Symbol;)D
-	public abstract fun get ([Ljava/lang/Object;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
-	public abstract fun getSymbols ()Ljava/util/List;
-	public abstract fun indexOf (Lspace/kscience/kmath/expressions/Symbol;)I
-	public abstract fun toDoubleArray (Ljava/util/Map;)[D
-	public abstract fun toList (Ljava/util/Map;)Ljava/util/List;
-	public abstract fun toMap ([D)Ljava/util/Map;
-	public abstract fun toPoint (Ljava/util/Map;Lkotlin/jvm/functions/Function2;)Lspace/kscience/kmath/structures/Buffer;
-}
-
 public final class space/kscience/kmath/expressions/SymbolIndexer$DefaultImpls {
-	public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Ljava/util/List;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
+	public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Ljava/util/List;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object;
-	public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/expressions/Symbol;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
+	public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/nd/Structure2D;Lspace/kscience/kmath/misc/Symbol;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object;
-	public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/structures/Buffer;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
+	public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/structures/Buffer;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object;
-	public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;[DLspace/kscience/kmath/expressions/Symbol;)D
+	public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;[DLspace/kscience/kmath/misc/Symbol;)D
-	public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;[Ljava/lang/Object;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
+	public static fun get (Lspace/kscience/kmath/expressions/SymbolIndexer;[Ljava/lang/Object;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object;
-	public static fun indexOf (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/expressions/Symbol;)I
+	public static fun indexOf (Lspace/kscience/kmath/expressions/SymbolIndexer;Lspace/kscience/kmath/misc/Symbol;)I
 	public static fun toDoubleArray (Lspace/kscience/kmath/expressions/SymbolIndexer;Ljava/util/Map;)[D
 	public static fun toList (Lspace/kscience/kmath/expressions/SymbolIndexer;Ljava/util/Map;)Ljava/util/List;
 	public static fun toMap (Lspace/kscience/kmath/expressions/SymbolIndexer;[D)Ljava/util/Map;
@@ -428,8 +376,6 @@ public final class space/kscience/kmath/expressions/SymbolIndexer$DefaultImpls {
 }
 
 public final class space/kscience/kmath/expressions/SymbolIndexerKt {
-	public static final fun withSymbols (Ljava/util/Collection;Lkotlin/jvm/functions/Function1;)Ljava/lang/Object;
-	public static final fun withSymbols ([Lspace/kscience/kmath/expressions/Symbol;Lkotlin/jvm/functions/Function1;)Ljava/lang/Object;
 }
 
 public final class space/kscience/kmath/linear/BufferedLinearSpace : space/kscience/kmath/linear/LinearSpace {
@@ -575,7 +521,7 @@ public final class space/kscience/kmath/linear/MatrixBuilderKt {
 	public static final fun row (Lspace/kscience/kmath/linear/LinearSpace;[Ljava/lang/Object;)Lspace/kscience/kmath/nd/Structure2D;
 }
 
-public abstract interface class space/kscience/kmath/linear/MatrixFeature {
+public abstract interface class space/kscience/kmath/linear/MatrixFeature : space/kscience/kmath/nd/StructureFeature {
 }
 
 public final class space/kscience/kmath/linear/MatrixFeaturesKt {
@@ -672,6 +618,35 @@ public final class space/kscience/kmath/misc/CumulativeKt {
 	public static final fun cumulativeSumOfLong (Lkotlin/sequences/Sequence;)Lkotlin/sequences/Sequence;
 }
 
+public final class space/kscience/kmath/misc/StringSymbol : space/kscience/kmath/misc/Symbol {
+	public static final synthetic fun box-impl (Ljava/lang/String;)Lspace/kscience/kmath/misc/StringSymbol;
+	public static fun constructor-impl (Ljava/lang/String;)Ljava/lang/String;
+	public fun equals (Ljava/lang/Object;)Z
+	public static fun equals-impl (Ljava/lang/String;Ljava/lang/Object;)Z
+	public static final fun equals-impl0 (Ljava/lang/String;Ljava/lang/String;)Z
+	public fun getIdentity ()Ljava/lang/String;
+	public fun hashCode ()I
+	public static fun hashCode-impl (Ljava/lang/String;)I
+	public fun toString ()Ljava/lang/String;
+	public static fun toString-impl (Ljava/lang/String;)Ljava/lang/String;
+	public final synthetic fun unbox-impl ()Ljava/lang/String;
+}
+
+public abstract interface class space/kscience/kmath/misc/Symbol {
+	public static final field Companion Lspace/kscience/kmath/misc/Symbol$Companion;
+	public abstract fun getIdentity ()Ljava/lang/String;
+}
+
+public final class space/kscience/kmath/misc/Symbol$Companion {
+	public final fun getX-tWtZOCg ()Ljava/lang/String;
+	public final fun getY-tWtZOCg ()Ljava/lang/String;
+	public final fun getZ-tWtZOCg ()Ljava/lang/String;
+}
+
+public final class space/kscience/kmath/misc/SymbolKt {
+	public static final fun getSymbol ()Lkotlin/properties/ReadOnlyProperty;
+}
+
 public abstract interface annotation class space/kscience/kmath/misc/UnstableKMathAPI : java/lang/annotation/Annotation {
 }
 
@@ -1086,11 +1061,15 @@ public final class space/kscience/kmath/nd/Strides$DefaultImpls {
 }
 
 public abstract interface class space/kscience/kmath/nd/Structure1D : space/kscience/kmath/nd/StructureND, space/kscience/kmath/structures/Buffer {
+	public static final field Companion Lspace/kscience/kmath/nd/Structure1D$Companion;
 	public abstract fun get ([I)Ljava/lang/Object;
 	public abstract fun getDimension ()I
 	public abstract fun iterator ()Ljava/util/Iterator;
 }
 
+public final class space/kscience/kmath/nd/Structure1D$Companion {
+}
+
 public final class space/kscience/kmath/nd/Structure1D$DefaultImpls {
 	public static fun get (Lspace/kscience/kmath/nd/Structure1D;[I)Ljava/lang/Object;
 	public static fun getDimension (Lspace/kscience/kmath/nd/Structure1D;)I
@@ -1132,6 +1111,9 @@ public final class space/kscience/kmath/nd/Structure2DKt {
 	public static final fun as2D (Lspace/kscience/kmath/nd/StructureND;)Lspace/kscience/kmath/nd/Structure2D;
 }
 
+public abstract interface class space/kscience/kmath/nd/StructureFeature {
+}
+
 public abstract interface class space/kscience/kmath/nd/StructureND {
 	public static final field Companion Lspace/kscience/kmath/nd/StructureND$Companion;
 	public abstract fun elements ()Lkotlin/sequences/Sequence;
@@ -1201,7 +1183,7 @@ public final class space/kscience/kmath/operations/AlgebraExtensionsKt {
 }
 
 public final class space/kscience/kmath/operations/AlgebraKt {
-	public static final fun bindSymbol (Lspace/kscience/kmath/operations/Algebra;Lspace/kscience/kmath/expressions/Symbol;)Ljava/lang/Object;
+	public static final fun bindSymbol (Lspace/kscience/kmath/operations/Algebra;Lspace/kscience/kmath/misc/Symbol;)Ljava/lang/Object;
 	public static final fun invoke (Lspace/kscience/kmath/operations/Algebra;Lkotlin/jvm/functions/Function1;)Ljava/lang/Object;
 }
 
@@ -2140,6 +2122,7 @@ public final class space/kscience/kmath/structures/BufferKt {
 public final class space/kscience/kmath/structures/BufferOperationKt {
 	public static final fun asIterable (Lspace/kscience/kmath/structures/Buffer;)Ljava/lang/Iterable;
 	public static final fun asSequence (Lspace/kscience/kmath/structures/Buffer;)Lkotlin/sequences/Sequence;
+	public static final fun fold (Lspace/kscience/kmath/structures/Buffer;Ljava/lang/Object;Lkotlin/jvm/functions/Function2;)Ljava/lang/Object;
 	public static final fun toList (Lspace/kscience/kmath/structures/Buffer;)Ljava/util/List;
 }
 

kmath-core/docs/README-TEMPLATE.md

@@ -1,6 +1,6 @@
-# The Core Module (`kmath-core`)
+# Module kmath-core
 
-The core features of KMath:
+The core interfaces of KMath.
 
 ${features}
 

kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/ColumnarData.kt

@@ -0,0 +1,34 @@
+package space.kscience.kmath.data
+
+import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.nd.Structure2D
+import space.kscience.kmath.structures.Buffer
+
+/**
+ * A column-based data set with all columns of the same size (not necessary fixed in time).
+ * The column could be retrieved by a [get] operation.
+ */
+@UnstableKMathAPI
+public interface ColumnarData<out T> {
+    public val size: Int
+
+    public operator fun get(symbol: Symbol): Buffer<T>
+}
+
+/**
+ * A zero-copy method to represent a [Structure2D] as a two-column x-y data.
+ * There could more than two columns in the structure.
+ */
+@UnstableKMathAPI
+public fun <T> Structure2D<T>.asColumnarData(mapping: Map<Symbol, Int>): ColumnarData<T> {
+    require(shape[1] >= mapping.maxOf { it.value }) { "Column index out of bounds" }
+    return object : ColumnarData<T> {
+        override val size: Int get() = shape[0]
+        override fun get(symbol: Symbol): Buffer<T> {
+            val index = mapping[symbol] ?: error("No column mapping for symbol $symbol")
+            return columns[index]
+        }
+    }
+}
+

kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYColumnarData.kt

@@ -0,0 +1,55 @@
+package space.kscience.kmath.data
+
+import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.nd.Structure2D
+import space.kscience.kmath.structures.Buffer
+import kotlin.math.max
+
+/**
+ * The buffer of X values.
+ */
+@UnstableKMathAPI
+public interface XYColumnarData<T, out X : T, out Y : T> : ColumnarData<T> {
+    /**
+     * The buffer of X values
+     */
+    public val x: Buffer<X>
+
+    /**
+     * The buffer of Y values.
+     */
+    public val y: Buffer<Y>
+
+    override fun get(symbol: Symbol): Buffer<T> = when (symbol) {
+        Symbol.x -> x
+        Symbol.y -> y
+        else -> error("A column for symbol $symbol not found")
+    }
+}
+
+@Suppress("FunctionName")
+@UnstableKMathAPI
+public fun <T, X : T, Y : T> XYColumnarData(x: Buffer<X>, y: Buffer<Y>): XYColumnarData<T, X, Y> {
+    require(x.size == y.size) { "Buffer size mismatch. x buffer size is ${x.size}, y buffer size is ${y.size}" }
+    return object : XYColumnarData<T, X, Y> {
+        override val size: Int = x.size
+        override val x: Buffer<X> = x
+        override val y: Buffer<Y> = y
+    }
+}
+
+
+/**
+ * A zero-copy method to represent a [Structure2D] as a two-column x-y data.
+ * There could more than two columns in the structure.
+ */
+@UnstableKMathAPI
+public fun <T> Structure2D<T>.asXYData(xIndex: Int = 0, yIndex: Int = 1): XYColumnarData<T, T, T> {
+    require(shape[1] >= max(xIndex, yIndex)) { "Column index out of bounds" }
+    return object : XYColumnarData<T, T, T> {
+        override val size: Int get() = this@asXYData.shape[0]
+        override val x: Buffer<T> get() = columns[xIndex]
+        override val y: Buffer<T> get() = columns[yIndex]
+    }
+}

kmath-core/src/commonMain/kotlin/space/kscience/kmath/data/XYZColumnarData.kt

@@ -0,0 +1,21 @@
+package space.kscience.kmath.data
+
+import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.structures.Buffer
+
+/**
+ * A [XYColumnarData] with guaranteed [x], [y] and [z] columns designated by corresponding symbols.
+ * Inherits [XYColumnarData].
+ */
+@UnstableKMathAPI
+public interface XYZColumnarData<T, out X : T, out Y : T, out Z : T> : XYColumnarData<T, X, Y> {
+    public val z: Buffer<Z>
+
+    override fun get(symbol: Symbol): Buffer<T> = when (symbol) {
+        Symbol.x -> x
+        Symbol.y -> y
+        Symbol.z -> z
+        else -> error("A column for symbol $symbol not found")
+    }
+}

kmath-core/src/commonMain/kotlin/space/kscience/kmath/domains/DoubleDomain.kt

@@ -24,7 +24,6 @@ import space.kscience.kmath.misc.UnstableKMathAPI
  */
 @UnstableKMathAPI
 public interface DoubleDomain : Domain<Double> {
-
     /**
      * Global lower edge
      * @param num axis number

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

@@ -1,5 +1,8 @@
 package space.kscience.kmath.expressions
 
+import space.kscience.kmath.misc.StringSymbol
+import space.kscience.kmath.misc.Symbol
+
 /**
  * Represents expression which structure can be differentiated.
  *

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

@@ -1,26 +1,11 @@
 package space.kscience.kmath.expressions
 
+import space.kscience.kmath.misc.StringSymbol
+import space.kscience.kmath.misc.Symbol
 import space.kscience.kmath.operations.Algebra
 import kotlin.jvm.JvmName
 import kotlin.properties.ReadOnlyProperty
 
-/**
- * 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.
  *
@@ -92,13 +77,6 @@ public interface ExpressionAlgebra<in T, E> : Algebra<E> {
 public fun <T, E> ExpressionAlgebra<T, E>.bindSymbol(symbol: Symbol): E =
     bindSymbolOrNull(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?, Symbol> = ReadOnlyProperty { _, property ->
-    StringSymbol(property.name)
-}
-
 /**
  * Bind a symbol by name inside the [ExpressionAlgebra]
  */

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

@@ -1,5 +1,6 @@
 package space.kscience.kmath.expressions
 
+import space.kscience.kmath.misc.Symbol
 import space.kscience.kmath.operations.*
 
 /**

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

@@ -1,6 +1,7 @@
 package space.kscience.kmath.expressions
 
 import space.kscience.kmath.linear.Point
+import space.kscience.kmath.misc.Symbol
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.*
 import space.kscience.kmath.structures.asBuffer

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

@@ -1,6 +1,8 @@
 package space.kscience.kmath.expressions
 
 import space.kscience.kmath.linear.Point
+import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.nd.Structure2D
 import space.kscience.kmath.structures.BufferFactory
 
@@ -8,6 +10,7 @@ import space.kscience.kmath.structures.BufferFactory
  * An environment to easy transform indexed variables to symbols and back.
  * TODO requires multi-receivers to be beautiful
  */
+@UnstableKMathAPI
 public interface SymbolIndexer {
     public val symbols: List<Symbol>
     public fun indexOf(symbol: Symbol): Int = symbols.indexOf(symbol)
@@ -49,13 +52,16 @@ public interface SymbolIndexer {
     public fun Map<Symbol, Double>.toDoubleArray(): DoubleArray = DoubleArray(symbols.size) { getValue(symbols[it]) }
 }
 
+@UnstableKMathAPI
 public inline class SimpleSymbolIndexer(override val symbols: List<Symbol>) : SymbolIndexer
 
 /**
  * Execute the block with symbol indexer based on given symbol order
  */
+@UnstableKMathAPI
 public inline fun <R> withSymbols(vararg symbols: Symbol, block: SymbolIndexer.() -> R): R =
     with(SimpleSymbolIndexer(symbols.toList()), block)
 
+@UnstableKMathAPI
 public inline fun <R> withSymbols(symbols: Collection<Symbol>, block: SymbolIndexer.() -> R): R =
     with(SimpleSymbolIndexer(symbols.toList()), block)

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSolver.kt

@@ -3,7 +3,10 @@ package space.kscience.kmath.linear
 import space.kscience.kmath.nd.as1D
 
 /**
- * A group of methods to resolve equation A dot X = B, where A and B are matrices or vectors
+ * A group of methods to solve for *X* in equation *X = A <sup>-1</sup> &middot; B*, where *A* and *B* are matrices or
+ * vectors.
+ *
+ * @param T the type of items.
  */
 public interface LinearSolver<T : Any> {
     /**
@@ -23,7 +26,7 @@ public interface LinearSolver<T : Any> {
 }
 
 /**
- * Convert matrix to vector if it is possible
+ * Convert matrix to vector if it is possible.
  */
 public fun <T : Any> Matrix<T>.asVector(): Point<T> =
     if (this.colNum == 1)
@@ -31,4 +34,11 @@ public fun <T : Any> Matrix<T>.asVector(): Point<T> =
     else
         error("Can't convert matrix with more than one column to vector")
 
+/**
+ * Creates an n &times; 1 [VirtualMatrix], where n is the size of the given buffer.
+ *
+ * @param T the type of elements contained in the buffer.
+ * @receiver a buffer.
+ * @return the new matrix.
+ */
 public fun <T : Any> Point<T>.asMatrix(): VirtualMatrix<T> = VirtualMatrix(size, 1) { i, _ -> get(i) }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/LinearSpace.kt

@@ -18,6 +18,8 @@ public typealias MutableMatrix<T> = MutableStructure2D<T>
 
 /**
  * Alias or using [Buffer] as a point/vector in a many-dimensional space.
+ *
+ * @param T the type of elements contained in the buffer.
  */
 public typealias Point<T> = Buffer<T>
 
@@ -165,7 +167,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
      * @return a feature object or `null` if it isn't present.
      */
     @UnstableKMathAPI
-    public fun <F : Any> getFeature(structure: Matrix<T>, type: KClass<F>): F? = structure.getFeature(type)
+    public fun <F : StructureFeature> getFeature(structure: Matrix<T>, type: KClass<out F>): F? = structure.getFeature(type)
 
     public companion object {
 
@@ -195,7 +197,7 @@ public interface LinearSpace<T : Any, out A : Ring<T>> {
  * @return a feature object or `null` if it isn't present.
  */
 @UnstableKMathAPI
-public inline fun <T : Any, reified F : Any> LinearSpace<T, *>.getFeature(structure: Matrix<T>): F? =
+public inline fun <T : Any, reified F : StructureFeature> LinearSpace<T, *>.getFeature(structure: Matrix<T>): F? =
     getFeature(structure, F::class)
 
 

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixFeatures.kt

@@ -1,10 +1,12 @@
 package space.kscience.kmath.linear
 
+import space.kscience.kmath.nd.StructureFeature
+
 /**
  * A marker interface representing some properties of matrices or additional transformations of them. Features are used
  * to optimize matrix operations performance in some cases or retrieve the APIs.
  */
-public interface MatrixFeature
+public interface MatrixFeature: StructureFeature
 
 /**
  * Matrices with this feature are considered to have only diagonal non-null elements.

kmath-core/src/commonMain/kotlin/space/kscience/kmath/linear/MatrixWrapper.kt

@@ -1,6 +1,7 @@
 package space.kscience.kmath.linear
 
 import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.nd.StructureFeature
 import space.kscience.kmath.nd.getFeature
 import space.kscience.kmath.operations.Ring
 import kotlin.reflect.KClass
@@ -20,7 +21,7 @@ public class MatrixWrapper<T : Any> internal constructor(
      */
     @UnstableKMathAPI
     @Suppress("UNCHECKED_CAST")
-    override fun <T : Any> getFeature(type: KClass<T>): T? = features.singleOrNull { type.isInstance(it) } as? T
+    override fun <F : StructureFeature> getFeature(type: KClass<out F>): F? = features.singleOrNull { type.isInstance(it) } as? F
         ?: origin.getFeature(type)
 
     override fun toString(): String {

kmath-core/src/commonMain/kotlin/space/kscience/kmath/misc/Symbol.kt

@@ -0,0 +1,34 @@
+package space.kscience.kmath.misc
+
+import kotlin.properties.ReadOnlyProperty
+
+/**
+ * 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
+
+    public companion object{
+        public val x: StringSymbol = StringSymbol("x")
+        public val y: StringSymbol = StringSymbol("y")
+        public val z: StringSymbol = StringSymbol("z")
+    }
+}
+
+/**
+ * A [Symbol] with a [String] identity
+ */
+public inline class StringSymbol(override val identity: String) : Symbol {
+    override fun toString(): String = identity
+}
+
+
+/**
+ * A delegate to create a symbol with a string identity in this scope
+ */
+public val symbol: ReadOnlyProperty<Any?, Symbol> = ReadOnlyProperty { _, property ->
+    StringSymbol(property.name)
+}

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt

@@ -67,7 +67,8 @@ public interface AlgebraND<T, C : Algebra<T>> {
      * @return a feature object or `null` if it isn't present.
      */
     @UnstableKMathAPI
-    public fun <F : Any> getFeature(structure: StructureND<T>, type: KClass<F>): F? = structure.getFeature(type)
+    public fun <F : StructureFeature> getFeature(structure: StructureND<T>, type: KClass<out F>): F? =
+        structure.getFeature(type)
 
     public companion object
 }
@@ -81,7 +82,7 @@ public interface AlgebraND<T, C : Algebra<T>> {
  * @return a feature object or `null` if it isn't present.
  */
 @UnstableKMathAPI
-public inline fun <T : Any, reified F : Any> AlgebraND<T, *>.getFeature(structure: StructureND<T>): F? =
+public inline fun <T : Any, reified F : StructureFeature> AlgebraND<T, *>.getFeature(structure: StructureND<T>): F? =
     getFeature(structure, F::class)
 
 /**

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure1D.kt

@@ -17,6 +17,8 @@ public interface Structure1D<T> : StructureND<T>, Buffer<T> {
     }
 
     public override operator fun iterator(): Iterator<T> = (0 until size).asSequence().map(::get).iterator()
+
+    public companion object
 }
 
 /**

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/Structure2D.kt

@@ -97,7 +97,7 @@ private inline class Structure2DWrapper<T>(val structure: StructureND<T>) : Stru
     override operator fun get(i: Int, j: Int): T = structure[i, j]
 
     @UnstableKMathAPI
-    override fun <F : Any> getFeature(type: KClass<F>): F? = structure.getFeature(type)
+    override fun <F : StructureFeature> getFeature(type: KClass<out F>): F? = structure.getFeature(type)
 
     override fun elements(): Sequence<Pair<IntArray, T>> = structure.elements()
 }

kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/StructureND.kt

@@ -7,6 +7,8 @@ import kotlin.jvm.JvmName
 import kotlin.native.concurrent.ThreadLocal
 import kotlin.reflect.KClass
 
+public interface StructureFeature
+
 /**
  * Represents n-dimensional structure, i.e. multidimensional container of items of the same type and size. The number
  * of dimensions and items in an array is defined by its shape, which is a sequence of non-negative integers that
@@ -48,7 +50,7 @@ public interface StructureND<T> {
      * If the feature is not present, null is returned.
      */
     @UnstableKMathAPI
-    public fun <F : Any> getFeature(type: KClass<F>): F? = null
+    public fun <F : StructureFeature> getFeature(type: KClass<out F>): F? = null
 
     public companion object {
         /**
@@ -144,7 +146,7 @@ public interface StructureND<T> {
 public operator fun <T> StructureND<T>.get(vararg index: Int): T = get(index)
 
 @UnstableKMathAPI
-public inline fun <reified T : Any> StructureND<*>.getFeature(): T? = getFeature(T::class)
+public inline fun <reified T : StructureFeature> StructureND<*>.getFeature(): T? = getFeature(T::class)
 
 /**
  * Represents mutable [StructureND].

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/Algebra.kt

@@ -1,6 +1,6 @@
 package space.kscience.kmath.operations
 
-import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.misc.Symbol
 
 /**
  * Stub for DSL the [Algebra] is.
@@ -100,8 +100,8 @@ public fun <T : Any> Algebra<T>.bindSymbol(symbol: Symbol): T = bindSymbol(symbo
 public inline operator fun <A : Algebra<*>, R> A.invoke(block: A.() -> R): R = run(block)
 
 /**
- * Represents linear space without neutral element, i.e. algebraic structure with associative, binary operation [add]
+ * Represents group without neutral element (also known as inverse semigroup), i.e. algebraic structure with
- * and scalar multiplication [multiply].
+ * associative, binary operation [add].
  *
  * @param T the type of element of this semispace.
  */
@@ -177,7 +177,7 @@ public interface GroupOperations<T> : Algebra<T> {
 }
 
 /**
- * Represents linear space with neutral element, i.e. algebraic structure with associative, binary operation [add].
+ * Represents group, i.e. algebraic structure with associative, binary operation [add].
  *
  * @param T the type of element of this semispace.
  */
@@ -189,8 +189,8 @@ public interface Group<T> : GroupOperations<T> {
 }
 
 /**
- * Represents rng, i.e. algebraic structure with associative, binary, commutative operation [add] and associative,
+ * Represents ring without multiplicative and additive identities, i.e. algebraic structure with
- * operation [multiply] distributive over [add].
+ * associative, binary, commutative operation [add] and associative, operation [multiply] distributive over [add].
  *
  * @param T the type of element of this semiring.
  */
@@ -238,7 +238,7 @@ public interface Ring<T> : Group<T>, RingOperations<T> {
 }
 
 /**
- * Represents field without identity elements, i.e. algebraic structure with associative, binary, commutative operations
+ * Represents field without without multiplicative and additive identities, i.e. algebraic structure with associative, binary, commutative operations
  * [add] and [multiply]; binary operation [divide] as multiplication of left operand by reciprocal of right one.
  *
  * @param T the type of element of this semifield.
@@ -276,10 +276,11 @@ public interface FieldOperations<T> : RingOperations<T> {
 }
 
 /**
- * Represents field, i.e. algebraic structure with three operations: associative "addition" and "multiplication",
+ * Represents field, i.e. algebraic structure with three operations: associative, commutative addition and
- * and "division" and their neutral elements.
+ * multiplication, and division. **This interface differs from the eponymous mathematical definition: fields in KMath
+ * also support associative multiplication by scalar.**
  *
- * @param T the type of element of this semifield.
+ * @param T the type of element of this field.
  */
 public interface Field<T> : Ring<T>, FieldOperations<T>, ScaleOperations<T>, NumericAlgebra<T> {
     override fun number(value: Number): T = scale(one, value.toDouble())

kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BigInt.kt

@@ -12,8 +12,8 @@ import kotlin.math.max
 import kotlin.math.min
 import kotlin.math.sign
 
-public typealias Magnitude = UIntArray
+private typealias Magnitude = UIntArray
-public typealias TBase = ULong
+private typealias TBase = ULong
 
 /**
  * Kotlin Multiplatform implementation of Big Integer numbers (KBigInteger).
@@ -241,18 +241,18 @@ public class BigInt internal constructor(
         )
 
         private fun compareMagnitudes(mag1: Magnitude, mag2: Magnitude): Int {
-            when {
+            return when {
-                mag1.size > mag2.size -> return 1
+                mag1.size > mag2.size -> 1
-                mag1.size < mag2.size -> return -1
+                mag1.size < mag2.size -> -1
+
                 else -> {
-                    for (i in mag1.size - 1 downTo 0) {
+                    for (i in mag1.size - 1 downTo 0) return when {
-                        if (mag1[i] > mag2[i]) {
+                        mag1[i] > mag2[i] -> 1
-                            return 1
+                        mag1[i] < mag2[i] -> -1
-                        } else if (mag1[i] < mag2[i]) {
+                        else -> continue
-                            return -1
                     }
-                    }
+
-                    return 0
+                    0
                 }
             }
         }
@@ -302,10 +302,11 @@ public class BigInt internal constructor(
             var carry = 0uL
 
             for (i in mag.indices) {
-                val cur: ULong = carry + mag[i].toULong() * x.toULong()
+                val cur = carry + mag[i].toULong() * x.toULong()
                 result[i] = (cur and BASE).toUInt()
                 carry = cur shr BASE_SIZE
             }
+
             result[resultLength - 1] = (carry and BASE).toUInt()
 
             return stripLeadingZeros(result)
@@ -358,6 +359,9 @@ private fun stripLeadingZeros(mag: Magnitude): Magnitude {
     return mag.sliceArray(IntRange(0, resSize))
 }
 
+/**
+ * Returns the absolute value of the given value [x].
+ */
 public fun abs(x: BigInt): BigInt = x.abs()
 
 /**

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/Buffer.kt

@@ -40,7 +40,6 @@ public interface Buffer<out T> {
     public operator fun iterator(): Iterator<T>
 
     public companion object {
-
         /**
          * Check the element-by-element match of content of two buffers.
          */
@@ -110,7 +109,6 @@ public interface MutableBuffer<T> : Buffer<T> {
     public fun copy(): MutableBuffer<T>
 
     public companion object {
-
         /**
          * Creates a [DoubleBuffer] with the specified [size], where each element is calculated by calling the specified
          * [initializer] function.

kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferOperation.kt

@@ -70,6 +70,15 @@ public inline fun <T : Any, reified R : Any> Buffer<T>.mapIndexed(
     crossinline block: (index: Int, value: T) -> R,
 ): Buffer<R> = bufferFactory(size) { block(it, get(it)) }
 
+/**
+ * Fold given buffer according to [operation]
+ */
+public inline fun <T : Any, R> Buffer<T>.fold(initial: R, operation: (acc: R, T) -> R): R {
+    var accumulator = initial
+    for (index in this.indices) accumulator = operation(accumulator, get(index))
+    return accumulator
+}
+
 /**
  * Zip two buffers using given [transform].
  */

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

@@ -1,5 +1,6 @@
 package space.kscience.kmath.expressions
 
+import space.kscience.kmath.misc.symbol
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.invoke
 import kotlin.test.Test

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

@@ -1,5 +1,7 @@
 package space.kscience.kmath.expressions
 
+import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.misc.symbol
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.asBuffer

kmath-core/src/commonTest/kotlin/space/kscience/kmath/structures/NumberNDFieldTest.kt

@@ -38,13 +38,12 @@ class NumberNDFieldTest {
             (i * 10 + j).toDouble()
         }
 
-        for (i in 0..2) {
+        for (i in 0..2)
             for (j in 0..2) {
                 val expected = (i * 10 + j).toDouble()
                 assertEquals(expected, array[i, j], "Error at index [$i, $j]")
             }
     }
-    }
 
     @Test
     fun testExternalFunction() {

kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingDoubleChain.kt

@@ -1,12 +1,13 @@
 package space.kscience.kmath.chains
 
 /**
- * Performance optimized chain for real values
+ * Chunked, specialized chain for real values.
  */
-public abstract class BlockingDoubleChain : Chain<Double> {
+public interface BlockingDoubleChain : Chain<Double> {
-    public abstract fun nextDouble(): Double
+    public override suspend fun next(): Double
 
-    override suspend fun next(): Double = nextDouble()
+    /**
-
+     * Returns an [DoubleArray] chunk of [size] values of [next].
-    public fun nextBlock(size: Int): DoubleArray = DoubleArray(size) { nextDouble() }
+     */
+    public suspend fun nextBlock(size: Int): DoubleArray = DoubleArray(size) { next() }
 }

kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/BlockingIntChain.kt

@@ -3,10 +3,7 @@ package space.kscience.kmath.chains
 /**
  * Performance optimized chain for integer values
  */
-public abstract class BlockingIntChain : Chain<Int> {
+public interface BlockingIntChain : Chain<Int> {
-    public abstract fun nextInt(): Int
+    public override suspend fun next(): Int
-
+    public suspend fun nextBlock(size: Int): IntArray = IntArray(size) { next() }
-    override suspend fun next(): Int = nextInt()
-
-    public fun nextBlock(size: Int): IntArray = IntArray(size) { nextInt() }
 }

kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/chains/Chain.kt

@@ -63,12 +63,10 @@ public class MarkovChain<out R : Any>(private val seed: suspend () -> R, private
 
     public fun value(): R? = value
 
-    public override suspend fun next(): R {
+    public override suspend fun next(): R = mutex.withLock {
-        mutex.withLock {
         val newValue = gen(value ?: seed())
         value = newValue
-            return newValue
+        newValue
-        }
     }
 
     public override fun fork(): Chain<R> = MarkovChain(seed = { value ?: seed() }, gen = gen)
@@ -90,12 +88,10 @@ public class StatefulChain<S, out R>(
 
     public fun value(): R? = value
 
-    public override suspend fun next(): R {
+    public override suspend fun next(): R = mutex.withLock {
-        mutex.withLock {
         val newValue = state.gen(value ?: state.seed())
         value = newValue
-            return newValue
+        newValue
-        }
     }
 
     public override fun fork(): Chain<R> = StatefulChain(forkState(state), seed, forkState, gen)

kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/BufferFlow.kt

@@ -28,7 +28,7 @@ public fun <T> Flow<T>.chunked(bufferSize: Int, bufferFactory: BufferFactory<T>)
     var counter = 0
 
     this@chunked.collect { element ->
-        list.add(element)
+        list += element
         counter++
 
         if (counter == bufferSize) {

kmath-coroutines/src/commonMain/kotlin/space/kscience/kmath/streaming/RingBuffer.kt

@@ -48,11 +48,9 @@ public class RingBuffer<T>(
     /**
      * A safe snapshot operation
      */
-    public suspend fun snapshot(): Buffer<T> {
+    public suspend fun snapshot(): Buffer<T> = mutex.withLock {
-        mutex.withLock {
         val copy = buffer.copy()
-            return VirtualBuffer(size) { i -> copy[startIndex.forward(i)] as T }
+        VirtualBuffer(size) { i -> copy[startIndex.forward(i)] as T }
-        }
     }
 
     public suspend fun push(element: T) {

kmath-dimensions/src/commonTest/kotlin/kscience/dimensions/DMatrixContextTest.kt

@@ -1,4 +1,4 @@
-package kscience.dimensions
+package space.kscience.dimensions
 
 import space.kscience.kmath.dimensions.D2
 import space.kscience.kmath.dimensions.D3

kmath-ejml/README.md

@@ -1,43 +1,37 @@
-# ejml-simple support (`kmath-ejml`)
+# Module kmath-ejml
 
-This subproject implements the following features:
+EJML based linear algebra implementation.
 
  - [ejml-vector](src/main/kotlin/space/kscience/kmath/ejml/EjmlVector.kt) : The Point implementation using SimpleMatrix.
  - [ejml-matrix](src/main/kotlin/space/kscience/kmath/ejml/EjmlMatrix.kt) : The Matrix implementation using SimpleMatrix.
  - [ejml-linear-space](src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt) : The LinearSpace implementation using SimpleMatrix.
 
 
-> #### Artifact:
+## Artifact:
->
+
-> This module artifact: `space.kscience:kmath-ejml:0.3.0-dev-3`.
+The Maven coordinates of this project are `space.kscience:kmath-ejml:0.3.0-dev-3`.
->
+
-> Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-ejml/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-ejml/_latestVersion)
+**Gradle:**
->
+```gradle
-> Bintray development version:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-ejml/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-ejml/_latestVersion)
+repositories {
->
+    maven { url 'https://repo.kotlin.link' }
-> **Gradle:**
+    maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
->
+    maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
-> ```gradle
+}
-> repositories {
+
->     maven { url 'https://repo.kotlin.link' }
+dependencies {
->     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+    implementation 'space.kscience:kmath-ejml:0.3.0-dev-3'
->     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
+}
-> }
+```
-> 
+**Gradle Kotlin DSL:**
-> dependencies {
+```kotlin
->     implementation 'space.kscience:kmath-ejml:0.3.0-dev-3'
+repositories {
-> }
+    maven("https://repo.kotlin.link")
-> ```
+    maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
-> **Gradle Kotlin DSL:**
+    maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
->
+}
-> ```kotlin
+
-> repositories {
+dependencies {
->     maven("https://repo.kotlin.link")
+    implementation("space.kscience:kmath-ejml:0.3.0-dev-3")
->     maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
+}
->     maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
+```
-> }
-> 
-> dependencies {
->     implementation("space.kscience:kmath-ejml:0.3.0-dev-3")
-> }
-> ```

kmath-ejml/docs/README-TEMPLATE.md

@@ -1,6 +1,6 @@
-# ejml-simple support (`kmath-ejml`)
+# Module kmath-ejml
 
-This subproject implements the following features:
+EJML based linear algebra implementation.
 
 ${features}
 

kmath-ejml/src/main/kotlin/space/kscience/kmath/ejml/EjmlLinearSpace.kt

@@ -4,6 +4,7 @@ import org.ejml.dense.row.factory.DecompositionFactory_DDRM
 import org.ejml.simple.SimpleMatrix
 import space.kscience.kmath.linear.*
 import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.nd.StructureFeature
 import space.kscience.kmath.nd.getFeature
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.structures.DoubleBuffer
@@ -95,7 +96,7 @@ public object EjmlLinearSpace : LinearSpace<Double, DoubleField> {
         v.toEjml().origin.scale(this).wrapVector()
 
     @UnstableKMathAPI
-    public override fun <F : Any> getFeature(structure: Matrix<Double>, type: KClass<F>): F? {
+    public override fun <F : StructureFeature> getFeature(structure: Matrix<Double>, type: KClass<out F>): F? {
         //Return the feature if it is intrinsic to the structure
         structure.getFeature(type)?.let { return it }
 

kmath-for-real/README.md

@@ -1,41 +1,37 @@
-# Real number specialization module (`kmath-for-real`)
+# Module kmath-for-real
+
+Specialization of KMath APIs for Double numbers.
 
  - [DoubleVector](src/commonMain/kotlin/space/kscience/kmath/real/DoubleVector.kt) : Numpy-like operations for Buffers/Points
  - [DoubleMatrix](src/commonMain/kotlin/space/kscience/kmath/real/DoubleMatrix.kt) : Numpy-like operations for 2d real structures
  - [grids](src/commonMain/kotlin/space/kscience/kmath/structures/grids.kt) : Uniform grid generators
 
 
-> #### Artifact:
+## Artifact:
->
+
-> This module artifact: `space.kscience:kmath-for-real:0.3.0-dev-3`.
+The Maven coordinates of this project are `space.kscience:kmath-for-real:0.3.0-dev-3`.
->
+
-> Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-for-real/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-for-real/_latestVersion)
+**Gradle:**
->
+```gradle
-> Bintray development version:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-for-real/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-for-real/_latestVersion)
+repositories {
->
+    maven { url 'https://repo.kotlin.link' }
-> **Gradle:**
+    maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
->
+    maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
-> ```gradle
+}
-> repositories {
+
->     maven { url 'https://repo.kotlin.link' }
+dependencies {
->     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+    implementation 'space.kscience:kmath-for-real:0.3.0-dev-3'
->     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
+}
-> }
+```
-> 
+**Gradle Kotlin DSL:**
-> dependencies {
+```kotlin
->     implementation 'space.kscience:kmath-for-real:0.3.0-dev-3'
+repositories {
-> }
+    maven("https://repo.kotlin.link")
-> ```
+    maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
-> **Gradle Kotlin DSL:**
+    maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
->
+}
-> ```kotlin
+
-> repositories {
+dependencies {
->     maven("https://repo.kotlin.link")
+    implementation("space.kscience:kmath-for-real:0.3.0-dev-3")
->     maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
+}
->     maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
+```
-> }
-> 
-> dependencies {
->     implementation("space.kscience:kmath-for-real:0.3.0-dev-3")
-> }
-> ```

kmath-for-real/docs/README-TEMPLATE.md

@@ -1,4 +1,6 @@
-# Real number specialization module (`kmath-for-real`)
+# Module kmath-for-real
+
+Specialization of KMath APIs for Double numbers.
 
 ${features}
 

kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/RealVector.kt

@@ -6,17 +6,13 @@ import space.kscience.kmath.operations.Norm
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.MutableBuffer.Companion.double
 import space.kscience.kmath.structures.asBuffer
-import space.kscience.kmath.structures.asIterable
+import space.kscience.kmath.structures.fold
 import space.kscience.kmath.structures.indices
 import kotlin.math.pow
 import kotlin.math.sqrt
 
 public typealias DoubleVector = Point<Double>
 
-public object VectorL2Norm : Norm<Point<out Number>, Double> {
-    override fun norm(arg: Point<out Number>): Double = sqrt(arg.asIterable().sumByDouble(Number::toDouble))
-}
-
 @Suppress("FunctionName")
 public fun DoubleVector(vararg doubles: Double): DoubleVector = doubles.asBuffer()
 
@@ -103,3 +99,9 @@ public fun DoubleVector.sum(): Double {
     }
     return res
 }
+
+public object VectorL2Norm : Norm<DoubleVector, Double> {
+    override fun norm(arg: DoubleVector): Double = sqrt(arg.fold(0.0) { acc: Double, d: Double -> acc + d.pow(2) })
+}
+
+public val DoubleVector.norm: Double get() = VectorL2Norm.norm(this)

kmath-for-real/src/commonMain/kotlin/space/kscience/kmath/real/grids.kt

@@ -1,7 +1,41 @@
 package space.kscience.kmath.real
 
-import space.kscience.kmath.structures.asBuffer
+import space.kscience.kmath.misc.UnstableKMathAPI
-import kotlin.math.abs
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.DoubleBuffer
+import kotlin.math.floor
+
+public val ClosedFloatingPointRange<Double>.length: Double get() = endInclusive - start
+
+/**
+ * Create a Buffer-based grid with equally distributed [numberOfPoints] points. The range could be increasing or  decreasing.
+ * If range has a zero size, then the buffer consisting of [numberOfPoints] equal values is returned.
+ */
+public fun Buffer.Companion.fromRange(range: ClosedFloatingPointRange<Double>, numberOfPoints: Int): DoubleBuffer {
+    require(numberOfPoints >= 2) { "Number of points in grid must be more than 1" }
+    val normalizedRange = when {
+        range.endInclusive > range.start -> range
+        range.endInclusive < range.start -> range.endInclusive..range.start
+        else -> return DoubleBuffer(numberOfPoints) { range.start }
+    }
+    val step = normalizedRange.length / (numberOfPoints - 1)
+    return DoubleBuffer(numberOfPoints) { normalizedRange.start + step * it }
+}
+
+/**
+ * Create a Buffer-based grid with equally distributed points with a fixed [step]. The range could be increasing or  decreasing.
+ * If the step is larger than the range size, single point is returned.
+ */
+public fun Buffer.Companion.withFixedStep(range: ClosedFloatingPointRange<Double>, step: Double): DoubleBuffer {
+    require(step > 0) { "The grid step must be positive" }
+    val normalizedRange = when {
+        range.endInclusive > range.start -> range
+        range.endInclusive < range.start -> range.endInclusive..range.start
+        else -> return DoubleBuffer(range.start)
+    }
+    val numberOfPoints = floor(normalizedRange.length / step).toInt() + 1
+    return DoubleBuffer(numberOfPoints) { normalizedRange.start + step * it  }
+}
 
 /**
  * Convert double range to sequence.
@@ -11,35 +45,5 @@ import kotlin.math.abs
  *
  * If step is negative, the same goes from upper boundary downwards
  */
-public fun ClosedFloatingPointRange<Double>.toSequenceWithStep(step: Double): Sequence<Double> = when {
+@UnstableKMathAPI
-    step == 0.0 -> error("Zero step in double progression")
+public infix fun ClosedFloatingPointRange<Double>.step(step: Double): DoubleBuffer = Buffer.withFixedStep(this, step)
-
-    step > 0 -> sequence {
-        var current = start
-
-        while (current <= endInclusive) {
-            yield(current)
-            current += step
-        }
-    }
-
-    else -> sequence {
-        var current = endInclusive
-
-        while (current >= start) {
-            yield(current)
-            current += step
-        }
-    }
-}
-
-public infix fun ClosedFloatingPointRange<Double>.step(step: Double): DoubleVector =
-    toSequenceWithStep(step).toList().asBuffer()
-
-/**
- * Convert double range to sequence with the fixed number of points
- */
-public fun ClosedFloatingPointRange<Double>.toSequenceWithPoints(numPoints: Int): Sequence<Double> {
-    require(numPoints > 1) { "The number of points should be more than 2" }
-    return toSequenceWithStep(abs(endInclusive - start) / (numPoints - 1))
-}

kmath-for-real/src/commonTest/kotlin/kaceince/kmath/real/GridTest.kt

@@ -1,13 +1,18 @@
 package kaceince.kmath.real
 
+import space.kscience.kmath.real.DoubleVector
+import space.kscience.kmath.real.minus
+import space.kscience.kmath.real.norm
 import space.kscience.kmath.real.step
 import kotlin.test.Test
 import kotlin.test.assertEquals
+import kotlin.test.assertTrue
 
 class GridTest {
     @Test
     fun testStepGrid() {
         val grid = 0.0..1.0 step 0.2
         assertEquals(6, grid.size)
+        assertTrue { (grid - DoubleVector(0.0, 0.2, 0.4, 0.6, 0.8, 1.0)).norm < 1e-4 }
     }
 }

kmath-functions/README.md

@@ -1,6 +1,6 @@
-# Functions (`kmath-functions`)
+# Module kmath-functions
 
-Functions and interpolations:
+Functions and interpolations.
 
  - [piecewise](Piecewise functions.) : src/commonMain/kotlin/space/kscience/kmath/functions/Piecewise.kt
  - [polynomials](Polynomial functions.) : src/commonMain/kotlin/space/kscience/kmath/functions/Polynomial.kt
@@ -8,37 +8,31 @@ Functions and interpolations:
  - [spline interpolation](Cubic spline XY interpolator.) : src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt
 
 
-> #### Artifact:
+## Artifact:
->
+
-> This module artifact: `space.kscience:kmath-functions:0.3.0-dev-3`.
+The Maven coordinates of this project are `space.kscience:kmath-functions:0.3.0-dev-3`.
->
+
-> Bintray release version:        [ ![Download](https://api.bintray.com/packages/mipt-npm/kscience/kmath-functions/images/download.svg) ](https://bintray.com/mipt-npm/kscience/kmath-functions/_latestVersion)
+**Gradle:**
->
+```gradle
-> Bintray development version:    [ ![Download](https://api.bintray.com/packages/mipt-npm/dev/kmath-functions/images/download.svg) ](https://bintray.com/mipt-npm/dev/kmath-functions/_latestVersion)
+repositories {
->
+    maven { url 'https://repo.kotlin.link' }
-> **Gradle:**
+    maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
->
+    maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
-> ```gradle
+}
-> repositories {
+
->     maven { url 'https://repo.kotlin.link' }
+dependencies {
->     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+    implementation 'space.kscience:kmath-functions:0.3.0-dev-3'
->     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
+}
-> }
+```
-> 
+**Gradle Kotlin DSL:**
-> dependencies {
+```kotlin
->     implementation 'space.kscience:kmath-functions:0.3.0-dev-3'
+repositories {
-> }
+    maven("https://repo.kotlin.link")
-> ```
+    maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
-> **Gradle Kotlin DSL:**
+    maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
->
+}
-> ```kotlin
+
-> repositories {
+dependencies {
->     maven("https://repo.kotlin.link")
+    implementation("space.kscience:kmath-functions:0.3.0-dev-3")
->     maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
+}
->     maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
+```
-> }
-> 
-> dependencies {
->     implementation("space.kscience:kmath-functions:0.3.0-dev-3")
-> }
-> ```

kmath-functions/docs/README-TEMPLATE.md

@@ -1,6 +1,6 @@
-# Functions (`kmath-functions`)
+# Module kmath-functions
 
-Functions and interpolations:
+Functions and interpolations.
 
 ${features}
 

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

@@ -2,14 +2,28 @@ package space.kscience.kmath.functions
 
 import space.kscience.kmath.operations.Ring
 
+/**
+ * Represents piecewise-defined function.
+ *
+ * @param T the piece key type.
+ * @param R the sub-function type.
+ */
 public fun interface Piecewise<T, R> {
+    /**
+     * Returns the appropriate sub-function for given piece key.
+     */
     public fun findPiece(arg: T): R?
 }
 
+/**
+ * Represents piecewise-defined function where all the sub-functions are polynomials.
+ */
 public fun interface PiecewisePolynomial<T : Any> : Piecewise<T, Polynomial<T>>
 
 /**
- * Ordered list of pieces in piecewise function
+ * Basic [Piecewise] implementation where all the pieces are ordered by the [Comparable] type instances.
+ *
+ * @param T the comparable piece key type.
  */
 public class OrderedPiecewisePolynomial<T : Comparable<T>>(delimiter: T) :
     PiecewisePolynomial<T> {
@@ -17,22 +31,30 @@ public class OrderedPiecewisePolynomial<T : Comparable<T>>(delimiter: T) :
     private val pieces: MutableList<Polynomial<T>> = arrayListOf()
 
     /**
-     * Dynamically add a piece to the "right" side (beyond maximum argument value of previous piece)
+     * Dynamically adds a piece to the right side (beyond maximum argument value of previous piece)
-     * @param right new rightmost position. If is less then current rightmost position, a error is thrown.
+     *
+     * @param right new rightmost position. If is less then current rightmost position, an error is thrown.
+     * @param piece the sub-function.
      */
     public fun putRight(right: T, piece: Polynomial<T>) {
         require(right > delimiters.last()) { "New delimiter should be to the right of old one" }
-        delimiters.add(right)
+        delimiters += right
-        pieces.add(piece)
+        pieces += piece
     }
 
+    /**
+     * Dynamically adds a piece to the left side (beyond maximum argument value of previous piece)
+     *
+     * @param left the new leftmost position. If is less then current rightmost position, an error is thrown.
+     * @param piece the sub-function.
+     */
     public fun putLeft(left: T, piece: Polynomial<T>) {
         require(left < delimiters.first()) { "New delimiter should be to the left of old one" }
         delimiters.add(0, left)
         pieces.add(0, piece)
     }
 
-    override fun findPiece(arg: T): Polynomial<T>? {
+    public override fun findPiece(arg: T): Polynomial<T>? {
         if (arg < delimiters.first() || arg >= delimiters.last())
             return null
         else {
@@ -46,7 +68,8 @@ public class OrderedPiecewisePolynomial<T : Comparable<T>>(delimiter: T) :
 }
 
 /**
- * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise definition.
+ * Return a value of polynomial function with given [ring] an given [arg] or null if argument is outside of piecewise
+ * definition.
  */
 public fun <T : Comparable<T>, C : Ring<T>> PiecewisePolynomial<T>.value(ring: C, arg: T): T? =
     findPiece(arg)?.value(ring, arg)

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

@@ -10,16 +10,26 @@ import kotlin.math.max
 import kotlin.math.pow
 
 /**
- * Polynomial coefficients without fixation on specific context they are applied to
+ * Polynomial coefficients model without fixation on specific context they are applied to.
- * @param coefficients constant is the leftmost coefficient
+ *
+ * @param coefficients constant is the leftmost coefficient.
  */
 public inline class Polynomial<T : Any>(public val coefficients: List<T>)
 
+/**
+ * Returns a [Polynomial] instance with given [coefficients].
+ */
 @Suppress("FunctionName")
 public fun <T : Any> Polynomial(vararg coefficients: T): Polynomial<T> = Polynomial(coefficients.toList())
 
+/**
+ * Evaluates the value of the given double polynomial for given double argument.
+ */
 public fun Polynomial<Double>.value(): Double = coefficients.reduceIndexed { index, acc, d -> acc + d.pow(index) }
 
+/**
+ * Evaluates the value of the given polynomial for given argument.
+ */
 public fun <T : Any, C : Ring<T>> Polynomial<T>.value(ring: C, arg: T): T = ring {
     if (coefficients.isEmpty()) return@ring zero
     var res = coefficients.first()
@@ -35,19 +45,23 @@ public fun <T : Any, C : Ring<T>> Polynomial<T>.value(ring: C, arg: T): T = ring
 }
 
 /**
- * Represent the polynomial as a regular context-less function
+ * Represent the polynomial as a regular context-less function.
  */
 public fun <T : Any, C : Ring<T>> Polynomial<T>.asFunction(ring: C): (T) -> T = { value(ring, it) }
 
 /**
- * An algebra for polynomials
+ * Space of polynomials.
+ *
+ * @param T the type of operated polynomials.
+ * @param C the intersection of [Ring] of [T] and [ScaleOperations] of [T].
+ * @param ring the [C] instance.
  */
 public class PolynomialSpace<T : Any, C>(
     private val ring: C,
 ) : Group<Polynomial<T>>, ScaleOperations<Polynomial<T>> where C : Ring<T>, C : ScaleOperations<T> {
     public override val zero: Polynomial<T> = Polynomial(emptyList())
 
-    override fun Polynomial<T>.unaryMinus(): Polynomial<T> = with(ring) {
+    override fun Polynomial<T>.unaryMinus(): Polynomial<T> = ring {
         Polynomial(coefficients.map { -it })
     }
 
@@ -64,6 +78,9 @@ public class PolynomialSpace<T : Any, C>(
     public override fun scale(a: Polynomial<T>, value: Double): Polynomial<T> =
         ring { Polynomial(List(a.coefficients.size) { index -> a.coefficients[index] * value }) }
 
+    /**
+     * Evaluates the polynomial for the given value [arg].
+     */
     public operator fun Polynomial<T>.invoke(arg: T): T = value(ring, arg)
 }
 

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/integration/Integrator.kt

@@ -1,11 +1,11 @@
 package space.kscience.kmath.integration
 
 /**
- * A general interface for all integrators
+ * A general interface for all integrators.
  */
 public interface Integrator<I : Integrand> {
     /**
-     * Run one integration pass and return a new [Integrand] with a new set of features
+     * Runs one integration pass and return a new [Integrand] with a new set of features.
      */
     public fun integrate(integrand: I): I
 }

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/Interpolator.kt

@@ -1,45 +1,50 @@
+@file:OptIn(UnstableKMathAPI::class)
+
 package space.kscience.kmath.interpolation
 
+import space.kscience.kmath.data.XYColumnarData
 import space.kscience.kmath.functions.PiecewisePolynomial
 import space.kscience.kmath.functions.value
+import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Ring
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.asBuffer
 
-public fun interface Interpolator<X, Y> {
+public fun interface Interpolator<T, X : T, Y : T> {
-    public fun interpolate(points: XYPointSet<X, Y>): (X) -> Y
+    public fun interpolate(points: XYColumnarData<T, X, Y>): (X) -> Y
 }
 
-public interface PolynomialInterpolator<T : Comparable<T>> : Interpolator<T, T> {
+public interface PolynomialInterpolator<T : Comparable<T>> : Interpolator<T, T, T> {
     public val algebra: Ring<T>
 
     public fun getDefaultValue(): T = error("Out of bounds")
 
-    public fun interpolatePolynomials(points: XYPointSet<T, T>): PiecewisePolynomial<T>
+    public fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T>
 
-    override fun interpolate(points: XYPointSet<T, T>): (T) -> T = { x ->
+    override fun interpolate(points: XYColumnarData<T, T, T>): (T) -> T = { x ->
         interpolatePolynomials(points).value(algebra, x) ?: getDefaultValue()
     }
 }
 
+
 public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
     x: Buffer<T>,
     y: Buffer<T>,
 ): PiecewisePolynomial<T> {
-    val pointSet = BufferXYPointSet(x, y)
+    val pointSet = XYColumnarData(x, y)
     return interpolatePolynomials(pointSet)
 }
 
 public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
     data: Map<T, T>,
 ): PiecewisePolynomial<T> {
-    val pointSet = BufferXYPointSet(data.keys.toList().asBuffer(), data.values.toList().asBuffer())
+    val pointSet = XYColumnarData(data.keys.toList().asBuffer(), data.values.toList().asBuffer())
     return interpolatePolynomials(pointSet)
 }
 
 public fun <T : Comparable<T>> PolynomialInterpolator<T>.interpolatePolynomials(
     data: List<Pair<T, T>>,
 ): PiecewisePolynomial<T> {
-    val pointSet = BufferXYPointSet(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer())
+    val pointSet = XYColumnarData(data.map { it.first }.asBuffer(), data.map { it.second }.asBuffer())
     return interpolatePolynomials(pointSet)
 }

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/LinearInterpolator.kt

@@ -1,16 +1,25 @@
 package space.kscience.kmath.interpolation
 
+import space.kscience.kmath.data.XYColumnarData
 import space.kscience.kmath.functions.OrderedPiecewisePolynomial
 import space.kscience.kmath.functions.PiecewisePolynomial
 import space.kscience.kmath.functions.Polynomial
+import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Field
 import space.kscience.kmath.operations.invoke
 
+@OptIn(UnstableKMathAPI::class)
+internal fun <T : Comparable<T>> insureSorted(points: XYColumnarData<*, T, *>) {
+    for (i in 0 until points.size - 1)
+        require(points.x[i + 1] > points.x[i]) { "Input data is not sorted at index $i" }
+}
+
 /**
  * Reference JVM implementation: https://github.com/apache/commons-math/blob/master/src/main/java/org/apache/commons/math4/analysis/interpolation/LinearInterpolator.java
  */
 public class LinearInterpolator<T : Comparable<T>>(public override val algebra: Field<T>) : PolynomialInterpolator<T> {
-    public override fun interpolatePolynomials(points: XYPointSet<T, T>): PiecewisePolynomial<T> = algebra {
+    @OptIn(UnstableKMathAPI::class)
+    public override fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T> = algebra {
         require(points.size > 0) { "Point array should not be empty" }
         insureSorted(points)
 

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/SplineInterpolator.kt

@@ -1,15 +1,20 @@
 package space.kscience.kmath.interpolation
 
+import space.kscience.kmath.data.XYColumnarData
 import space.kscience.kmath.functions.OrderedPiecewisePolynomial
 import space.kscience.kmath.functions.PiecewisePolynomial
 import space.kscience.kmath.functions.Polynomial
+import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Field
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.MutableBufferFactory
 
 /**
- * Generic spline interpolator. Not recommended for performance critical places, use platform-specific and type specific ones.
+ * Generic spline interpolator. Not recommended for performance critical places, use platform-specific and type
- * Based on https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/SplineInterpolator.java
+ * specific ones.
+ *
+ * Based on
+ * https://github.com/apache/commons-math/blob/eb57d6d457002a0bb5336d789a3381a24599affe/src/main/java/org/apache/commons/math4/analysis/interpolation/SplineInterpolator.java
  */
 public class SplineInterpolator<T : Comparable<T>>(
     public override val algebra: Field<T>,
@@ -17,7 +22,8 @@ public class SplineInterpolator<T : Comparable<T>>(
 ) : PolynomialInterpolator<T> {
     //TODO possibly optimize zeroed buffers
 
-    public override fun interpolatePolynomials(points: XYPointSet<T, T>): PiecewisePolynomial<T> = algebra {
+    @OptIn(UnstableKMathAPI::class)
+    public override fun interpolatePolynomials(points: XYColumnarData<T, T, T>): PiecewisePolynomial<T> = algebra {
         require(points.size >= 3) { "Can't use spline interpolator with less than 3 points" }
         insureSorted(points)
         // Number of intervals.  The number of data points is n + 1.

kmath-functions/src/commonMain/kotlin/space/kscience/kmath/interpolation/XYPointSet.kt

@@ -1,53 +0,0 @@
-package space.kscience.kmath.interpolation
-
-import space.kscience.kmath.nd.Structure2D
-import space.kscience.kmath.structures.Buffer
-
-public interface XYPointSet<X, Y> {
-    public val size: Int
-    public val x: Buffer<X>
-    public val y: Buffer<Y>
-}
-
-public interface XYZPointSet<X, Y, Z> : XYPointSet<X, Y> {
-    public val z: Buffer<Z>
-}
-
-internal fun <T : Comparable<T>> insureSorted(points: XYPointSet<T, *>) {
-    for (i in 0 until points.size - 1)
-        require(points.x[i + 1] > points.x[i]) { "Input data is not sorted at index $i" }
-}
-
-public class NDStructureColumn<T>(public val structure: Structure2D<T>, public val column: Int) : Buffer<T> {
-    public override val size: Int
-        get() = structure.rowNum
-
-    init {
-        require(column < structure.colNum) { "Column index is outside of structure column range" }
-    }
-
-    public override operator fun get(index: Int): T = structure[index, column]
-    public override operator fun iterator(): Iterator<T> = sequence { repeat(size) { yield(get(it)) } }.iterator()
-}
-
-public class BufferXYPointSet<X, Y>(
-    public override val x: Buffer<X>,
-    public override val y: Buffer<Y>,
-) : XYPointSet<X, Y> {
-    public override val size: Int
-        get() = x.size
-
-    init {
-        require(x.size == y.size) { "Sizes of x and y buffers should be the same" }
-    }
-}
-
-public fun <T> Structure2D<T>.asXYPointSet(): XYPointSet<T, T> {
-    require(shape[1] == 2) { "Structure second dimension should be of size 2" }
-
-    return object : XYPointSet<T, T> {
-        override val size: Int get() = this@asXYPointSet.shape[0]
-        override val x: Buffer<T> get() = NDStructureColumn(this@asXYPointSet, 0)
-        override val y: Buffer<T> get() = NDStructureColumn(this@asXYPointSet, 1)
-    }
-}

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

@@ -5,7 +5,7 @@ import space.kscience.kmath.ast.MST
 import space.kscience.kmath.ast.MstAlgebra
 import space.kscience.kmath.ast.MstExpression
 import space.kscience.kmath.expressions.DifferentiableExpression
-import space.kscience.kmath.expressions.Symbol
+import space.kscience.kmath.misc.Symbol
 import space.kscience.kmath.operations.NumericAlgebra
 
 /**
@@ -18,8 +18,10 @@ import space.kscience.kmath.operations.NumericAlgebra
  * @param A the [NumericAlgebra] of [T].
  * @property expr the underlying [MstExpression].
  */
-public inline class DifferentiableMstExpression<T, A>(public val expr: MstExpression<T, A>) :
+public inline class DifferentiableMstExpression<T: Number, A>(
-    DifferentiableExpression<T, MstExpression<T, A>> where A : NumericAlgebra<T>, T : Number {
+    public val expr: MstExpression<T, A>,
+) : DifferentiableExpression<T, MstExpression<T, A>> where A : NumericAlgebra<T> {
+
     public constructor(algebra: A, mst: MST) : this(MstExpression(algebra, mst))
 
     /**

kmath-nd4j/README.md

@@ -1,46 +1,40 @@
-# ND4J NDStructure implementation (`kmath-nd4j`)
+# Module kmath-nd4j
 
-This subproject implements the following features:
+ND4J based implementations of KMath abstractions.
 
  - [nd4jarraystructure](#) : NDStructure wrapper for INDArray
  - [nd4jarrayrings](#) : Rings over Nd4jArrayStructure of Int and Long
  - [nd4jarrayfields](#) : Fields over Nd4jArrayStructure of Float and Double
 
 
-> #### Artifact:
+## Artifact:
->
+
-> This module artifact: `space.kscience:kmath-nd4j:0.3.0-dev-3`.
+The Maven coordinates of this project are `space.kscience:kmath-nd4j:0.3.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)
+**Gradle:**
->
+```gradle
-> 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)
+repositories {
->
+    maven { url 'https://repo.kotlin.link' }
-> **Gradle:**
+    maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
->
+    maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
-> ```gradle
+}
-> repositories {
+
->     maven { url 'https://repo.kotlin.link' }
+dependencies {
->     maven { url 'https://dl.bintray.com/hotkeytlt/maven' }
+    implementation 'space.kscience:kmath-nd4j:0.3.0-dev-3'
->     maven { url "https://dl.bintray.com/kotlin/kotlin-eap" } // include for builds based on kotlin-eap
+}
-> }
+```
-> 
+**Gradle Kotlin DSL:**
-> dependencies {
+```kotlin
->     implementation 'space.kscience:kmath-nd4j:0.3.0-dev-3'
+repositories {
-> }
+    maven("https://repo.kotlin.link")
-> ```
+    maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
-> **Gradle Kotlin DSL:**
+    maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
->
+}
-> ```kotlin
+
-> repositories {
+dependencies {
->     maven("https://repo.kotlin.link")
+    implementation("space.kscience:kmath-nd4j:0.3.0-dev-3")
->     maven("https://dl.bintray.com/kotlin/kotlin-eap") // include for builds based on kotlin-eap
+}
->     maven("https://dl.bintray.com/hotkeytlt/maven") // required for a
+```
-> }
-> 
-> dependencies {
->     implementation("space.kscience:kmath-nd4j:0.3.0-dev-3")
-> }
-> ```
 
 ## Examples
 

kmath-nd4j/docs/README-TEMPLATE.md

@@ -1,6 +1,6 @@
-# ND4J NDStructure implementation (`kmath-nd4j`)
+# Module kmath-nd4j
 
-This subproject implements the following features:
+ND4J based implementations of KMath abstractions.
 
 ${features}
 

kmath-stat/build.gradle.kts

@@ -3,14 +3,6 @@ plugins {
 }
 
 kotlin.sourceSets {
-    all {
-        languageSettings.apply {
-            useExperimentalAnnotation("kotlinx.coroutines.FlowPreview")
-            useExperimentalAnnotation("kotlinx.coroutines.ExperimentalCoroutinesApi")
-            useExperimentalAnnotation("kotlinx.coroutines.ObsoleteCoroutinesApi")
-        }
-    }
-
     commonMain {
         dependencies {
             api(project(":kmath-coroutines"))

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/FunctionOptimization.kt

@@ -0,0 +1,89 @@
+package space.kscience.kmath.optimization
+
+import space.kscience.kmath.expressions.AutoDiffProcessor
+import space.kscience.kmath.expressions.DifferentiableExpression
+import space.kscience.kmath.expressions.Expression
+import space.kscience.kmath.expressions.ExpressionAlgebra
+import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.operations.ExtendedField
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.indices
+
+/**
+ * A likelihood function optimization problem with provided derivatives
+ */
+public interface FunctionOptimization<T : Any> : Optimization<T> {
+    /**
+     * The optimization direction. If true search for function maximum, if false, search for the minimum
+     */
+    public var maximize: Boolean
+
+    /**
+     * Define the initial guess for the optimization problem
+     */
+    public fun initialGuess(map: Map<Symbol, T>)
+
+    /**
+     * Set a differentiable expression as objective function as function and gradient provider
+     */
+    public fun diffFunction(expression: DifferentiableExpression<T, Expression<T>>)
+
+    public companion object {
+        /**
+         * 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
+            }
+        }
+    }
+}
+
+/**
+ * Define a chi-squared-based objective function
+ */
+public fun <T: Any, I : Any, A> FunctionOptimization<T>.chiSquared(
+    autoDiff: AutoDiffProcessor<T, I, A, Expression<T>>,
+    x: Buffer<T>,
+    y: Buffer<T>,
+    yErr: Buffer<T>,
+    model: A.(I) -> I,
+) where A : ExtendedField<I>, A : ExpressionAlgebra<T, I> {
+    val chiSquared = FunctionOptimization.chiSquared(autoDiff, x, y, yErr, model)
+    diffFunction(chiSquared)
+    maximize = false
+}
+
+/**
+ * Optimize differentiable expression using specific [OptimizationProblemFactory]
+ */
+public fun <T : Any, F : FunctionOptimization<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.diffFunction(this)
+    return problem.optimize()
+}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/NoDerivFunctionOptimization.kt

@@ -0,0 +1,69 @@
+package space.kscience.kmath.optimization
+
+import space.kscience.kmath.expressions.Expression
+import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.indices
+import kotlin.math.pow
+
+/**
+ * A likelihood function optimization problem
+ */
+public interface NoDerivFunctionOptimization<T : Any> : Optimization<T> {
+    /**
+     * The optimization direction. If true search for function maximum, if false, search for the minimum
+     */
+    public var maximize: Boolean
+
+    /**
+     * Define the initial guess for the optimization problem
+     */
+    public fun initialGuess(map: Map<Symbol, T>)
+
+    /**
+     * Set an objective function expression
+     */
+    public fun function(expression: Expression<T>)
+
+    public companion object {
+        /**
+         * 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 = Symbol.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)
+                }
+            }
+        }
+    }
+}
+
+
+/**
+ * Optimize expression without derivatives using specific [OptimizationProblemFactory]
+ */
+public fun <T : Any, F : NoDerivFunctionOptimization<T>> Expression<T>.noDerivOptimizeWith(
+    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.function(this)
+    return problem.optimize()
+}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/Optimization.kt

@@ -0,0 +1,44 @@
+package space.kscience.kmath.optimization
+
+import space.kscience.kmath.misc.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)
+
+/**
+ * An optimization problem builder over [T] variables
+ */
+public interface Optimization<T : Any> {
+
+    /**
+     * 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 : Optimization<T>> {
+    public fun build(symbols: List<Symbol>): P
+}
+
+public operator fun <T : Any, P : Optimization<T>> OptimizationProblemFactory<T, P>.invoke(
+    symbols: List<Symbol>,
+    block: P.() -> Unit,
+): P = build(symbols).apply(block)

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/optimization/XYFit.kt

@@ -0,0 +1,40 @@
+package space.kscience.kmath.optimization
+
+import space.kscience.kmath.data.ColumnarData
+import space.kscience.kmath.expressions.AutoDiffProcessor
+import space.kscience.kmath.expressions.DifferentiableExpression
+import space.kscience.kmath.expressions.Expression
+import space.kscience.kmath.expressions.ExpressionAlgebra
+import space.kscience.kmath.misc.Symbol
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.ExtendedField
+import space.kscience.kmath.operations.Field
+
+@UnstableKMathAPI
+public interface XYFit<T : Any> : Optimization<T> {
+
+    public val algebra: Field<T>
+
+    /**
+     * Set X-Y data for this fit optionally including x and y errors
+     */
+    public fun data(
+        dataSet: ColumnarData<T>,
+        xSymbol: Symbol,
+        ySymbol: Symbol,
+        xErrSymbol: Symbol? = null,
+        yErrSymbol: Symbol? = null,
+    )
+
+    public fun model(model: (T) -> DifferentiableExpression<T, *>)
+
+    /**
+     * Set the differentiable model for this fit
+     */
+    public fun <I : Any, A> model(
+        autoDiff: AutoDiffProcessor<T, I, A, Expression<T>>,
+        modelFunction: A.(I) -> I,
+    ): Unit where A : ExtendedField<I>, A : ExpressionAlgebra<T, I> = model { arg ->
+        autoDiff.process { modelFunction(const(arg)) }
+    }
+}

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

@@ -1,17 +1,29 @@
 package space.kscience.kmath.stat
 
+import kotlinx.coroutines.flow.first
 import space.kscience.kmath.chains.Chain
 import space.kscience.kmath.chains.collect
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.BufferFactory
-import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.IntBuffer
+import space.kscience.kmath.structures.MutableBuffer
+import kotlin.jvm.JvmName
 
-public interface Sampler<T : Any> {
+/**
+ * Sampler that generates chains of values of type [T].
+ */
+public fun interface Sampler<T : Any> {
+    /**
+     * Generates a chain of samples.
+     *
+     * @param generator the randomness provider.
+     * @return the new chain.
+     */
     public fun sample(generator: RandomGenerator): Chain<T>
 }
 
 /**
- * A distribution of typed objects
+ * A distribution of typed objects.
  */
 public interface Distribution<T : Any> : Sampler<T> {
     /**
@@ -20,11 +32,7 @@ public interface Distribution<T : Any> : Sampler<T> {
      */
     public fun probability(arg: T): Double
 
-    /**
+    public override fun sample(generator: RandomGenerator): Chain<T>
-     * Create a chain of samples from this distribution.
-     * The chain is not guaranteed to be stateless, but different sample chains should be independent.
-     */
-    override fun sample(generator: RandomGenerator): Chain<T>
 
     /**
      * An empty companion. Distribution factories should be written as its extensions
@@ -63,16 +71,27 @@ public fun <T : Any> Sampler<T>.sampleBuffer(
         //clear list from previous run
         tmp.clear()
         //Fill list
-        repeat(size) {
+        repeat(size) { tmp += chain.next() }
-            tmp.add(chain.next())
-        }
         //return new buffer with elements from tmp
         bufferFactory(size) { tmp[it] }
     }
 }
 
 /**
- * Generate a bunch of samples from real distributions
+ * Samples one value from this [Sampler].
  */
+public suspend fun <T : Any> Sampler<T>.next(generator: RandomGenerator): T = sample(generator).first()
+
+/**
+ * Generates [size] real samples and chunks them into some buffers.
+ */
+@JvmName("sampleRealBuffer")
 public fun Sampler<Double>.sampleBuffer(generator: RandomGenerator, size: Int): Chain<Buffer<Double>> =
-    sampleBuffer(generator, size, ::DoubleBuffer)
+    sampleBuffer(generator, size, MutableBuffer.Companion::double)
+
+/**
+ * Generates [size] integer samples and chunks them into some buffers.
+ */
+@JvmName("sampleIntBuffer")
+public fun Sampler<Int>.sampleBuffer(generator: RandomGenerator, size: Int): Chain<Buffer<Int>> =
+    sampleBuffer(generator, size, ::IntBuffer)

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

@@ -14,7 +14,7 @@ public interface NamedDistribution<T> : Distribution<Map<String, T>>
 public class FactorizedDistribution<T>(public val distributions: Collection<NamedDistribution<T>>) :
     NamedDistribution<T> {
     override fun probability(arg: Map<String, T>): Double =
-        distributions.fold(1.0) { acc, distr -> acc * distr.probability(arg) }
+        distributions.fold(1.0) { acc, dist -> acc * dist.probability(arg) }
 
     override fun sample(generator: RandomGenerator): Chain<Map<String, T>> {
         val chains = distributions.map { it.sample(generator) }
@@ -38,6 +38,6 @@ public class DistributionBuilder<T : Any> {
     private val distributions = ArrayList<NamedDistribution<T>>()
 
     public infix fun String.to(distribution: Distribution<T>) {
-        distributions.add(NamedDistributionWrapper(this, distribution))
+        distributions += NamedDistributionWrapper(this, distribution)
     }
 }

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

@@ -1,63 +0,0 @@
-package space.kscience.kmath.stat
-
-import space.kscience.kmath.expressions.*
-import space.kscience.kmath.operations.ExtendedField
-import space.kscience.kmath.structures.Buffer
-import space.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/space/kscience/kmath/stat/OptimizationProblem.kt

@@ -1,88 +0,0 @@
-package space.kscience.kmath.stat
-
-import space.kscience.kmath.expressions.DifferentiableExpression
-import space.kscience.kmath.expressions.Expression
-import space.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/space/kscience/kmath/stat/RandomChain.kt

@@ -1,17 +1,22 @@
 package space.kscience.kmath.stat
 
+import space.kscience.kmath.chains.BlockingDoubleChain
+import space.kscience.kmath.chains.BlockingIntChain
 import space.kscience.kmath.chains.Chain
 
 /**
  * A possibly stateful chain producing random values.
+ *
+ * @property generator the underlying [RandomGenerator] instance.
  */
 public class RandomChain<out R>(
     public val generator: RandomGenerator,
-    private val gen: suspend RandomGenerator.() -> R,
+    private val gen: suspend RandomGenerator.() -> R
 ) : Chain<R> {
     override suspend fun next(): R = generator.gen()
-
     override fun fork(): Chain<R> = RandomChain(generator.fork(), gen)
 }
 
 public fun <R> RandomGenerator.chain(gen: suspend RandomGenerator.() -> R): RandomChain<R> = RandomChain(this, gen)
+public fun Chain<Double>.blocking(): BlockingDoubleChain = object : Chain<Double> by this, BlockingDoubleChain {}
+public fun Chain<Int>.blocking(): BlockingIntChain = object : Chain<Int> by this, BlockingIntChain {}

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

@@ -82,6 +82,8 @@ public interface RandomGenerator {
 
 /**
  * Implements [RandomGenerator] by delegating all operations to [Random].
+ *
+ * @property random the underlying [Random] object.
  */
 public class DefaultGenerator(public val random: Random = Random) : RandomGenerator {
     public override fun nextBoolean(): Boolean = random.nextBoolean()

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

@@ -8,16 +8,28 @@ import space.kscience.kmath.operations.Group
 import space.kscience.kmath.operations.ScaleOperations
 import space.kscience.kmath.operations.invoke
 
-public class BasicSampler<T : Any>(public val chainBuilder: (RandomGenerator) -> Chain<T>) : Sampler<T> {
+/**
-    public override fun sample(generator: RandomGenerator): Chain<T> = chainBuilder(generator)
+ * Implements [Sampler] by sampling only certain [value].
-}
+ *
-
+ * @property value the value to sample.
+ */
 public class ConstantSampler<T : Any>(public val value: T) : Sampler<T> {
     public override fun sample(generator: RandomGenerator): Chain<T> = ConstantChain(value)
 }
 
 /**
- * A space for samplers. Allows to perform simple operations on distributions
+ * Implements [Sampler] by delegating sampling to value of [chainBuilder].
+ *
+ * @property chainBuilder the provider of [Chain].
+ */
+public class BasicSampler<T : Any>(public val chainBuilder: (RandomGenerator) -> Chain<T>) : Sampler<T> {
+    public override fun sample(generator: RandomGenerator): Chain<T> = chainBuilder(generator)
+}
+
+/**
+ * A space of samplers. Allows to perform simple operations on distributions.
+ *
+ * @property algebra the space to provide addition and scalar multiplication for [T].
  */
 public class SamplerSpace<T : Any, S>(public val algebra: S) : Group<Sampler<T>>,
     ScaleOperations<Sampler<T>> where S : Group<T>, S : ScaleOperations<T> {
@@ -29,8 +41,10 @@ public class SamplerSpace<T : Any, S>(public val algebra: S) : Group<Sampler<T>>
     }
 
     public override fun scale(a: Sampler<T>, value: Double): Sampler<T> = BasicSampler { generator ->
-        a.sample(generator).map { algebra { it * value } }
+        a.sample(generator).map { a ->
+            algebra { a * value }
+        }
     }
 
-    override fun Sampler<T>.unaryMinus(): Sampler<T> = scale(this, -1.0)
+    public override fun Sampler<T>.unaryMinus(): Sampler<T> = scale(this, -1.0)
 }

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/distributions/NormalDistribution.kt

@@ -0,0 +1,41 @@
+package space.kscience.kmath.stat.distributions
+
+import space.kscience.kmath.chains.Chain
+import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.stat.UnivariateDistribution
+import space.kscience.kmath.stat.internal.InternalErf
+import space.kscience.kmath.stat.samplers.GaussianSampler
+import space.kscience.kmath.stat.samplers.NormalizedGaussianSampler
+import space.kscience.kmath.stat.samplers.ZigguratNormalizedGaussianSampler
+import kotlin.math.*
+
+/**
+ * Implements [UnivariateDistribution] for the normal (gaussian) distribution.
+ */
+public inline class NormalDistribution(public val sampler: GaussianSampler) : UnivariateDistribution<Double> {
+    public constructor(
+        mean: Double,
+        standardDeviation: Double,
+        normalized: NormalizedGaussianSampler = ZigguratNormalizedGaussianSampler.of(),
+    ) : this(GaussianSampler.of(mean, standardDeviation, normalized))
+
+    public override fun probability(arg: Double): Double {
+        val x1 = (arg - sampler.mean) / sampler.standardDeviation
+        return exp(-0.5 * x1 * x1 - (ln(sampler.standardDeviation) + 0.5 * ln(2 * PI)))
+    }
+
+    public override fun sample(generator: RandomGenerator): Chain<Double> = sampler.sample(generator)
+
+    public override fun cumulative(arg: Double): Double {
+        val dev = arg - sampler.mean
+
+        return when {
+            abs(dev) > 40 * sampler.standardDeviation -> if (dev < 0) 0.0 else 1.0
+            else -> 0.5 * InternalErf.erfc(-dev / (sampler.standardDeviation * SQRT2))
+        }
+    }
+
+    private companion object {
+        private val SQRT2 = sqrt(2.0)
+    }
+}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalErf.kt

@@ -0,0 +1,15 @@
+package space.kscience.kmath.stat.internal
+
+import kotlin.math.abs
+
+/**
+ * Based on Commons Math implementation.
+ * See [https://commons.apache.org/proper/commons-math/javadocs/api-3.3/org/apache/commons/math3/special/Erf.html].
+ */
+internal object InternalErf {
+    fun erfc(x: Double): Double {
+        if (abs(x) > 40) return if (x > 0) 0.0 else 2.0
+        val ret = InternalGamma.regularizedGammaQ(0.5, x * x, 10000)
+        return if (x < 0) 2 - ret else ret
+    }
+}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalGamma.kt

@@ -0,0 +1,238 @@
+package space.kscience.kmath.stat.internal
+
+import kotlin.math.*
+
+private abstract class ContinuedFraction protected constructor() {
+    protected abstract fun getA(n: Int, x: Double): Double
+    protected abstract fun getB(n: Int, x: Double): Double
+
+    fun evaluate(x: Double, maxIterations: Int): Double {
+        val small = 1e-50
+        var hPrev = getA(0, x)
+        if (hPrev == 0.0 || abs(0.0 - hPrev) <= small) hPrev = small
+        var n = 1
+        var dPrev = 0.0
+        var cPrev = hPrev
+        var hN = hPrev
+
+        while (n < maxIterations) {
+            val a = getA(n, x)
+            val b = getB(n, x)
+            var dN = a + b * dPrev
+            if (dN == 0.0 || abs(0.0 - dN) <= small) dN = small
+            var cN = a + b / cPrev
+            if (cN == 0.0 || abs(0.0 - cN) <= small) cN = small
+            dN = 1 / dN
+            val deltaN = cN * dN
+            hN = hPrev * deltaN
+            check(!hN.isInfinite()) { "hN is infinite" }
+            check(!hN.isNaN()) { "hN is NaN" }
+            if (abs(deltaN - 1.0) < 10e-9) break
+            dPrev = dN
+            cPrev = cN
+            hPrev = hN
+            n++
+        }
+
+        check(n < maxIterations) { "n is more than maxIterations" }
+        return hN
+    }
+}
+
+internal object InternalGamma {
+    const val LANCZOS_G = 607.0 / 128.0
+
+    private val LANCZOS = doubleArrayOf(
+        0.99999999999999709182,
+        57.156235665862923517,
+        -59.597960355475491248,
+        14.136097974741747174,
+        -0.49191381609762019978,
+        .33994649984811888699e-4,
+        .46523628927048575665e-4,
+        -.98374475304879564677e-4,
+        .15808870322491248884e-3,
+        -.21026444172410488319e-3,
+        .21743961811521264320e-3,
+        -.16431810653676389022e-3,
+        .84418223983852743293e-4,
+        -.26190838401581408670e-4,
+        .36899182659531622704e-5
+    )
+
+    private val HALF_LOG_2_PI = 0.5 * ln(2.0 * PI)
+    private const val INV_GAMMA1P_M1_A0 = .611609510448141581788E-08
+    private const val INV_GAMMA1P_M1_A1 = .624730830116465516210E-08
+    private const val INV_GAMMA1P_M1_B1 = .203610414066806987300E+00
+    private const val INV_GAMMA1P_M1_B2 = .266205348428949217746E-01
+    private const val INV_GAMMA1P_M1_B3 = .493944979382446875238E-03
+    private const val INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05
+    private const val INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05
+    private const val INV_GAMMA1P_M1_B6 = .992641840672773722196E-06
+    private const val INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07
+    private const val INV_GAMMA1P_M1_B8 = .195755836614639731882E-09
+    private const val INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08
+    private const val INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08
+    private const val INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09
+    private const val INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10
+    private const val INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11
+    private const val INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12
+    private const val INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14
+    private const val INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00
+    private const val INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01
+    private const val INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02
+    private const val INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03
+    private const val INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00
+    private const val INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00
+    private const val INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00
+    private const val INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01
+    private const val INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00
+    private const val INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01
+    private const val INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02
+    private const val INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02
+    private const val INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02
+    private const val INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03
+    private const val INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03
+    private const val INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04
+    private const val INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05
+    private const val INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05
+    private const val INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06
+
+    fun logGamma(x: Double): Double = when {
+        x.isNaN() || x <= 0.0 -> Double.NaN
+        x < 0.5 -> logGamma1p(x) - ln(x)
+        x <= 2.5 -> logGamma1p(x - 0.5 - 0.5)
+
+        x <= 8.0 -> {
+            val n = floor(x - 1.5).toInt()
+            val prod = (1..n).fold(1.0, { prod, i -> prod * (x - i) })
+            logGamma1p(x - (n + 1)) + ln(prod)
+        }
+
+        else -> {
+            val tmp = x + LANCZOS_G + .5
+            (x + .5) * ln(tmp) - tmp + HALF_LOG_2_PI + ln(lanczos(x) / x)
+        }
+    }
+
+    private fun regularizedGammaP(
+        a: Double,
+        x: Double,
+        maxIterations: Int = Int.MAX_VALUE
+    ): Double = when {
+        a.isNaN() || x.isNaN() || a <= 0.0 || x < 0.0 -> Double.NaN
+        x == 0.0 -> 0.0
+        x >= a + 1 -> 1.0 - regularizedGammaQ(a, x, maxIterations)
+
+        else -> {
+            // calculate series
+            var n = 0.0 // current element index
+            var an = 1.0 / a // n-th element in the series
+            var sum = an // partial sum
+
+            while (abs(an / sum) > 10e-15 && n < maxIterations && sum < Double.POSITIVE_INFINITY) {
+                // compute next element in the series
+                n += 1.0
+                an *= x / (a + n)
+
+                // update partial sum
+                sum += an
+            }
+
+            when {
+                n >= maxIterations -> throw error("Maximal iterations is exceeded $maxIterations")
+                sum.isInfinite() -> 1.0
+                else -> exp(-x + a * ln(x) - logGamma(a)) * sum
+            }
+        }
+    }
+
+    fun regularizedGammaQ(
+        a: Double,
+        x: Double,
+        maxIterations: Int = Int.MAX_VALUE
+    ): Double = when {
+        a.isNaN() || x.isNaN() || a <= 0.0 || x < 0.0 -> Double.NaN
+        x == 0.0 -> 1.0
+        x < a + 1.0 -> 1.0 - regularizedGammaP(a, x, maxIterations)
+
+        else -> 1.0 / object : ContinuedFraction() {
+            override fun getA(n: Int, x: Double): Double = 2.0 * n + 1.0 - a + x
+            override fun getB(n: Int, x: Double): Double = n * (a - n)
+        }.evaluate(x, maxIterations) * exp(-x + a * ln(x) - logGamma(a))
+    }
+
+    private fun lanczos(x: Double): Double =
+        (LANCZOS.size - 1 downTo 1).sumByDouble { LANCZOS[it] / (x + it) } + LANCZOS[0]
+
+    private fun invGamma1pm1(x: Double): Double {
+        require(x >= -0.5)
+        require(x <= 1.5)
+        val ret: Double
+        val t = if (x <= 0.5) x else x - 0.5 - 0.5
+
+        if (t < 0.0) {
+            val a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1
+            var b = INV_GAMMA1P_M1_B8
+            b = INV_GAMMA1P_M1_B7 + t * b
+            b = INV_GAMMA1P_M1_B6 + t * b
+            b = INV_GAMMA1P_M1_B5 + t * b
+            b = INV_GAMMA1P_M1_B4 + t * b
+            b = INV_GAMMA1P_M1_B3 + t * b
+            b = INV_GAMMA1P_M1_B2 + t * b
+            b = INV_GAMMA1P_M1_B1 + t * b
+            b = 1.0 + t * b
+            var c = INV_GAMMA1P_M1_C13 + t * (a / b)
+            c = INV_GAMMA1P_M1_C12 + t * c
+            c = INV_GAMMA1P_M1_C11 + t * c
+            c = INV_GAMMA1P_M1_C10 + t * c
+            c = INV_GAMMA1P_M1_C9 + t * c
+            c = INV_GAMMA1P_M1_C8 + t * c
+            c = INV_GAMMA1P_M1_C7 + t * c
+            c = INV_GAMMA1P_M1_C6 + t * c
+            c = INV_GAMMA1P_M1_C5 + t * c
+            c = INV_GAMMA1P_M1_C4 + t * c
+            c = INV_GAMMA1P_M1_C3 + t * c
+            c = INV_GAMMA1P_M1_C2 + t * c
+            c = INV_GAMMA1P_M1_C1 + t * c
+            c = INV_GAMMA1P_M1_C + t * c
+            ret = (if (x > 0.5) t * c / x else x * (c + 0.5 + 0.5))
+        } else {
+            var p = INV_GAMMA1P_M1_P6
+            p = INV_GAMMA1P_M1_P5 + t * p
+            p = INV_GAMMA1P_M1_P4 + t * p
+            p = INV_GAMMA1P_M1_P3 + t * p
+            p = INV_GAMMA1P_M1_P2 + t * p
+            p = INV_GAMMA1P_M1_P1 + t * p
+            p = INV_GAMMA1P_M1_P0 + t * p
+            var q = INV_GAMMA1P_M1_Q4
+            q = INV_GAMMA1P_M1_Q3 + t * q
+            q = INV_GAMMA1P_M1_Q2 + t * q
+            q = INV_GAMMA1P_M1_Q1 + t * q
+            q = 1.0 + t * q
+            var c = INV_GAMMA1P_M1_C13 + p / q * t
+            c = INV_GAMMA1P_M1_C12 + t * c
+            c = INV_GAMMA1P_M1_C11 + t * c
+            c = INV_GAMMA1P_M1_C10 + t * c
+            c = INV_GAMMA1P_M1_C9 + t * c
+            c = INV_GAMMA1P_M1_C8 + t * c
+            c = INV_GAMMA1P_M1_C7 + t * c
+            c = INV_GAMMA1P_M1_C6 + t * c
+            c = INV_GAMMA1P_M1_C5 + t * c
+            c = INV_GAMMA1P_M1_C4 + t * c
+            c = INV_GAMMA1P_M1_C3 + t * c
+            c = INV_GAMMA1P_M1_C2 + t * c
+            c = INV_GAMMA1P_M1_C1 + t * c
+            c = INV_GAMMA1P_M1_C0 + t * c
+            ret = (if (x > 0.5) t / x * (c - 0.5 - 0.5) else x * c)
+        }
+
+        return ret
+    }
+
+    private fun logGamma1p(x: Double): Double {
+        require(x >= -0.5)
+        require(x <= 1.5)
+        return -ln1p(invGamma1pm1(x))
+    }
+}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/internal/InternalUtils.kt

@@ -0,0 +1,70 @@
+package space.kscience.kmath.stat.internal
+
+import kotlin.math.ln
+import kotlin.math.min
+
+internal object InternalUtils {
+    private val FACTORIALS = longArrayOf(
+        1L, 1L, 2L,
+        6L, 24L, 120L,
+        720L, 5040L, 40320L,
+        362880L, 3628800L, 39916800L,
+        479001600L, 6227020800L, 87178291200L,
+        1307674368000L, 20922789888000L, 355687428096000L,
+        6402373705728000L, 121645100408832000L, 2432902008176640000L
+    )
+
+    private const val BEGIN_LOG_FACTORIALS = 2
+
+    fun factorial(n: Int): Long = FACTORIALS[n]
+
+    fun validateProbabilities(probabilities: DoubleArray?): Double {
+        require(!(probabilities == null || probabilities.isEmpty())) { "Probabilities must not be empty." }
+
+        val sumProb = probabilities.sumByDouble { prob ->
+            require(!(prob < 0 || prob.isInfinite() || prob.isNaN())) { "Invalid probability: $prob" }
+            prob
+        }
+
+        require(!(sumProb.isInfinite() || sumProb <= 0)) { "Invalid sum of probabilities: $sumProb" }
+        return sumProb
+    }
+
+    class FactorialLog private constructor(numValues: Int, cache: DoubleArray?) {
+        private val logFactorials: DoubleArray = DoubleArray(numValues)
+
+        init {
+            val endCopy: Int
+
+            if (cache != null && cache.size > BEGIN_LOG_FACTORIALS) {
+                // Copy available values.
+                endCopy = min(cache.size, numValues)
+
+                cache.copyInto(
+                    logFactorials,
+                    BEGIN_LOG_FACTORIALS,
+                    BEGIN_LOG_FACTORIALS, endCopy
+                )
+            } else
+            // All values to be computed
+                endCopy = BEGIN_LOG_FACTORIALS
+
+            // Compute remaining values.
+            (endCopy until numValues).forEach { i ->
+                if (i < FACTORIALS.size)
+                    logFactorials[i] = ln(FACTORIALS[i].toDouble())
+                else
+                    logFactorials[i] = logFactorials[i - 1] + ln(i.toDouble())
+            }
+        }
+
+        fun value(n: Int): Double {
+            if (n < logFactorials.size) return logFactorials[n]
+            return if (n < FACTORIALS.size) ln(FACTORIALS[n].toDouble()) else InternalGamma.logGamma(n + 1.0)
+        }
+
+        companion object {
+            fun create(): FactorialLog = FactorialLog(0, null)
+        }
+    }
+}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AhrensDieterExponentialSampler.kt

@@ -0,0 +1,73 @@
+package space.kscience.kmath.stat.samplers
+
+import space.kscience.kmath.chains.Chain
+import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.stat.Sampler
+import space.kscience.kmath.stat.chain
+import space.kscience.kmath.stat.internal.InternalUtils
+import kotlin.math.ln
+import kotlin.math.pow
+
+/**
+ * Sampling from an [exponential distribution](http://mathworld.wolfram.com/ExponentialDistribution.html).
+ *
+ * Based on Commons RNG implementation.
+ * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/AhrensDieterExponentialSampler.html].
+ */
+public class AhrensDieterExponentialSampler private constructor(public val mean: Double) : Sampler<Double> {
+    public override fun sample(generator: RandomGenerator): Chain<Double> = generator.chain {
+        // Step 1:
+        var a = 0.0
+        var u = nextDouble()
+
+        // Step 2 and 3:
+        while (u < 0.5) {
+            a += EXPONENTIAL_SA_QI[0]
+            u *= 2.0
+        }
+
+        // Step 4 (now u >= 0.5):
+        u += u - 1
+        // Step 5:
+        if (u <= EXPONENTIAL_SA_QI[0]) return@chain mean * (a + u)
+        // Step 6:
+        var i = 0 // Should be 1, be we iterate before it in while using 0.
+        var u2 = nextDouble()
+        var umin = u2
+
+        // Step 7 and 8:
+        do {
+            ++i
+            u2 = nextDouble()
+            if (u2 < umin) umin = u2
+            // Step 8:
+        } while (u > EXPONENTIAL_SA_QI[i]) // Ensured to exit since EXPONENTIAL_SA_QI[MAX] = 1.
+
+        mean * (a + umin * EXPONENTIAL_SA_QI[0])
+    }
+
+    override fun toString(): String = "Ahrens-Dieter Exponential deviate"
+
+    public companion object {
+        private val EXPONENTIAL_SA_QI by lazy { DoubleArray(16) }
+
+        init {
+            /**
+             * Filling EXPONENTIAL_SA_QI table.
+             * Note that we don't want qi = 0 in the table.
+             */
+            val ln2 = ln(2.0)
+            var qi = 0.0
+
+            EXPONENTIAL_SA_QI.indices.forEach { i ->
+                qi += ln2.pow(i + 1.0) / InternalUtils.factorial(i + 1)
+                EXPONENTIAL_SA_QI[i] = qi
+            }
+        }
+
+        public fun of(mean: Double): AhrensDieterExponentialSampler {
+            require(mean > 0) { "mean is not strictly positive: $mean" }
+            return AhrensDieterExponentialSampler(mean)
+        }
+    }
+}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AhrensDieterMarsagliaTsangGammaSampler.kt

@@ -0,0 +1,120 @@
+package space.kscience.kmath.stat.samplers
+
+import space.kscience.kmath.chains.Chain
+import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.stat.Sampler
+import space.kscience.kmath.stat.chain
+import space.kscience.kmath.stat.next
+import kotlin.math.*
+
+/**
+ * Sampling from the [gamma distribution](http://mathworld.wolfram.com/GammaDistribution.html).
+ *  - For 0 < alpha < 1:
+ *    Ahrens, J. H. and Dieter, U., Computer methods for sampling from gamma, beta, Poisson and binomial distributions, Computing, 12, 223-246, 1974.
+ *  - For alpha >= 1:
+ *    Marsaglia and Tsang, A Simple Method for Generating Gamma Variables. ACM Transactions on Mathematical Software, Volume 26 Issue 3, September, 2000.
+ *
+ * Based on Commons RNG implementation.
+ *
+ * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/AhrensDieterMarsagliaTsangGammaSampler.html].
+ */
+public class AhrensDieterMarsagliaTsangGammaSampler private constructor(
+    alpha: Double,
+    theta: Double
+) : Sampler<Double> {
+    private val delegate: BaseGammaSampler =
+        if (alpha < 1) AhrensDieterGammaSampler(alpha, theta) else MarsagliaTsangGammaSampler(alpha, theta)
+
+    private abstract class BaseGammaSampler internal constructor(
+        protected val alpha: Double,
+        protected val theta: Double
+    ) : Sampler<Double> {
+        init {
+            require(alpha > 0) { "alpha is not strictly positive: $alpha" }
+            require(theta > 0) { "theta is not strictly positive: $theta" }
+        }
+
+        override fun toString(): String = "Ahrens-Dieter-Marsaglia-Tsang Gamma deviate"
+    }
+
+    private class AhrensDieterGammaSampler(alpha: Double, theta: Double) :
+        BaseGammaSampler(alpha, theta) {
+        private val oneOverAlpha: Double = 1.0 / alpha
+        private val bGSOptim: Double = 1.0 + alpha / E
+
+        override fun sample(generator: RandomGenerator): Chain<Double> = generator.chain {
+            var x: Double
+
+            // [1]: p. 228, Algorithm GS.
+            while (true) {
+                // Step 1:
+                val u = generator.nextDouble()
+                val p = bGSOptim * u
+
+                if (p <= 1) {
+                    // Step 2:
+                    x = p.pow(oneOverAlpha)
+                    val u2 = generator.nextDouble()
+
+                    if (u2 > exp(-x)) // Reject.
+                        continue
+
+                    break
+                }
+
+                // Step 3:
+                x = -ln((bGSOptim - p) * oneOverAlpha)
+                val u2: Double = generator.nextDouble()
+                if (u2 <= x.pow(alpha - 1.0)) break
+                // Reject and continue.
+            }
+
+            x * theta
+        }
+    }
+
+    private class MarsagliaTsangGammaSampler(alpha: Double, theta: Double) :
+        BaseGammaSampler(alpha, theta) {
+        private val dOptim: Double
+        private val cOptim: Double
+        private val gaussian: NormalizedGaussianSampler
+
+        init {
+            gaussian = ZigguratNormalizedGaussianSampler.of()
+            dOptim = alpha - ONE_THIRD
+            cOptim = ONE_THIRD / sqrt(dOptim)
+        }
+
+        override fun sample(generator: RandomGenerator): Chain<Double> = generator.chain {
+            var v: Double
+
+            while (true) {
+                val x = gaussian.next(generator)
+                val oPcTx = 1 + cOptim * x
+                v = oPcTx * oPcTx * oPcTx
+                if (v <= 0) continue
+                val x2 = x * x
+                val u = generator.nextDouble()
+                // Squeeze.
+                if (u < 1 - 0.0331 * x2 * x2) break
+                if (ln(u) < 0.5 * x2 + dOptim * (1 - v + ln(v))) break
+            }
+
+            theta * dOptim * v
+        }
+
+        companion object {
+            private const val ONE_THIRD = 1.0 / 3.0
+        }
+    }
+
+    public override fun sample(generator: RandomGenerator): Chain<Double> = delegate.sample(generator)
+    public override fun toString(): String = delegate.toString()
+
+    public companion object {
+        public fun of(
+            alpha: Double,
+            theta: Double
+        ): Sampler<Double> = AhrensDieterMarsagliaTsangGammaSampler(alpha, theta)
+    }
+}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/AliasMethodDiscreteSampler.kt

@@ -0,0 +1,286 @@
+package space.kscience.kmath.stat.samplers
+
+import space.kscience.kmath.chains.Chain
+import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.stat.Sampler
+import space.kscience.kmath.stat.chain
+import space.kscience.kmath.stat.internal.InternalUtils
+import kotlin.math.ceil
+import kotlin.math.max
+import kotlin.math.min
+
+/**
+ * Distribution sampler that uses the Alias method. It can be used to sample from n values each with an associated
+ * probability. This implementation is based on the detailed explanation of the alias method by Keith Schartz and
+ * implements Vose's algorithm.
+ *
+ * Vose, M.D., A linear algorithm for generating random numbers with a given distribution, IEEE Transactions on
+ * Software Engineering, 17, 972-975, 1991. he algorithm will sample values in O(1) time after a pre-processing step
+ * of O(n) time.
+ *
+ * The alias tables are constructed using fraction probabilities with an assumed denominator of 253. In the generic
+ * case sampling uses UniformRandomProvider.nextInt(int) and the upper 53-bits from UniformRandomProvider.nextLong().
+ *
+ * Zero padding the input probabilities can be used to make more sampling more efficient. Any zero entry will always be
+ * aliased removing the requirement to compute a long. Increased sampling speed comes at the cost of increased storage
+ * space. The algorithm requires approximately 12 bytes of storage per input probability, that is n * 12 for size n.
+ * Zero-padding only requires 4 bytes of storage per padded value as the probability is known to be zero.
+ *
+ * An optimisation is performed for small table sizes that are a power of 2. In this case the sampling uses 1 or 2
+ * calls from UniformRandomProvider.nextInt() to generate up to 64-bits for creation of an 11-bit index and 53-bits
+ * for the long. This optimisation requires a generator with a high cycle length for the lower order bits.
+ *
+ * Larger table sizes that are a power of 2 will benefit from fast algorithms for UniformRandomProvider.nextInt(int)
+ * that exploit the power of 2.
+ *
+ * Based on Commons RNG implementation.
+ * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/AliasMethodDiscreteSampler.html].
+ */
+public open class AliasMethodDiscreteSampler private constructor(
+    // Deliberate direct storage of input arrays
+    protected val probability: LongArray,
+    protected val alias: IntArray
+) : Sampler<Int> {
+
+    private class SmallTableAliasMethodDiscreteSampler(
+        probability: LongArray,
+        alias: IntArray
+    ) : AliasMethodDiscreteSampler(probability, alias) {
+        // Assume the table size is a power of 2 and create the mask
+        private val mask: Int = alias.size - 1
+
+        override fun sample(generator: RandomGenerator): Chain<Int> = generator.chain {
+            val bits = generator.nextInt()
+            // Isolate lower bits
+            val j = bits and mask
+
+            // Optimisation for zero-padded input tables
+            if (j >= probability.size)
+            // No probability must use the alias
+                return@chain alias[j]
+
+            // Create a uniform random deviate as a long.
+            // This replicates functionality from the o.a.c.rng.core.utils.NumberFactory.makeLong
+            val longBits = generator.nextInt().toLong() shl 32 or (bits.toLong() and hex_ffffffff)
+            // Choose between the two. Use a 53-bit long for the probability.
+            if (longBits ushr 11 < probability[j]) j else alias[j]
+        }
+
+        private companion object {
+            private const val hex_ffffffff = 4294967295L
+        }
+    }
+
+    public override fun sample(generator: RandomGenerator): Chain<Int> = generator.chain {
+        // This implements the algorithm as per Vose (1991):
+        // v = uniform()  in [0, 1)
+        // j = uniform(n) in [0, n)
+        // if v < prob[j] then
+        //   return j
+        // else
+        //   return alias[j]
+        val j = generator.nextInt(alias.size)
+
+        // Optimisation for zero-padded input tables
+        // No probability must use the alias
+        if (j >= probability.size) return@chain alias[j]
+
+        // Note: We could check the probability before computing a deviate.
+        // p(j) == 0  => alias[j]
+        // p(j) == 1  => j
+        // However it is assumed these edge cases are rare:
+        //
+        // The probability table will be 1 for approximately 1/n samples, i.e. only the
+        // last unpaired probability. This is only worth checking for when the table size (n)
+        // is small. But in that case the user should zero-pad the table for performance.
+        //
+        // The probability table will be 0 when an input probability was zero. We
+        // will assume this is also rare if modelling a discrete distribution where
+        // all samples are possible. The edge case for zero-padded tables is handled above.
+
+        // Choose between the two. Use a 53-bit long for the probability.
+        if (generator.nextLong() ushr 11 < probability[j]) j else alias[j]
+    }
+
+    public override fun toString(): String = "Alias method"
+
+    public companion object {
+        private const val DEFAULT_ALPHA = 0
+        private const val ZERO = 0.0
+        private const val ONE_AS_NUMERATOR = 1L shl 53
+        private const val CONVERT_TO_NUMERATOR: Double = ONE_AS_NUMERATOR.toDouble()
+        private const val MAX_SMALL_POWER_2_SIZE = 1 shl 11
+
+        public fun of(
+            probabilities: DoubleArray,
+            alpha: Int = DEFAULT_ALPHA
+        ): Sampler<Int> {
+            // The Alias method balances N categories with counts around the mean into N sections,
+            // each allocated 'mean' observations.
+            //
+            // Consider 4 categories with counts 6,3,2,1. The histogram can be balanced into a
+            // 2D array as 4 sections with a height of the mean:
+            //
+            // 6
+            // 6
+            // 6
+            // 63   => 6366   --
+            // 632     6326    |-- mean
+            // 6321    6321   --
+            //
+            // section abcd
+            //
+            // Each section is divided as:
+            // a: 6=1/1
+            // b: 3=1/1
+            // c: 2=2/3; 6=1/3   (6 is the alias)
+            // d: 1=1/3; 6=2/3   (6 is the alias)
+            //
+            // The sample is obtained by randomly selecting a section, then choosing which category
+            // from the pair based on a uniform random deviate.
+            val sumProb = InternalUtils.validateProbabilities(probabilities)
+            // Allow zero-padding
+            val n = computeSize(probabilities.size, alpha)
+            // Partition into small and large by splitting on the average.
+            val mean = sumProb / n
+            // The cardinality of smallSize + largeSize = n.
+            // So fill the same array from either end.
+            val indices = IntArray(n)
+            var large = n
+            var small = 0
+
+            probabilities.indices.forEach { i ->
+                if (probabilities[i] >= mean) indices[--large] = i else indices[small++] = i
+            }
+
+            small = fillRemainingIndices(probabilities.size, indices, small)
+            // This may be smaller than the input length if the probabilities were already padded.
+            val nonZeroIndex = findLastNonZeroIndex(probabilities)
+            // The probabilities are modified so use a copy.
+            // Note: probabilities are required only up to last nonZeroIndex
+            val remainingProbabilities = probabilities.copyOf(nonZeroIndex + 1)
+            // Allocate the final tables.
+            // Probability table may be truncated (when zero padded).
+            // The alias table is full length.
+            val probability = LongArray(remainingProbabilities.size)
+            val alias = IntArray(n)
+
+            // This loop uses each large in turn to fill the alias table for small probabilities that
+            // do not reach the requirement to fill an entire section alone (i.e. p < mean).
+            // Since the sum of the small should be less than the sum of the large it should use up
+            // all the small first. However floating point round-off can result in
+            // misclassification of items as small or large. The Vose algorithm handles this using
+            // a while loop conditioned on the size of both sets and a subsequent loop to use
+            // unpaired items.
+            while (large != n && small != 0) {
+                // Index of the small and the large probabilities.
+                val j = indices[--small]
+                val k = indices[large++]
+
+                // Optimisation for zero-padded input:
+                // p(j) = 0 above the last nonZeroIndex
+                if (j > nonZeroIndex)
+                // The entire amount for the section is taken from the alias.
+                    remainingProbabilities[k] -= mean
+                else {
+                    val pj = remainingProbabilities[j]
+                    // Item j is a small probability that is below the mean.
+                    // Compute the weight of the section for item j: pj / mean.
+                    // This is scaled by 2^53 and the ceiling function used to round-up
+                    // the probability to a numerator of a fraction in the range [1,2^53].
+                    // Ceiling ensures non-zero values.
+                    probability[j] = ceil(CONVERT_TO_NUMERATOR * (pj / mean)).toLong()
+                    // The remaining amount for the section is taken from the alias.
+                    // Effectively: probabilities[k] -= (mean - pj)
+                    remainingProbabilities[k] += pj - mean
+                }
+
+                // If not j then the alias is k
+                alias[j] = k
+
+                // Add the remaining probability from large to the appropriate list.
+                if (remainingProbabilities[k] >= mean) indices[--large] = k else indices[small++] = k
+            }
+
+            // Final loop conditions to consume unpaired items.
+            // Note: The large set should never be non-empty but this can occur due to round-off
+            // error so consume from both.
+            fillTable(probability, alias, indices, 0, small)
+            fillTable(probability, alias, indices, large, n)
+
+            // Change the algorithm for small power of 2 sized tables
+            return if (isSmallPowerOf2(n))
+                SmallTableAliasMethodDiscreteSampler(probability, alias)
+            else
+                AliasMethodDiscreteSampler(probability, alias)
+        }
+
+        private fun fillRemainingIndices(length: Int, indices: IntArray, small: Int): Int {
+            var updatedSmall = small
+            (length until indices.size).forEach { i -> indices[updatedSmall++] = i }
+            return updatedSmall
+        }
+
+        private fun findLastNonZeroIndex(probabilities: DoubleArray): Int {
+            // No bounds check is performed when decrementing as the array contains at least one
+            // value above zero.
+            var nonZeroIndex = probabilities.size - 1
+            while (probabilities[nonZeroIndex] == ZERO) nonZeroIndex--
+            return nonZeroIndex
+        }
+
+        private fun computeSize(length: Int, alpha: Int): Int {
+            // If No padding
+            if (alpha < 0) return length
+            // Use the number of leading zeros function to find the next power of 2,
+            // i.e. ceil(log2(x))
+            var pow2 = 32 - numberOfLeadingZeros(length - 1)
+            // Increase by the alpha. Clip this to limit to a positive integer (2^30)
+            pow2 = min(30, pow2 + alpha)
+            // Use max to handle a length above the highest possible power of 2
+            return max(length, 1 shl pow2)
+        }
+
+        private fun fillTable(
+            probability: LongArray,
+            alias: IntArray,
+            indices: IntArray,
+            start: Int,
+            end: Int
+        ) = (start until end).forEach { i ->
+            val index = indices[i]
+            probability[index] = ONE_AS_NUMERATOR
+            alias[index] = index
+        }
+
+        private fun isSmallPowerOf2(n: Int): Boolean = n <= MAX_SMALL_POWER_2_SIZE && n and n - 1 == 0
+
+        private fun numberOfLeadingZeros(i: Int): Int {
+            var mutI = i
+            if (mutI <= 0) return if (mutI == 0) 32 else 0
+            var n = 31
+
+            if (mutI >= 1 shl 16) {
+                n -= 16
+                mutI = mutI ushr 16
+            }
+
+            if (mutI >= 1 shl 8) {
+                n -= 8
+                mutI = mutI ushr 8
+            }
+
+            if (mutI >= 1 shl 4) {
+                n -= 4
+                mutI = mutI ushr 4
+            }
+
+            if (mutI >= 1 shl 2) {
+                n -= 2
+                mutI = mutI ushr 2
+            }
+
+            return n - (mutI ushr 1)
+        }
+    }
+}

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/stat/samplers/BoxMullerNormalizedGaussianSampler.kt

@@ -0,0 +1,48 @@
+package space.kscience.kmath.stat.samplers
+
+import space.kscience.kmath.chains.Chain
+import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.stat.Sampler
+import space.kscience.kmath.stat.chain
+import kotlin.math.*
+
+/**
+ * [Box-Muller algorithm](https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform) for sampling from a Gaussian
+ * distribution.
+ *
+ * Based on Commons RNG implementation.
+ * See [https://commons.apache.org/proper/commons-rng/commons-rng-sampling/apidocs/org/apache/commons/rng/sampling/distribution/BoxMullerNormalizedGaussianSampler.html].
+ */
+public class BoxMullerNormalizedGaussianSampler private constructor() : NormalizedGaussianSampler, Sampler<Double> {
+    private var nextGaussian: Double = Double.NaN
+
+    public override fun sample(generator: RandomGenerator): Chain<Double> = generator.chain {
+        val random: Double
+
+        if (nextGaussian.isNaN()) {
+            // Generate a pair of Gaussian numbers.
+            val x = nextDouble()
+            val y = nextDouble()
+            val alpha = 2 * PI * x
+            val r = sqrt(-2 * ln(y))
+            // Return the first element of the generated pair.
+            random = r * cos(alpha)
+            // Keep second element of the pair for next invocation.
+            nextGaussian = r * sin(alpha)
+        } else {
+            // Use the second element of the pair (generated at the
+            // previous invocation).
+            random = nextGaussian
+            // Both elements of the pair have been used.
+            nextGaussian = Double.NaN
+        }
+
+        random
+    }
+
+    public override fun toString(): String = "Box-Muller normalized Gaussian deviate"
+
+    public companion object {
+        public fun of(): BoxMullerNormalizedGaussianSampler = BoxMullerNormalizedGaussianSampler()
+    }
+}