benchmarks/build.gradle

@@ -1,39 +0,0 @@
-plugins {
-    id "java"
-    id "me.champeau.gradle.jmh" version "0.4.8"
-    id 'org.jetbrains.kotlin.jvm'
-}
-
-repositories {
-    maven { url 'https://dl.bintray.com/kotlin/kotlin-eap' }
-    maven{ url "http://dl.bintray.com/kyonifer/maven"}
-    mavenCentral()
-}
-
-dependencies {
-    implementation project(":kmath-core")
-    implementation project(":kmath-coroutines")
-    implementation project(":kmath-commons")
-    implementation project(":kmath-koma")
-    implementation group: "com.kyonifer", name:"koma-core-ejml", version: "0.12"
-    implementation "org.jetbrains.kotlinx:kotlinx-io-jvm:0.1.5"
-    //compile "org.jetbrains.kotlin:kotlin-stdlib-jdk8"
-    //jmh project(':kmath-core')
-}
-
-jmh {
-    warmupIterations = 1
-}
-
-jmhClasses.dependsOn(compileKotlin)
-
-compileKotlin {
-    kotlinOptions {
-        jvmTarget = "1.8"
-    }
-}
-compileTestKotlin {
-    kotlinOptions {
-        jvmTarget = "1.8"
-    }
-}

build.gradle.kts

@@ -1,4 +1,4 @@
-val kmathVersion by extra("0.1.2")
+val kmathVersion by extra("0.1.3-dev-1")
 
 allprojects {
     repositories {

doc/nd-structure.md

@@ -9,17 +9,12 @@ structures. In `kmath` performance depends on which particular context was used
 
 Let us consider following contexts:
 ```kotlin
-    // specialized nd-field for Double. It works as generic Double field as well
-    val specializedField = NDField.real(intArrayOf(dim, dim))
-
     // automatically build context most suited for given type.
-    val autoField = NDField.auto(intArrayOf(dim, dim), RealField)
-
-    //A field implementing lazy computations. All elements are computed on-demand
-    val lazyField = NDField.lazy(intArrayOf(dim, dim), RealField)
-
+    val autoField = NDField.auto(RealField, dim, dim)
+    // specialized nd-field for Double. It works as generic Double field as well
+    val specializedField = NDField.real(dim, dim)
     //A generic boxing field. It should be used for objects, not primitives.
-    val genericField = NDField.buffered(intArrayOf(dim, dim), RealField)
+    val genericField = NDField.buffered(RealField, dim, dim)
 ```
 Now let us perform several tests and see which implementation is best suited for each case:
 
@@ -32,7 +27,7 @@ to it `n = 1000` times.
 The code to run this looks like:
 ```kotlin
     specializedField.run {
-        var res = one
+        var res: NDBuffer<Double> = one
         repeat(n) {
             res += 1.0
         }
@@ -93,7 +88,7 @@ In this case it completes in about `4x-5x` time due to boxing.
 The boxing field produced by
 ```kotlin
     genericField.run {
-        var res = one
+        var res: NDBuffer<Double> = one
         repeat(n) {
             res += 1.0
         }

examples/build.gradle.kts

@@ -0,0 +1,67 @@
+import org.jetbrains.gradle.benchmarks.JvmBenchmarkTarget
+import org.jetbrains.kotlin.allopen.gradle.AllOpenExtension
+import org.jetbrains.kotlin.gradle.tasks.KotlinCompile
+
+plugins {
+    java
+    kotlin("jvm")
+    kotlin("plugin.allopen") version "1.3.31"
+    id("org.jetbrains.gradle.benchmarks.plugin") version "0.1.7-dev-24"
+}
+
+configure<AllOpenExtension> {
+    annotation("org.openjdk.jmh.annotations.State")
+}
+
+repositories {
+    maven("https://dl.bintray.com/kotlin/kotlin-eap")
+    maven("http://dl.bintray.com/kyonifer/maven")
+    maven("https://dl.bintray.com/orangy/maven")
+    mavenCentral()
+}
+
+sourceSets {
+    register("benchmarks")
+}
+
+dependencies {
+    implementation(project(":kmath-core"))
+    implementation(project(":kmath-coroutines"))
+    implementation(project(":kmath-commons"))
+    implementation(project(":kmath-koma"))
+    implementation("com.kyonifer:koma-core-ejml:0.12")
+    implementation("org.jetbrains.kotlinx:kotlinx-io-jvm:0.1.5")
+
+    implementation("org.jetbrains.gradle.benchmarks:runtime:0.1.7-dev-24")
+
+
+    "benchmarksCompile"(sourceSets.main.get().compileClasspath)
+}
+
+// Configure benchmark
+benchmark {
+    // Setup configurations
+    targets {
+        // This one matches sourceSet name above
+        register("benchmarks") {
+            this as JvmBenchmarkTarget
+            jmhVersion = "1.21"
+        }
+    }
+
+    configurations {
+        register("fast") {
+            warmups = 5 // number of warmup iterations
+            iterations = 3 // number of iterations
+            iterationTime = 500 // time in seconds per iteration
+            iterationTimeUnit = "ms" // time unity for iterationTime, default is seconds
+        }
+    }
+}
+
+
+tasks.withType<KotlinCompile> {
+    kotlinOptions {
+        jvmTarget = "1.8"
+    }
+}

examples/src/benchmarks/kotlin/scientifik/kmath/structures/ArrayBenchmark.kt

@@ -7,7 +7,7 @@ import java.nio.IntBuffer
 
 
 @State(Scope.Benchmark)
-open class ArrayBenchmark {
+class ArrayBenchmark {
 
     @Benchmark
     fun benchmarkArrayRead() {

examples/src/benchmarks/kotlin/scientifik/kmath/structures/BufferBenchmark.kt

@@ -7,7 +7,7 @@ import scientifik.kmath.operations.Complex
 import scientifik.kmath.operations.complex
 
 @State(Scope.Benchmark)
-open class BufferBenchmark {
+class BufferBenchmark {
 
     @Benchmark
     fun genericDoubleBufferReadWrite() {

examples/src/benchmarks/kotlin/scientifik/kmath/structures/NDFieldBenchmark.kt

@@ -1,9 +1,12 @@
 package scientifik.kmath.structures
 
 import org.openjdk.jmh.annotations.Benchmark
+import org.openjdk.jmh.annotations.Scope
+import org.openjdk.jmh.annotations.State
 import scientifik.kmath.operations.RealField
 
-open class NDFieldBenchmark {
+@State(Scope.Benchmark)
+class NDFieldBenchmark {
 
     @Benchmark
     fun autoFieldAdd() {
@@ -50,6 +53,6 @@ open class NDFieldBenchmark {
 
         val bufferedField = NDField.auto(RealField, dim, dim)
         val specializedField = NDField.real(dim, dim)
-        val genericField = NDField.buffered(intArrayOf(dim, dim), RealField)
+        val genericField = NDField.buffered(RealField, dim, dim)
     }
 }

examples/src/main/kotlin/scientifik/kmath/commons/prob/DistributionDemo.kt

@@ -0,0 +1,31 @@
+package scientifik.kmath.commons.prob
+
+import kotlinx.coroutines.runBlocking
+import scientifik.kmath.chains.Chain
+import scientifik.kmath.chains.StatefulChain
+import scientifik.kmath.prob.Distribution
+import scientifik.kmath.prob.RandomGenerator
+
+data class AveragingChainState(var num: Int = 0, var value: Double = 0.0)
+
+fun Chain<Double>.mean(): Chain<Double> = StatefulChain(AveragingChainState(), 0.0) {
+    val next = this@mean.next()
+    num++
+    value += next
+    return@StatefulChain value / num
+}
+
+
+fun main() {
+    val normal = Distribution.normal()
+    val chain = normal.sample(RandomGenerator.default).mean()
+
+    runBlocking {
+        repeat(10001) { counter ->
+            val mean = chain.next()
+            if(counter % 1000 ==0){
+                println("[$counter] Average value is $mean")
+            }
+        }
+    }
+}

examples/src/main/kotlin/scientifik/kmath/linear/LinearAlgebraBenchmark.kt

@@ -1,6 +1,9 @@
 package scientifik.kmath.linear
 
 import koma.matrix.ejml.EJMLMatrixFactory
+import scientifik.kmath.commons.linear.CMMatrixContext
+import scientifik.kmath.commons.linear.inverse
+import scientifik.kmath.commons.linear.toCM
 import scientifik.kmath.operations.RealField
 import scientifik.kmath.structures.Matrix
 import kotlin.contracts.ExperimentalContracts

examples/src/main/kotlin/scientifik/kmath/linear/MultiplicationBenchmark.kt

@@ -1,6 +1,8 @@
 package scientifik.kmath.linear
 
 import koma.matrix.ejml.EJMLMatrixFactory
+import scientifik.kmath.commons.linear.CMMatrixContext
+import scientifik.kmath.commons.linear.toCM
 import scientifik.kmath.operations.RealField
 import scientifik.kmath.structures.Matrix
 import kotlin.random.Random

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

@@ -0,0 +1,10 @@
+package scientifik.kmath.operations
+
+import scientifik.kmath.structures.NDElement
+import scientifik.kmath.structures.complex
+
+fun main() {
+    val element = NDElement.complex(2, 2) { index: IntArray ->
+        Complex(index[0].toDouble() - index[1].toDouble(), index[0].toDouble() + index[1].toDouble())
+    }
+}

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

@@ -13,7 +13,7 @@ fun main(args: Array<String>) {
     // specialized nd-field for Double. It works as generic Double field as well
     val specializedField = NDField.real(dim, dim)
     //A generic boxing field. It should be used for objects, not primitives.
-    val genericField = NDField.buffered(intArrayOf(dim, dim), RealField)
+    val genericField = NDField.buffered(RealField, dim, dim)
 
 
     val autoTime = measureTimeMillis {

kmath-commons/build.gradle.kts

@@ -8,6 +8,7 @@ description = "Commons math binding for kmath"
 dependencies {
     api(project(":kmath-core"))
     api(project(":kmath-coroutines"))
+    api(project(":kmath-prob"))
     api("org.apache.commons:commons-math3:3.6.1")
     testImplementation("org.jetbrains.kotlin:kotlin-test")
     testImplementation("org.jetbrains.kotlin:kotlin-test-junit")

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

@@ -1,6 +1,8 @@
-package scientifik.kmath.expressions
+package scientifik.kmath.commons.expressions
 
 import org.apache.commons.math3.analysis.differentiation.DerivativeStructure
+import scientifik.kmath.expressions.Expression
+import scientifik.kmath.expressions.ExpressionContext
 import scientifik.kmath.operations.ExtendedField
 import scientifik.kmath.operations.Field
 import kotlin.properties.ReadOnlyProperty
@@ -81,8 +83,12 @@ class DerivativeStructureField(
 /**
  * A constructs that creates a derivative structure with required order on-demand
  */
-class DiffExpression(val function: DerivativeStructureField.() -> DerivativeStructure) : Expression<Double> {
-    override fun invoke(arguments: Map<String, Double>): Double = DerivativeStructureField(0, arguments)
+class DiffExpression(val function: DerivativeStructureField.() -> DerivativeStructure) :
+    Expression<Double> {
+    override fun invoke(arguments: Map<String, Double>): Double = DerivativeStructureField(
+        0,
+        arguments
+    )
         .run(function).value
 
     /**
@@ -109,21 +115,27 @@ fun DiffExpression.derivative(name: String) = derivative(name to 1)
  * A context for [DiffExpression] (not to be confused with [DerivativeStructure])
  */
 object DiffExpressionContext : ExpressionContext<Double>, Field<DiffExpression> {
-    override fun variable(name: String, default: Double?) = DiffExpression { variable(name, default?.const()) }
+    override fun variable(name: String, default: Double?) =
+        DiffExpression { variable(name, default?.const()) }
 
-    override fun const(value: Double): DiffExpression = DiffExpression { value.const() }
+    override fun const(value: Double): DiffExpression =
+        DiffExpression { value.const() }
 
-    override fun add(a: DiffExpression, b: DiffExpression) = DiffExpression { a.function(this) + b.function(this) }
+    override fun add(a: DiffExpression, b: DiffExpression) =
+        DiffExpression { a.function(this) + b.function(this) }
 
     override val zero = DiffExpression { 0.0.const() }
 
-    override fun multiply(a: DiffExpression, k: Number) = DiffExpression { a.function(this) * k }
+    override fun multiply(a: DiffExpression, k: Number) =
+        DiffExpression { a.function(this) * k }
 
     override val one = DiffExpression { 1.0.const() }
 
-    override fun multiply(a: DiffExpression, b: DiffExpression) = DiffExpression { a.function(this) * b.function(this) }
+    override fun multiply(a: DiffExpression, b: DiffExpression) =
+        DiffExpression { a.function(this) * b.function(this) }
 
-    override fun divide(a: DiffExpression, b: DiffExpression) = DiffExpression { a.function(this) / b.function(this) }
+    override fun divide(a: DiffExpression, b: DiffExpression) =
+        DiffExpression { a.function(this) / b.function(this) }
 }
 
 

kmath-commons/src/main/kotlin/scientifik/kmath/commons/linear/CMMatrix.kt

@@ -1,11 +1,13 @@
-package scientifik.kmath.linear
+package scientifik.kmath.commons.linear
 
 import org.apache.commons.math3.linear.*
 import org.apache.commons.math3.linear.RealMatrix
 import org.apache.commons.math3.linear.RealVector
+import scientifik.kmath.linear.*
 import scientifik.kmath.structures.Matrix
 
-class CMMatrix(val origin: RealMatrix, features: Set<MatrixFeature>? = null) : FeaturedMatrix<Double> {
+class CMMatrix(val origin: RealMatrix, features: Set<MatrixFeature>? = null) :
+    FeaturedMatrix<Double> {
     override val rowNum: Int get() = origin.rowDimension
     override val colNum: Int get() = origin.columnDimension
 
@@ -70,10 +72,14 @@ object CMMatrixContext : MatrixContext<Double> {
     override fun multiply(a: Matrix<Double>, k: Number) =
         CMMatrix(a.toCM().origin.scalarMultiply(k.toDouble()))
 
-    override fun Matrix<Double>.times(value: Double): Matrix<Double>  = produce(rowNum,colNum){i,j-> get(i,j)*value}
+    override fun Matrix<Double>.times(value: Double): Matrix<Double>  =
+        produce(rowNum, colNum) { i, j -> get(i, j) * value }
 }
 
-operator fun CMMatrix.plus(other: CMMatrix): CMMatrix = CMMatrix(this.origin.add(other.origin))
-operator fun CMMatrix.minus(other: CMMatrix): CMMatrix = CMMatrix(this.origin.subtract(other.origin))
+operator fun CMMatrix.plus(other: CMMatrix): CMMatrix =
+    CMMatrix(this.origin.add(other.origin))
+operator fun CMMatrix.minus(other: CMMatrix): CMMatrix =
+    CMMatrix(this.origin.subtract(other.origin))
 
-infix fun CMMatrix.dot(other: CMMatrix): CMMatrix = CMMatrix(this.origin.multiply(other.origin))
+infix fun CMMatrix.dot(other: CMMatrix): CMMatrix =
+    CMMatrix(this.origin.multiply(other.origin))

kmath-commons/src/main/kotlin/scientifik/kmath/commons/linear/CMSolver.kt

@@ -1,6 +1,7 @@
-package scientifik.kmath.linear
+package scientifik.kmath.commons.linear
 
 import org.apache.commons.math3.linear.*
+import scientifik.kmath.linear.Point
 import scientifik.kmath.structures.Matrix
 
 enum class CMDecomposition {

kmath-commons/src/main/kotlin/scientifik/kmath/commons/prob/CMRandomGeneratorWrapper.kt

@@ -0,0 +1,28 @@
+package scientifik.kmath.commons.prob
+
+import org.apache.commons.math3.random.JDKRandomGenerator
+import scientifik.kmath.prob.RandomGenerator
+import org.apache.commons.math3.random.RandomGenerator as CMRandom
+
+inline class CMRandomGeneratorWrapper(val generator: CMRandom) : RandomGenerator {
+    override fun nextDouble(): Double = generator.nextDouble()
+
+    override fun nextInt(): Int = generator.nextInt()
+
+    override fun nextLong(): Long = generator.nextLong()
+
+    override fun nextBlock(size: Int): ByteArray = ByteArray(size).apply { generator.nextBytes(this) }
+}
+
+fun CMRandom.asKmathGenerator(): RandomGenerator = CMRandomGeneratorWrapper(this)
+
+fun RandomGenerator.asCMGenerator(): CMRandom =
+    (this as? CMRandomGeneratorWrapper)?.generator ?: TODO("Implement reverse CM wrapper")
+
+val RandomGenerator.Companion.default: RandomGenerator by lazy { JDKRandomGenerator().asKmathGenerator() }
+
+fun RandomGenerator.Companion.jdk(seed: Int? = null): RandomGenerator = if (seed == null) {
+    JDKRandomGenerator()
+} else {
+    JDKRandomGenerator(seed)
+}.asKmathGenerator()

kmath-commons/src/main/kotlin/scientifik/kmath/commons/prob/CommonsDistribution.kt

@@ -0,0 +1,82 @@
+package scientifik.kmath.commons.prob
+
+import org.apache.commons.math3.distribution.*
+import scientifik.kmath.prob.Distribution
+import scientifik.kmath.prob.RandomChain
+import scientifik.kmath.prob.RandomGenerator
+import scientifik.kmath.prob.UnivariateDistribution
+import org.apache.commons.math3.random.RandomGenerator as CMRandom
+
+class CMRealDistributionWrapper(val builder: (CMRandom?) -> RealDistribution) : UnivariateDistribution<Double> {
+
+    private val defaultDistribution by lazy { builder(null) }
+
+    override fun probability(arg: Double): Double = defaultDistribution.probability(arg)
+
+    override fun cumulative(arg: Double): Double = defaultDistribution.cumulativeProbability(arg)
+
+    override fun sample(generator: RandomGenerator): RandomChain<Double> {
+        val distribution = builder(generator.asCMGenerator())
+        return RandomChain(generator) { distribution.sample() }
+    }
+}
+
+class CMIntDistributionWrapper(val builder: (CMRandom?) -> IntegerDistribution) : UnivariateDistribution<Int> {
+
+    private val defaultDistribution by lazy { builder(null) }
+
+    override fun probability(arg: Int): Double = defaultDistribution.probability(arg)
+
+    override fun cumulative(arg: Int): Double = defaultDistribution.cumulativeProbability(arg)
+
+    override fun sample(generator: RandomGenerator): RandomChain<Int> {
+        val distribution = builder(generator.asCMGenerator())
+        return RandomChain(generator) { distribution.sample() }
+    }
+}
+
+
+fun Distribution.Companion.normal(mean: Double = 0.0, sigma: Double = 1.0): UnivariateDistribution<Double> =
+    CMRealDistributionWrapper { generator -> NormalDistribution(generator, mean, sigma) }
+
+fun Distribution.Companion.poisson(mean: Double): UnivariateDistribution<Int> = CMIntDistributionWrapper { generator ->
+    PoissonDistribution(
+        generator,
+        mean,
+        PoissonDistribution.DEFAULT_EPSILON,
+        PoissonDistribution.DEFAULT_MAX_ITERATIONS
+    )
+}
+
+fun Distribution.Companion.binomial(trials: Int, p: Double): UnivariateDistribution<Int> =
+    CMIntDistributionWrapper { generator ->
+        BinomialDistribution(generator, trials, p)
+    }
+
+fun Distribution.Companion.student(degreesOfFreedom: Double): UnivariateDistribution<Double> =
+    CMRealDistributionWrapper { generator ->
+        TDistribution(generator, degreesOfFreedom, TDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY)
+    }
+
+fun Distribution.Companion.chi2(degreesOfFreedom: Double): UnivariateDistribution<Double> =
+    CMRealDistributionWrapper { generator ->
+        ChiSquaredDistribution(generator, degreesOfFreedom)
+    }
+
+fun Distribution.Companion.fisher(
+    numeratorDegreesOfFreedom: Double,
+    denominatorDegreesOfFreedom: Double
+): UnivariateDistribution<Double> =
+    CMRealDistributionWrapper { generator ->
+        FDistribution(generator, numeratorDegreesOfFreedom, denominatorDegreesOfFreedom)
+    }
+
+fun Distribution.Companion.exponential(mean: Double): UnivariateDistribution<Double> =
+    CMRealDistributionWrapper { generator ->
+        ExponentialDistribution(generator, mean)
+    }
+
+fun Distribution.Companion.uniform(a: Double, b: Double): UnivariateDistribution<Double> =
+    CMRealDistributionWrapper { generator ->
+        UniformRealDistribution(generator, a, b)
+    }

kmath-commons/src/main/kotlin/scientifik/kmath/commons/transform/Transformations.kt

@@ -1,4 +1,4 @@
-package scientifik.kmath.transform
+package scientifik.kmath.commons.transform
 
 import kotlinx.coroutines.FlowPreview
 import kotlinx.coroutines.flow.Flow

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

@@ -1,6 +1,9 @@
-package scientifik.kmath.expressions
+package scientifik.kmath.commons.expressions
 
 import org.junit.Test
+import scientifik.kmath.commons.expressions.DerivativeStructureField
+import scientifik.kmath.commons.expressions.DiffExpression
+import scientifik.kmath.commons.expressions.derivative
 import kotlin.test.assertEquals
 
 inline fun <R> diff(order: Int, vararg parameters: Pair<String, Double>, block: DerivativeStructureField.() -> R) =

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

@@ -131,4 +131,7 @@ operator fun ComplexNDElement.plus(arg: Double) =
 operator fun ComplexNDElement.minus(arg: Double) =
     map { it - arg }
 
-fun NDField.Companion.complex(vararg shape: Int) = ComplexNDField(shape)
+fun NDField.Companion.complex(vararg shape: Int): ComplexNDField = ComplexNDField(shape)
+
+fun NDElement.Companion.complex(vararg shape: Int, initializer: ComplexField.(IntArray) -> Complex): ComplexNDElement =
+    NDField.complex(*shape).produce(initializer)

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

@@ -1,9 +1,9 @@
 package scientifik.kmath.structures
 
+import scientifik.kmath.operations.Complex
 import scientifik.kmath.operations.Field
 import scientifik.kmath.operations.Ring
 import scientifik.kmath.operations.Space
-import kotlin.jvm.JvmName
 
 
 /**
@@ -57,6 +57,8 @@ interface NDAlgebra<T, C, N : NDStructure<T>> {
      * element-by-element invoke a function working on [T] on a [NDStructure]
      */
     operator fun Function1<T, T>.invoke(structure: N) = map(structure) { value -> this@invoke(value) }
+
+    companion object
 }
 
 /**
@@ -75,6 +77,7 @@ interface NDSpace<T, S : Space<T>, N : NDStructure<T>> : Space<N>, NDAlgebra<T,
 
     //TODO move to extensions after KEEP-176
     operator fun N.plus(arg: T) = map(this) { value -> add(arg, value) }
+
     operator fun N.minus(arg: T) = map(this) { value -> add(arg, -value) }
 
     operator fun T.plus(arg: N) = map(arg) { value -> add(this@plus, value) }
@@ -93,6 +96,7 @@ interface NDRing<T, R : Ring<T>, N : NDStructure<T>> : Ring<N>, NDSpace<T, R, N>
 
     //TODO move to extensions after KEEP-176
     operator fun N.times(arg: T) = map(this) { value -> multiply(arg, value) }
+
     operator fun T.times(arg: N) = map(arg) { value -> multiply(this@times, value) }
 }
 
@@ -113,6 +117,7 @@ interface NDField<T, F : Field<T>, N : NDStructure<T>> : Field<N>, NDRing<T, F,
 
     //TODO move to extensions after KEEP-176
     operator fun N.div(arg: T) = map(this) { value -> divide(arg, value) }
+
     operator fun T.div(arg: N) = map(arg) { divide(it, this@div) }
 
     companion object {
@@ -128,11 +133,10 @@ interface NDField<T, F : Field<T>, N : NDStructure<T>> : Field<N>, NDRing<T, F,
          * Create a nd-field with boxing generic buffer
          */
         fun <T : Any, F : Field<T>> buffered(
-            shape: IntArray,
             field: F,
+            vararg shape: Int,
             bufferFactory: BufferFactory<T> = Buffer.Companion::boxing
-        ) =
-            BoxingNDField(shape, field, bufferFactory)
+        ) = BoxingNDField(shape, field, bufferFactory)
 
         /**
          * Create a most suitable implementation for nd-field using reified class.
@@ -141,6 +145,7 @@ interface NDField<T, F : Field<T>, N : NDStructure<T>> : Field<N>, NDRing<T, F,
         inline fun <reified T : Any, F : Field<T>> auto(field: F, vararg shape: Int): BufferedNDField<T, F> =
             when {
                 T::class == Double::class -> real(*shape) as BufferedNDField<T, F>
+                T::class == Complex::class -> complex(*shape) as BufferedNDField<T, F>
                 else -> BoxingNDField(shape, field, Buffer.Companion::auto)
             }
     }

kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/chains/Chain.kt

@@ -18,6 +18,7 @@ package scientifik.kmath.chains
 
 import kotlinx.atomicfu.atomic
 import kotlinx.coroutines.FlowPreview
+import kotlinx.coroutines.flow.Flow
 
 
 /**
@@ -28,7 +29,7 @@ interface Chain<out R> {
     /**
      * Last cached value of the chain. Returns null if [next] was not called
      */
-    val value: R?
+    fun peek(): R?
 
     /**
      * Generate next value, changing state if needed
@@ -46,13 +47,12 @@ interface Chain<out R> {
  * Chain as a coroutine flow. The flow emit affects chain state and vice versa
  */
 @FlowPreview
-val <R> Chain<R>.flow
+val <R> Chain<R>.flow: Flow<R>
     get() = kotlinx.coroutines.flow.flow { while (true) emit(next()) }
 
 fun <T> Iterator<T>.asChain(): Chain<T> = SimpleChain { next() }
 fun <T> Sequence<T>.asChain(): Chain<T> = iterator().asChain()
 
-
 /**
  * Map the chain result using suspended transformation. Initial chain result can no longer be safely consumed
  * since mapped chain consumes tokens. Accepts regular transformation function
@@ -60,7 +60,7 @@ fun <T> Sequence<T>.asChain(): Chain<T> = iterator().asChain()
 fun <T, R> Chain<T>.map(func: (T) -> R): Chain<R> {
     val parent = this;
     return object : Chain<R> {
-        override val value: R? get() = parent.value?.let(func)
+        override fun peek(): R? = parent.peek()?.let(func)
 
         override suspend fun next(): R {
             return func(parent.next())
@@ -77,7 +77,7 @@ fun <T, R> Chain<T>.map(func: (T) -> R): Chain<R> {
  */
 class SimpleChain<out R>(private val gen: suspend () -> R) : Chain<R> {
     private val atomicValue = atomic<R?>(null)
-    override val value: R? get() = atomicValue.value
+    override fun peek(): R? = atomicValue.value
 
     override suspend fun next(): R = gen().also { atomicValue.lazySet(it) }
 
@@ -95,16 +95,16 @@ class MarkovChain<out R : Any>(private val seed: () -> R, private val gen: suspe
     constructor(seed: R, gen: suspend (R) -> R) : this({ seed }, gen)
 
     private val atomicValue = atomic<R?>(null)
-    override val value: R get() = atomicValue.value ?: seed()
+    override fun peek(): R = atomicValue.value ?: seed()
 
     override suspend fun next(): R {
-        val newValue = gen(value)
+        val newValue = gen(peek())
         atomicValue.lazySet(newValue)
-        return value
+        return peek()
     }
 
     override fun fork(): Chain<R> {
-        return MarkovChain(value, gen)
+        return MarkovChain(peek(), gen)
     }
 }
 
@@ -121,12 +121,12 @@ class StatefulChain<S, out R>(
     constructor(state: S, seed: R, gen: suspend S.(R) -> R) : this(state, { seed }, gen)
 
     private val atomicValue = atomic<R?>(null)
-    override val value: R get() = atomicValue.value ?: seed(state)
+    override fun peek(): R = atomicValue.value ?: seed(state)
 
     override suspend fun next(): R {
-        val newValue = gen(state, value)
+        val newValue = gen(state, peek())
         atomicValue.lazySet(newValue)
-        return value
+        return peek()
     }
 
     override fun fork(): Chain<R> {
@@ -137,10 +137,10 @@ class StatefulChain<S, out R>(
 /**
  * A chain that repeats the same value
  */
-class ConstantChain<out T>(override val value: T) : Chain<T> {
-    override suspend fun next(): T {
-        return value
-    }
+class ConstantChain<out T>(val value: T) : Chain<T> {
+    override fun peek(): T? = value
+
+    override suspend fun next(): T = value
 
     override fun fork(): Chain<T> {
         return this

kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/coroutines/CoroutinesExtra.kt

@@ -1,4 +1,4 @@
-package scientifik.kmath
+package scientifik.kmath.coroutines
 
 import kotlinx.coroutines.*
 import kotlinx.coroutines.channels.produce
@@ -42,7 +42,8 @@ fun <T, R> Flow<T>.async(
 }
 
 @FlowPreview
-fun <T, R> AsyncFlow<T>.map(action: (T) -> R) = AsyncFlow(deferredFlow.map { input ->
+fun <T, R> AsyncFlow<T>.map(action: (T) -> R) =
+    AsyncFlow(deferredFlow.map { input ->
         //TODO add function composition
         LazyDeferred(input.dispatcher) {
             input.start(this)

kmath-coroutines/src/commonMain/kotlin/scientifik/kmath/streaming/BufferFlow.kt

@@ -22,7 +22,7 @@ fun <T> Flow<Buffer<out T>>.spread(): Flow<T> = flatMapConcat { it.asFlow() }
  * Collect incoming flow into fixed size chunks
  */
 @FlowPreview
-fun <T> Flow<T>.chunked(bufferSize: Int, bufferFactory: BufferFactory<T>) = flow {
+fun <T> Flow<T>.chunked(bufferSize: Int, bufferFactory: BufferFactory<T>): Flow<Buffer<T>> = flow {
     require(bufferSize > 0) { "Resulting chunk size must be more than zero" }
     val list = ArrayList<T>(bufferSize)
     var counter = 0
@@ -46,7 +46,7 @@ fun <T> Flow<T>.chunked(bufferSize: Int, bufferFactory: BufferFactory<T>) = flow
  * Specialized flow chunker for real buffer
  */
 @FlowPreview
-fun Flow<Double>.chunked(bufferSize: Int) = flow {
+fun Flow<Double>.chunked(bufferSize: Int): Flow<DoubleBuffer> = flow {
     require(bufferSize > 0) { "Resulting chunk size must be more than zero" }
     val array = DoubleArray(bufferSize)
     var counter = 0

kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/chains/ChainExt.kt

@@ -27,7 +27,7 @@ fun <R> Chain<R>.asSequence(): Sequence<R> = object : Sequence<R> {
 fun <T, R> Chain<T>.map(func: suspend (T) -> R): Chain<R> {
     val parent = this;
     return object : Chain<R> {
-        override val value: R? get() = runBlocking { parent.value?.let { func(it) } }
+        override fun peek(): R? = runBlocking { parent.peek()?.let { func(it) } }
 
         override suspend fun next(): R {
             return func(parent.next())

kmath-coroutines/src/jvmMain/kotlin/scientifik/kmath/structures/LazyNDStructure.kt

@@ -1,7 +1,7 @@
 package scientifik.kmath.structures
 
 import kotlinx.coroutines.*
-import scientifik.kmath.Math
+import scientifik.kmath.coroutines.Math
 
 class LazyNDStructure<T>(
     val scope: CoroutineScope,

kmath-coroutines/src/jvmTest/kotlin/scientifik/kmath/streaming/BufferFlowTest.kt

@@ -4,9 +4,9 @@ import kotlinx.coroutines.*
 import kotlinx.coroutines.flow.asFlow
 import kotlinx.coroutines.flow.collect
 import org.junit.Test
-import scientifik.kmath.async
-import scientifik.kmath.collect
-import scientifik.kmath.map
+import scientifik.kmath.coroutines.async
+import scientifik.kmath.coroutines.collect
+import scientifik.kmath.coroutines.map
 import java.util.concurrent.Executors
 
 

kmath-prob/build.gradle.kts

@@ -6,13 +6,14 @@ plugins {
 kotlin.sourceSets {
     commonMain {
         dependencies {
-            api(project(":kmath-core"))
             api(project(":kmath-coroutines"))
             compileOnly("org.jetbrains.kotlinx:atomicfu-common:${Versions.atomicfuVersion}")
         }
     }
     jvmMain {
         dependencies {
+            // https://mvnrepository.com/artifact/org.apache.commons/commons-rng-simple
+            //api("org.apache.commons:commons-rng-sampling:1.2")
             compileOnly("org.jetbrains.kotlinx:atomicfu:${Versions.atomicfuVersion}")
         }
     }

kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/Distribution.kt

@@ -0,0 +1,33 @@
+package scientifik.kmath.prob
+
+import scientifik.kmath.chains.Chain
+
+/**
+ * A distribution of typed objects
+ */
+interface Distribution<T : Any> {
+    /**
+     * A probability value for given argument [arg].
+     * For continuous distributions returns PDF
+     */
+    fun probability(arg: T): Double
+
+    /**
+     * Create a chain of samples from this distribution.
+     * The chain is not guaranteed to be stateless.
+     */
+    fun sample(generator: RandomGenerator): RandomChain<T>
+    //TODO add sample bunch generator
+
+    /**
+     * An empty companion. Distribution factories should be written as its extensions
+     */
+    companion object
+}
+
+interface UnivariateDistribution<T : Comparable<T>> : Distribution<T> {
+    /**
+     * Cumulative distribution for ordered parameter
+     */
+    fun cumulative(arg: T): Double
+}

kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomChain.kt

@@ -0,0 +1,17 @@
+package scientifik.kmath.prob
+
+import kotlinx.atomicfu.atomic
+import scientifik.kmath.chains.Chain
+
+/**
+ * A possibly stateful chain producing random values.
+ * TODO make random chain properly fork generator
+ */
+class RandomChain<out R>(val generator: RandomGenerator, private val gen: suspend RandomGenerator.() -> R) : Chain<R> {
+    private val atomicValue = atomic<R?>(null)
+    override fun peek(): R? = atomicValue.value
+
+    override suspend fun next(): R = generator.gen().also { atomicValue.lazySet(it) }
+
+    override fun fork(): Chain<R> = RandomChain(generator, gen)
+}

kmath-prob/src/commonMain/kotlin/scientifik/kmath/prob/RandomGenerator.kt

@@ -0,0 +1,13 @@
+package scientifik.kmath.prob
+
+/**
+ * A basic generator
+ */
+interface RandomGenerator {
+    fun nextDouble(): Double
+    fun nextInt(): Int
+    fun nextLong(): Long
+    fun nextBlock(size: Int): ByteArray
+
+    companion object
+}

settings.gradle.kts

@@ -3,6 +3,7 @@ pluginManagement {
         jcenter()
         gradlePluginPortal()
         maven("https://dl.bintray.com/kotlin/kotlin-eap")
+        maven("https://dl.bintray.com/orangy/maven")
     }
     resolutionStrategy {
         eachPlugin {
@@ -28,5 +29,5 @@ include(
     ":kmath-commons",
     ":kmath-koma",
     ":kmath-prob",
-    ":benchmarks"
+    ":examples"
 )