Merge remote-tracking branch 'mipt-npm/dev' into foreign-memory

build.gradle.kts

@@ -1,8 +1,8 @@
 plugins {
-    id("scientifik.publish") version "0.4.2" apply false
+    id("scientifik.publish") apply false
 }
 
-val kmathVersion by extra("0.1.4-dev-4")
+val kmathVersion by extra("0.1.4-dev-7")
 
 val bintrayRepo by extra("scientifik")
 val githubProject by extra("kmath")

docs/codestyle.md

@@ -0,0 +1,22 @@
+# Local coding conventions
+
+Kmath and other `scientifik` projects use general [kotlin code conventions](https://kotlinlang.org/docs/reference/coding-conventions.html), but with a number of small changes and clarifications.
+
+## Utility class names
+File name should coincide with a name of one of the classes contained in the file or start with small letter and describe its contents.
+
+The code convention [here](https://kotlinlang.org/docs/reference/coding-conventions.html#source-file-names) says that file names should start with capital letter even if file does not contain classes. Yet starting utility classes and aggregators with a small letter seems to be a good way to clearly visually separate those files.
+
+This convention could be changed in future in a non-breaking way.
+
+## Private variable names
+Private variable names could start with underscore `_` in case the private mutable variable is shadowed by the public read-only value with the same meaning.
+
+Code convention do not permit underscores in names, but is is sometimes useful to "underscore" the fact that public and private versions define the same entity. It is allowed only for private variables.
+
+This convention could be changed in future in a non-breaking way.
+
+## Functions and properties one-liners
+Use one-liners when they occupy single code window line both for functions and properties with getters like `val b: String get() = "fff"`. The same should be done with multiline expressions when they could be cleanly separated.
+
+There is not general consensus whenever use `fun a() = {}` or `fun a(){return}`. Yet from reader perspective one-lines seem to better show that the property or function is easily calculated.

examples/build.gradle.kts

@@ -27,6 +27,7 @@ dependencies {
     implementation(project(":kmath-core"))
     implementation(project(":kmath-coroutines"))
     implementation(project(":kmath-commons"))
+    implementation(project(":kmath-prob"))
     implementation(project(":kmath-koma"))
     implementation(project(":kmath-viktor"))
     implementation(project(":kmath-dimensions"))
@@ -57,6 +58,6 @@ benchmark {
 
 tasks.withType<KotlinCompile> {
     kotlinOptions {
-        jvmTarget = Scientifik.JVM_VERSION
+        jvmTarget = Scientifik.JVM_TARGET.toString()
     }
 }

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

@@ -0,0 +1,71 @@
+package scientifik.kmath.commons.prob
+
+import kotlinx.coroutines.Dispatchers
+import kotlinx.coroutines.async
+import kotlinx.coroutines.runBlocking
+import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler
+import org.apache.commons.rng.simple.RandomSource
+import scientifik.kmath.chains.BlockingRealChain
+import scientifik.kmath.prob.*
+import java.time.Duration
+import java.time.Instant
+
+
+private suspend fun runChain(): Duration {
+    val generator = RandomGenerator.fromSource(RandomSource.MT, 123L)
+
+    val normal = Distribution.normal(NormalSamplerMethod.Ziggurat)
+    val chain = normal.sample(generator) as BlockingRealChain
+
+    val startTime = Instant.now()
+    var sum = 0.0
+    repeat(10000001) { counter ->
+
+        sum += chain.nextDouble()
+
+        if (counter % 100000 == 0) {
+            val duration = Duration.between(startTime, Instant.now())
+            val meanValue = sum / counter
+            println("Chain sampler completed $counter elements in $duration: $meanValue")
+        }
+    }
+    return Duration.between(startTime, Instant.now())
+}
+
+private fun runDirect(): Duration {
+    val provider = RandomSource.create(RandomSource.MT, 123L)
+    val sampler = ZigguratNormalizedGaussianSampler(provider)
+    val startTime = Instant.now()
+
+    var sum = 0.0
+    repeat(10000001) { counter ->
+
+        sum += sampler.sample()
+
+        if (counter % 100000 == 0) {
+            val duration = Duration.between(startTime, Instant.now())
+            val meanValue = sum / counter
+            println("Direct sampler completed $counter elements in $duration: $meanValue")
+        }
+    }
+    return Duration.between(startTime, Instant.now())
+}
+
+/**
+ * Comparing chain sampling performance with direct sampling performance
+ */
+fun main() {
+    runBlocking(Dispatchers.Default) {
+        val chainJob = async {
+            runChain()
+        }
+
+        val directJob = async {
+            runDirect()
+        }
+
+        println("Chain: ${chainJob.await()}")
+        println("Direct: ${directJob.await()}")
+    }
+
+}

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

@@ -5,10 +5,11 @@ import scientifik.kmath.chains.Chain
 import scientifik.kmath.chains.collectWithState
 import scientifik.kmath.prob.Distribution
 import scientifik.kmath.prob.RandomGenerator
+import scientifik.kmath.prob.normal
 
 data class AveragingChainState(var num: Int = 0, var value: Double = 0.0)
 
-fun Chain<Double>.mean(): Chain<Double> = collectWithState(AveragingChainState(),{it.copy()}){ chain->
+fun Chain<Double>.mean(): Chain<Double> = collectWithState(AveragingChainState(), { it.copy() }) { chain ->
     val next = chain.next()
     num++
     value += next

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

@@ -5,6 +5,7 @@ import org.apache.commons.math3.linear.RealMatrix
 import org.apache.commons.math3.linear.RealVector
 import scientifik.kmath.linear.*
 import scientifik.kmath.structures.Matrix
+import scientifik.kmath.structures.NDStructure
 
 class CMMatrix(val origin: RealMatrix, features: Set<MatrixFeature>? = null) :
     FeaturedMatrix<Double> {
@@ -19,6 +20,16 @@ class CMMatrix(val origin: RealMatrix, features: Set<MatrixFeature>? = null) :
         CMMatrix(origin, this.features + features)
 
     override fun get(i: Int, j: Int): Double = origin.getEntry(i, j)
+
+    override fun equals(other: Any?): Boolean {
+        return NDStructure.equals(this, other as? NDStructure<*> ?: return false)
+    }
+
+    override fun hashCode(): Int {
+        var result = origin.hashCode()
+        result = 31 * result + features.hashCode()
+        return result
+    }
 }
 
 fun Matrix<Double>.toCM(): CMMatrix = if (this is CMMatrix) {

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

@@ -1,32 +0,0 @@
-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) }
-
-    override fun fork(): RandomGenerator {
-        TODO("not implemented") //To change body of created functions use File | Settings | File Templates.
-    }
-}
-
-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

@@ -1,82 +0,0 @@
-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-core/src/commonMain/kotlin/scientifik/kmath/linear/BufferMatrix.kt

@@ -27,7 +27,7 @@ class BufferMatrixContext<T : Any, R : Ring<T>>(
 @Suppress("OVERRIDE_BY_INLINE")
 object RealMatrixContext : GenericMatrixContext<Double, RealField> {
 
-    override val elementContext = RealField
+    override val elementContext get() = RealField
 
     override inline fun produce(rows: Int, columns: Int, initializer: (i: Int, j: Int) -> Double): Matrix<Double> {
         val buffer = DoubleBuffer(rows * columns) { offset -> initializer(offset / columns, offset % columns) }
@@ -101,8 +101,15 @@ infix fun BufferMatrix<Double>.dot(other: BufferMatrix<Double>): BufferMatrix<Do
 
     val array = DoubleArray(this.rowNum * other.colNum)
 
-    val a = this.buffer.array
+    //convert to array to insure there is not memory indirection
-    val b = other.buffer.array
+    fun Buffer<out Double>.unsafeArray(): DoubleArray = if (this is DoubleBuffer) {
+        array
+    } else {
+        DoubleArray(size) { get(it) }
+    }
+
+    val a = this.buffer.unsafeArray()
+    val b = other.buffer.unsafeArray()
 
     for (i in (0 until rowNum)) {
         for (j in (0 until other.colNum)) {

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

@@ -1,12 +1,8 @@
 package scientifik.kmath.linear
 
-import scientifik.kmath.operations.Field
-import scientifik.kmath.operations.Norm
-import scientifik.kmath.operations.RealField
 import scientifik.kmath.structures.Buffer
 import scientifik.kmath.structures.Matrix
 import scientifik.kmath.structures.VirtualBuffer
-import scientifik.kmath.structures.asSequence
 
 typealias Point<T> = Buffer<T>
 
@@ -19,8 +15,6 @@ interface LinearSolver<T : Any> {
     fun inverse(a: Matrix<T>): Matrix<T>
 }
 
-typealias RealMatrix = Matrix<Double>
-
 /**
  * Convert matrix to vector if it is possible
  */

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

@@ -1,14 +1,46 @@
 package scientifik.kmath.linear
 
+import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.BufferFactory
 import scientifik.kmath.structures.Structure2D
 import scientifik.kmath.structures.asBuffer
 
-class MatrixBuilder<T : Any>(val rows: Int, val columns: Int) {
+class MatrixBuilder(val rows: Int, val columns: Int) {
-    operator fun invoke(vararg elements: T): FeaturedMatrix<T> {
+    operator fun <T : Any> invoke(vararg elements: T): FeaturedMatrix<T> {
         if (rows * columns != elements.size) error("The number of elements ${elements.size} is not equal $rows * $columns")
         val buffer = elements.asBuffer()
         return BufferMatrix(rows, columns, buffer)
     }
+
+    //TODO add specific matrix builder functions like diagonal, etc
 }
 
-fun <T : Any> Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder<T> = MatrixBuilder(rows, columns)
+fun Structure2D.Companion.build(rows: Int, columns: Int): MatrixBuilder = MatrixBuilder(rows, columns)
+
+fun <T : Any> Structure2D.Companion.row(vararg values: T): FeaturedMatrix<T> {
+    val buffer = values.asBuffer()
+    return BufferMatrix(1, values.size, buffer)
+}
+
+inline fun <reified T : Any> Structure2D.Companion.row(
+    size: Int,
+    factory: BufferFactory<T> = Buffer.Companion::auto,
+    noinline builder: (Int) -> T
+): FeaturedMatrix<T> {
+    val buffer = factory(size, builder)
+    return BufferMatrix(1, size, buffer)
+}
+
+fun <T : Any> Structure2D.Companion.column(vararg values: T): FeaturedMatrix<T> {
+    val buffer = values.asBuffer()
+    return BufferMatrix(values.size, 1, buffer)
+}
+
+inline fun <reified T : Any> Structure2D.Companion.column(
+    size: Int,
+    factory: BufferFactory<T> = Buffer.Companion::auto,
+    noinline builder: (Int) -> T
+): FeaturedMatrix<T> {
+    val buffer = factory(size, builder)
+    return BufferMatrix(size, 1, buffer)
+}

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

@@ -2,6 +2,7 @@ package scientifik.kmath.operations
 
 import scientifik.kmath.operations.BigInt.Companion.BASE
 import scientifik.kmath.operations.BigInt.Companion.BASE_SIZE
+import scientifik.kmath.structures.*
 import kotlin.math.log2
 import kotlin.math.max
 import kotlin.math.min
@@ -482,3 +483,18 @@ fun String.parseBigInteger(): BigInt? {
     }
     return res * sign
 }
+
+inline fun Buffer.Companion.bigInt(size: Int, initializer: (Int) -> BigInt): Buffer<BigInt> =
+    boxing(size, initializer)
+
+inline fun MutableBuffer.Companion.bigInt(size: Int, initializer: (Int) -> BigInt): MutableBuffer<BigInt> =
+    boxing(size, initializer)
+
+fun NDAlgebra.Companion.bigInt(vararg shape: Int): BoxingNDRing<BigInt, BigIntField> =
+    BoxingNDRing(shape, BigIntField, Buffer.Companion::bigInt)
+
+fun NDElement.Companion.bigInt(
+    vararg shape: Int,
+    initializer: BigIntField.(IntArray) -> BigInt
+): BufferedNDRingElement<BigInt, BigIntField> =
+    NDAlgebra.bigInt(*shape).produce(initializer)

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

@@ -102,7 +102,7 @@ object IntRing : Ring<Int>, Norm<Int, Int> {
     override val zero: Int = 0
     override inline fun add(a: Int, b: Int) = a + b
     override inline fun multiply(a: Int, b: Int) = a * b
-    override inline fun multiply(a: Int, k: Number) = (k * a)
+    override inline fun multiply(a: Int, k: Number) = k.toInt() * a
     override val one: Int = 1
 
     override inline fun norm(arg: Int) = abs(arg)
@@ -124,7 +124,7 @@ object ShortRing : Ring<Short>, Norm<Short, Short> {
     override val zero: Short = 0
     override inline fun add(a: Short, b: Short) = (a + b).toShort()
     override inline fun multiply(a: Short, b: Short) = (a * b).toShort()
-    override inline fun multiply(a: Short, k: Number) = (a * k)
+    override inline fun multiply(a: Short, k: Number) = (a * k.toShort()).toShort()
     override val one: Short = 1
 
     override fun norm(arg: Short): Short = if (arg > 0) arg else (-arg).toShort()
@@ -146,7 +146,7 @@ object ByteRing : Ring<Byte>, Norm<Byte, Byte> {
     override val zero: Byte = 0
     override inline fun add(a: Byte, b: Byte) = (a + b).toByte()
     override inline fun multiply(a: Byte, b: Byte) = (a * b).toByte()
-    override inline fun multiply(a: Byte, k: Number) = (a * k)
+    override inline fun multiply(a: Byte, k: Number) = (a * k.toByte()).toByte()
     override val one: Byte = 1
 
     override fun norm(arg: Byte): Byte = if (arg > 0) arg else (-arg).toByte()
@@ -168,7 +168,7 @@ object LongRing : Ring<Long>, Norm<Long, Long> {
     override val zero: Long = 0
     override inline fun add(a: Long, b: Long) = (a + b)
     override inline fun multiply(a: Long, b: Long) = (a * b)
-    override inline fun multiply(a: Long, k: Number) = (a * k)
+    override inline fun multiply(a: Long, k: Number) = a * k.toLong()
     override val one: Long = 1
 
     override fun norm(arg: Long): Long = abs(arg)

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

@@ -3,10 +3,10 @@ package scientifik.kmath.structures
 import scientifik.kmath.operations.*
 
 /**
- * Base interface for an element with context, containing strides
+ * Base class for an element with context, containing strides
  */
-interface BufferedNDElement<T, C> : NDBuffer<T>, NDElement<T, C, NDBuffer<T>> {
+abstract class BufferedNDElement<T, C> : NDBuffer<T>(), NDElement<T, C, NDBuffer<T>> {
-    override val context: BufferedNDAlgebra<T, C>
+    abstract override val context: BufferedNDAlgebra<T, C>
 
     override val strides get() = context.strides
 
@@ -16,7 +16,7 @@ interface BufferedNDElement<T, C> : NDBuffer<T>, NDElement<T, C, NDBuffer<T>> {
 class BufferedNDSpaceElement<T, S : Space<T>>(
     override val context: BufferedNDSpace<T, S>,
     override val buffer: Buffer<T>
-) : BufferedNDElement<T, S>, SpaceElement<NDBuffer<T>, BufferedNDSpaceElement<T, S>, BufferedNDSpace<T, S>> {
+) : BufferedNDElement<T, S>(), SpaceElement<NDBuffer<T>, BufferedNDSpaceElement<T, S>, BufferedNDSpace<T, S>> {
 
     override fun unwrap(): NDBuffer<T> = this
 
@@ -29,7 +29,7 @@ class BufferedNDSpaceElement<T, S : Space<T>>(
 class BufferedNDRingElement<T, R : Ring<T>>(
     override val context: BufferedNDRing<T, R>,
     override val buffer: Buffer<T>
-) : BufferedNDElement<T, R>, RingElement<NDBuffer<T>, BufferedNDRingElement<T, R>, BufferedNDRing<T, R>> {
+) : BufferedNDElement<T, R>(), RingElement<NDBuffer<T>, BufferedNDRingElement<T, R>, BufferedNDRing<T, R>> {
 
     override fun unwrap(): NDBuffer<T> = this
 
@@ -42,7 +42,7 @@ class BufferedNDRingElement<T, R : Ring<T>>(
 class BufferedNDFieldElement<T, F : Field<T>>(
     override val context: BufferedNDField<T, F>,
     override val buffer: Buffer<T>
-) : BufferedNDElement<T, F>, FieldElement<NDBuffer<T>, BufferedNDFieldElement<T, F>, BufferedNDField<T, F>> {
+) : BufferedNDElement<T, F>(), FieldElement<NDBuffer<T>, BufferedNDFieldElement<T, F>, BufferedNDField<T, F>> {
 
     override fun unwrap(): NDBuffer<T> = this
 

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

@@ -71,7 +71,7 @@ interface Buffer<T> {
 
 fun <T> Buffer<T>.asSequence(): Sequence<T> = Sequence(::iterator)
 
-fun <T> Buffer<T>.asIterable(): Iterable<T> = asSequence().asIterable()
+fun <T> Buffer<T>.asIterable(): Iterable<T> = Iterable(::iterator)
 
 val Buffer<*>.indices: IntRange get() = IntRange(0, size - 1)
 
@@ -161,36 +161,7 @@ class ArrayBuffer<T>(private val array: Array<T>) : MutableBuffer<T> {
     override fun copy(): MutableBuffer<T> = ArrayBuffer(array.copyOf())
 }
 
-fun <T> Array<T>.asBuffer() = ArrayBuffer(this)
+fun <T> Array<T>.asBuffer(): ArrayBuffer<T> = ArrayBuffer(this)
-
-inline class DoubleBuffer(val array: DoubleArray) : MutableBuffer<Double> {
-    override val size: Int get() = array.size
-
-    override fun get(index: Int): Double = array[index]
-
-    override fun set(index: Int, value: Double) {
-        array[index] = value
-    }
-
-    override fun iterator() = array.iterator()
-
-    override fun copy(): MutableBuffer<Double> = DoubleBuffer(array.copyOf())
-}
-
-@Suppress("FunctionName")
-inline fun DoubleBuffer(size: Int, init: (Int) -> Double) = DoubleBuffer(DoubleArray(size) { init(it) })
-
-/**
- * Transform buffer of doubles into array for high performance operations
- */
-val Buffer<out Double>.array: DoubleArray
-    get() = if (this is DoubleBuffer) {
-        array
-    } else {
-        DoubleArray(size) { get(it) }
-    }
-
-fun DoubleArray.asBuffer() = DoubleBuffer(this)
 
 inline class ShortBuffer(val array: ShortArray) : MutableBuffer<Short> {
     override val size: Int get() = array.size

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

@@ -0,0 +1,34 @@
+package scientifik.kmath.structures
+
+inline class DoubleBuffer(val array: DoubleArray) : MutableBuffer<Double> {
+    override val size: Int get() = array.size
+
+    override fun get(index: Int): Double = array[index]
+
+    override fun set(index: Int, value: Double) {
+        array[index] = value
+    }
+
+    override fun iterator() = array.iterator()
+
+    override fun copy(): MutableBuffer<Double> =
+        DoubleBuffer(array.copyOf())
+}
+
+@Suppress("FunctionName")
+inline fun DoubleBuffer(size: Int, init: (Int) -> Double): DoubleBuffer = DoubleBuffer(DoubleArray(size) { init(it) })
+
+@Suppress("FunctionName")
+fun DoubleBuffer(vararg doubles: Double): DoubleBuffer = DoubleBuffer(doubles)
+
+/**
+ * Transform buffer of doubles into array for high performance operations
+ */
+val MutableBuffer<out Double>.array: DoubleArray
+    get() = if (this is DoubleBuffer) {
+        array
+    } else {
+        DoubleArray(size) { get(it) }
+    }
+
+fun DoubleArray.asBuffer() = DoubleBuffer(this)

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

@@ -14,15 +14,25 @@ interface NDStructure<T> {
 
     fun elements(): Sequence<Pair<IntArray, T>>
 
+    override fun equals(other: Any?): Boolean
+
+    override fun hashCode(): Int
+
     companion object {
         fun equals(st1: NDStructure<*>, st2: NDStructure<*>): Boolean {
-            return when {
+            if (st1 === st2) return true
-                st1 === st2 -> true
+
-                st1 is BufferNDStructure<*> && st2 is BufferNDStructure<*> && st1.strides == st2.strides -> st1.buffer.contentEquals(
+            // fast comparison of buffers if possible
-                    st2.buffer
+            if (
-                )
+                st1 is NDBuffer &&
-                else -> st1.elements().all { (index, value) -> value == st2[index] }
+                st2 is NDBuffer &&
+                st1.strides == st2.strides
+            ) {
+                return st1.buffer.contentEquals(st2.buffer)
             }
+
+            //element by element comparison if it could not be avoided
+            return st1.elements().all { (index, value) -> value == st2[index] }
         }
 
         /**
@@ -177,15 +187,25 @@ class DefaultStrides private constructor(override val shape: IntArray) : Strides
     }
 }
 
-interface NDBuffer<T> : NDStructure<T> {
+abstract class NDBuffer<T> : NDStructure<T> {
-    val buffer: Buffer<T>
+    abstract val buffer: Buffer<T>
-    val strides: Strides
+    abstract val strides: Strides
 
     override fun get(index: IntArray): T = buffer[strides.offset(index)]
 
     override val shape: IntArray get() = strides.shape
 
-    override fun elements() = strides.indices().map { it to this[it] }
+    override fun elements(): Sequence<Pair<IntArray, T>> = strides.indices().map { it to this[it] }
+
+    override fun equals(other: Any?): Boolean {
+        return NDStructure.equals(this, other as? NDStructure<*> ?: return false)
+    }
+
+    override fun hashCode(): Int {
+        var result = strides.hashCode()
+        result = 31 * result + buffer.hashCode()
+        return result
+    }
 }
 
 /**
@@ -194,34 +214,12 @@ interface NDBuffer<T> : NDStructure<T> {
 class BufferNDStructure<T>(
     override val strides: Strides,
     override val buffer: Buffer<T>
-) : NDBuffer<T> {
+) : NDBuffer<T>() {
-
     init {
         if (strides.linearSize != buffer.size) {
             error("Expected buffer side of ${strides.linearSize}, but found ${buffer.size}")
         }
     }
-
-    override fun get(index: IntArray): T = buffer[strides.offset(index)]
-
-    override val shape: IntArray get() = strides.shape
-
-    override fun elements() = strides.indices().map { it to this[it] }
-
-    override fun equals(other: Any?): Boolean {
-        return when {
-            this === other -> true
-            other is BufferNDStructure<*> && this.strides == other.strides -> this.buffer.contentEquals(other.buffer)
-            other is NDStructure<*> -> elements().all { (index, value) -> value == other[index] }
-            else -> false
-        }
-    }
-
-    override fun hashCode(): Int {
-        var result = strides.hashCode()
-        result = 31 * result + buffer.hashCode()
-        return result
-    }
 }
 
 /**
@@ -245,7 +243,7 @@ inline fun <T, reified R : Any> NDStructure<T>.mapToBuffer(
 class MutableBufferNDStructure<T>(
     override val strides: Strides,
     override val buffer: MutableBuffer<T>
-) : NDBuffer<T>, MutableNDStructure<T> {
+) : NDBuffer<T>(), MutableNDStructure<T> {
 
     init {
         if (strides.linearSize != buffer.size) {

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

@@ -17,7 +17,7 @@ interface Structure1D<T> : NDStructure<T>, Buffer<T> {
 /**
  * A 1D wrapper for nd-structure
  */
-private inline class Structure1DWrapper<T>(val structure: NDStructure<T>) : Structure1D<T> {
+private inline class Structure1DWrapper<T>(val structure: NDStructure<T>) : Structure1D<T>{
 
     override val shape: IntArray get() = structure.shape
     override val size: Int get() = structure.shape[0]

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

@@ -14,7 +14,6 @@ interface Structure2D<T> : NDStructure<T> {
         return get(index[0], index[1])
     }
 
-
     val rows: Buffer<Buffer<T>>
         get() = VirtualBuffer(rowNum) { i ->
             VirtualBuffer(colNum) { j -> get(i, j) }
@@ -58,22 +57,4 @@ fun <T> NDStructure<T>.as2D(): Structure2D<T> = if (shape.size == 2) {
     error("Can't create 2d-structure from ${shape.size}d-structure")
 }
 
-/**
- * Represent this 2D structure as 1D if it has exactly one column. Throw error otherwise.
- */
-fun <T> Structure2D<T>.as1D() = if (colNum == 1) {
-    object : Structure1D<T> {
-        override fun get(index: Int): T = get(index, 0)
-
-        override val shape: IntArray get() = intArrayOf(rowNum)
-
-        override fun elements(): Sequence<Pair<IntArray, T>> = elements()
-
-        override val size: Int get() = rowNum
-    }
-} else {
-    error("Can't convert matrix with more than one column to vector")
-}
-
-
 typealias Matrix<T> = Structure2D<T>

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

@@ -17,7 +17,7 @@ class MatrixTest {
 
     @Test
     fun testBuilder() {
-        val matrix = Matrix.build<Double>(2, 3)(
+        val matrix = Matrix.build(2, 3)(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0
         )

kmath-core/src/jvmMain/kotlin/scientifik/kmath/operations/bigNumbers.kt

@@ -5,7 +5,7 @@ import java.math.BigDecimal
 import java.math.BigInteger
 import java.math.MathContext
 
-object BigIntegerRing : Ring<BigInteger> {
+object JBigIntegerField : Field<BigInteger> {
     override val zero: BigInteger = BigInteger.ZERO
     override val one: BigInteger = BigInteger.ONE
 
@@ -14,9 +14,11 @@ object BigIntegerRing : Ring<BigInteger> {
     override fun multiply(a: BigInteger, k: Number): BigInteger = a.multiply(k.toInt().toBigInteger())
 
     override fun multiply(a: BigInteger, b: BigInteger): BigInteger = a.multiply(b)
+
+    override fun divide(a: BigInteger, b: BigInteger): BigInteger = a.div(b)
 }
 
-class BigDecimalField(val mathContext: MathContext = MathContext.DECIMAL64) : Field<BigDecimal> {
+class JBigDecimalField(val mathContext: MathContext = MathContext.DECIMAL64) : Field<BigDecimal> {
     override val zero: BigDecimal = BigDecimal.ZERO
     override val one: BigDecimal = BigDecimal.ONE
 
@@ -28,18 +30,3 @@ class BigDecimalField(val mathContext: MathContext = MathContext.DECIMAL64) : Fi
     override fun multiply(a: BigDecimal, b: BigDecimal): BigDecimal = a.multiply(b, mathContext)
     override fun divide(a: BigDecimal, b: BigDecimal): BigDecimal = a.divide(b, mathContext)
 }
-
-inline fun Buffer.Companion.bigInt(size: Int, initializer: (Int) -> BigInteger): Buffer<BigInteger> =
-    boxing(size, initializer)
-
-inline fun MutableBuffer.Companion.bigInt(size: Int, initializer: (Int) -> BigInteger): MutableBuffer<BigInteger> =
-    boxing(size, initializer)
-
-fun NDAlgebra.Companion.bigInt(vararg shape: Int): BoxingNDRing<BigInteger, BigIntegerRing> =
-    BoxingNDRing(shape, BigIntegerRing, Buffer.Companion::bigInt)
-
-fun NDElement.Companion.bigInt(
-    vararg shape: Int,
-    initializer: BigIntegerRing.(IntArray) -> BigInteger
-): BufferedNDRingElement<BigInteger, BigIntegerRing> =
-    NDAlgebra.bigInt(*shape).produce(initializer)

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

@@ -0,0 +1,12 @@
+package scientifik.kmath.chains
+
+/**
+ * Performance optimized chain for integer values
+ */
+abstract class BlockingIntChain : Chain<Int> {
+    abstract fun nextInt(): Int
+
+    override suspend fun next(): Int = nextInt()
+
+    fun nextBlock(size: Int): IntArray = IntArray(size) { nextInt() }
+}

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

@@ -0,0 +1,12 @@
+package scientifik.kmath.chains
+
+/**
+ * Performance optimized chain for real values
+ */
+abstract class BlockingRealChain : Chain<Double> {
+    abstract fun nextDouble(): Double
+
+    override suspend fun next(): Double = nextDouble()
+
+    fun nextBlock(size: Int): DoubleArray = DoubleArray(size) { nextDouble() }
+}

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

@@ -38,13 +38,16 @@ interface Chain<out R>: Flow<R> {
      */
     fun fork(): Chain<R>
 
-    @InternalCoroutinesApi
+    @OptIn(InternalCoroutinesApi::class)
     override suspend fun collect(collector: FlowCollector<R>) {
-        kotlinx.coroutines.flow.flow { while (true) emit(next()) }.collect(collector)
+        kotlinx.coroutines.flow.flow {
+            while (true){
+                emit(next())
+            }
+        }.collect(collector)
     }
 
     companion object
-
 }
 
 
@@ -68,6 +71,8 @@ class MarkovChain<out R : Any>(private val seed: suspend () -> R, private val ge
 
     private var value: R? = null
 
+    fun value() = value
+
     override suspend fun next(): R {
         mutex.withLock {
             val newValue = gen(value ?: seed())
@@ -97,6 +102,8 @@ class StatefulChain<S, out R>(
 
     private var value: R? = null
 
+    fun value() = value
+
     override suspend fun next(): R {
         mutex.withLock {
             val newValue = state.gen(value ?: state.seed())

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

@@ -2,9 +2,11 @@ package scientifik.kmath.streaming
 
 import kotlinx.coroutines.FlowPreview
 import kotlinx.coroutines.flow.*
+import scientifik.kmath.chains.BlockingRealChain
 import scientifik.kmath.structures.Buffer
 import scientifik.kmath.structures.BufferFactory
 import scientifik.kmath.structures.DoubleBuffer
+import scientifik.kmath.structures.asBuffer
 
 /**
  * Create a [Flow] from buffer
@@ -45,6 +47,13 @@ fun <T> Flow<T>.chunked(bufferSize: Int, bufferFactory: BufferFactory<T>): Flow<
  */
 fun Flow<Double>.chunked(bufferSize: Int): Flow<DoubleBuffer> = flow {
     require(bufferSize > 0) { "Resulting chunk size must be more than zero" }
+
+    if (this@chunked is BlockingRealChain) {
+        //performance optimization for blocking primitive chain
+        while (true) {
+            emit(nextBlock(bufferSize).asBuffer())
+        }
+    } else {
         val array = DoubleArray(bufferSize)
         var counter = 0
 
@@ -60,6 +69,7 @@ fun Flow<Double>.chunked(bufferSize: Int): Flow<DoubleBuffer> = flow {
         if (counter > 0) {
             emit(DoubleBuffer(counter) { array[it] })
         }
+    }
 }
 
 /**

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

@@ -30,6 +30,20 @@ class LazyNDStructure<T>(
         }
         return res.asSequence()
     }
+
+    override fun equals(other: Any?): Boolean {
+        return NDStructure.equals(this, other as? NDStructure<*> ?: return false)
+    }
+
+    override fun hashCode(): Int {
+        var result = scope.hashCode()
+        result = 31 * result + shape.contentHashCode()
+        result = 31 * result + function.hashCode()
+        result = 31 * result + cache.hashCode()
+        return result
+    }
+
+
 }
 
 fun <T> NDStructure<T>.deferred(index: IntArray) =

kmath-for-real/src/commonMain/kotlin/scientifik/kmath/linear/RealVector.kt

@@ -1,48 +0,0 @@
-package scientifik.kmath.linear
-
-import scientifik.kmath.operations.Field
-import scientifik.kmath.operations.Norm
-import scientifik.kmath.operations.RealField
-import scientifik.kmath.operations.SpaceElement
-import scientifik.kmath.structures.Buffer
-import scientifik.kmath.structures.DoubleBuffer
-import scientifik.kmath.structures.asBuffer
-import scientifik.kmath.structures.asSequence
-
-fun DoubleArray.asVector() = RealVector(this.asBuffer())
-fun List<Double>.asVector() = RealVector(this.asBuffer())
-
-
-object VectorL2Norm : Norm<Point<out Number>, Double> {
-    override fun norm(arg: Point<out Number>): Double =
-        kotlin.math.sqrt(arg.asSequence().sumByDouble { it.toDouble() })
-}
-
-inline class RealVector(val point: Point<Double>) :
-    SpaceElement<Point<Double>, RealVector, VectorSpace<Double, RealField>>, Point<Double> {
-    override val context: VectorSpace<Double, RealField> get() = space(point.size)
-
-    override fun unwrap(): Point<Double> = point
-
-    override fun Point<Double>.wrap(): RealVector = RealVector(this)
-
-    override val size: Int get() = point.size
-
-    override fun get(index: Int): Double = point[index]
-
-    override fun iterator(): Iterator<Double> = point.iterator()
-
-    companion object {
-
-        private val spaceCache = HashMap<Int, BufferVectorSpace<Double, RealField>>()
-
-        inline operator fun invoke(dim:Int, initalizer: (Int)-> Double) = RealVector(DoubleBuffer(dim, initalizer))
-
-        operator fun invoke(vararg values: Double) = values.asVector()
-
-        fun space(dim: Int) =
-            spaceCache.getOrPut(dim) {
-                BufferVectorSpace(dim, RealField) { size, init -> Buffer.real(size, init) }
-            }
-    }
-}

kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/DoubleMatrixOperations.kt

@@ -1,146 +0,0 @@
-package scientifik.kmath.real
-
-import scientifik.kmath.linear.MatrixContext
-import scientifik.kmath.linear.RealMatrixContext.elementContext
-import scientifik.kmath.linear.VirtualMatrix
-import scientifik.kmath.operations.sum
-import scientifik.kmath.structures.Buffer
-import scientifik.kmath.structures.Matrix
-import scientifik.kmath.structures.asSequence
-import kotlin.math.pow
-
-/*
- *  Functions for convenient "numpy-like" operations with Double matrices.
- *
- *  Initial implementation of these functions is taken from:
- *    https://github.com/thomasnield/numky/blob/master/src/main/kotlin/org/nield/numky/linear/DoubleOperators.kt
- *
- */
-
-/*
- *  Functions that help create a real (Double) matrix
- */
-
-fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double) =
-        MatrixContext.real.produce(rowNum, colNum, initializer)
-
-fun Sequence<DoubleArray>.toMatrix() = toList().let {
-    MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] }
-}
-
-fun Matrix<Double>.repeatStackVertical(n: Int) = VirtualMatrix(rowNum*n, colNum) {
-    row, col -> get(if (row == 0) 0 else row % rowNum, col)
-}
-
-/*
- *  Operations for matrix and real number
- */
-
-operator fun Matrix<Double>.times(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
-    row, col -> this[row, col] * double
-}
-
-operator fun Matrix<Double>.plus(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
-    row, col -> this[row, col] + double
-}
-
-operator fun Matrix<Double>.minus(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
-    row, col -> this[row, col] - double
-}
-
-operator fun Matrix<Double>.div(double: Double) = MatrixContext.real.produce(rowNum, colNum) {
-    row, col -> this[row, col] / double
-}
-
-operator fun Double.times(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
-    row, col -> this * matrix[row, col]
-}
-
-operator fun Double.plus(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
-    row, col -> this + matrix[row, col]
-}
-
-operator fun Double.minus(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
-    row, col -> this - matrix[row, col]
-}
-
-// TODO: does this operation make sense? Should it be 'this/matrix[row, col]'?
-//operator fun Double.div(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
-//    row, col -> matrix[row, col] / this
-//}
-
-/*
- *  Per-element (!) square and power operations
- */
-
-fun Matrix<Double>.square() = MatrixContext.real.produce(rowNum, colNum) {
-    row, col -> this[row, col].pow(2)
-}
-
-fun Matrix<Double>.pow(n: Int) = MatrixContext.real.produce(rowNum, colNum) {
-    i, j -> this[i, j].pow(n)
-}
-
-/*
- * Operations on two matrices (per-element!)
- */
-
-operator fun Matrix<Double>.times(other: Matrix<Double>) = MatrixContext.real.produce(rowNum, colNum) {
-    row, col -> this[row, col] * other[row, col]
-}
-
-operator fun Matrix<Double>.plus(other: Matrix<Double>) = MatrixContext.real.add(this, other)
-
-operator fun Matrix<Double>.minus(other: Matrix<Double>) = MatrixContext.real.produce(rowNum, colNum) {
-    row, col -> this[row,col] - other[row,col]
-}
-
-/*
- *  Operations on columns
- */
-
-inline fun Matrix<Double>.appendColumn(crossinline mapper: (Buffer<Double>) -> Double) =
-    MatrixContext.real.produce(rowNum,colNum+1) {
-        row, col ->
-            if (col < colNum)
-                this[row, col]
-            else
-                mapper(rows[row])
-    }
-
-fun Matrix<Double>.extractColumns(columnRange: IntRange) = MatrixContext.real.produce(rowNum, columnRange.count()) {
-    row, col -> this[row, columnRange.first + col]
-}
-
-fun Matrix<Double>.extractColumn(columnIndex: Int) = extractColumns(columnIndex..columnIndex)
-
-fun Matrix<Double>.sumByColumn() = MatrixContext.real.produce(1, colNum) { _, j ->
-    val column = columns[j]
-    with(elementContext) {
-        sum(column.asSequence())
-    }
-}
-
-fun Matrix<Double>.minByColumn() = MatrixContext.real.produce(1, colNum) {
-    _, j -> columns[j].asSequence().min() ?: throw Exception("Cannot produce min on empty column")
-}
-
-fun Matrix<Double>.maxByColumn() = MatrixContext.real.produce(1, colNum) {
-    _, j -> columns[j].asSequence().max() ?: throw Exception("Cannot produce min on empty column")
-}
-
-fun Matrix<Double>.averageByColumn() = MatrixContext.real.produce(1, colNum) {
-    _, j -> columns[j].asSequence().average()
-}
-
-/*
- * Operations processing all elements
- */
-
-fun Matrix<Double>.sum() = elements().map { (_, value) -> value }.sum()
-
-fun Matrix<Double>.min() = elements().map { (_, value) -> value }.min()
-
-fun Matrix<Double>.max() = elements().map { (_, value) -> value }.max()
-
-fun Matrix<Double>.average() = elements().map { (_, value) -> value }.average()

kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/RealVector.kt

@@ -0,0 +1,59 @@
+package scientifik.kmath.real
+
+import scientifik.kmath.linear.BufferVectorSpace
+import scientifik.kmath.linear.Point
+import scientifik.kmath.linear.VectorSpace
+import scientifik.kmath.operations.Norm
+import scientifik.kmath.operations.RealField
+import scientifik.kmath.operations.SpaceElement
+import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.DoubleBuffer
+import scientifik.kmath.structures.asBuffer
+import scientifik.kmath.structures.asIterable
+import kotlin.math.sqrt
+
+fun DoubleArray.asVector() = RealVector(this.asBuffer())
+fun List<Double>.asVector() = RealVector(this.asBuffer())
+
+
+object VectorL2Norm : Norm<Point<out Number>, Double> {
+    override fun norm(arg: Point<out Number>): Double = sqrt(arg.asIterable().sumByDouble { it.toDouble() })
+}
+
+inline class RealVector(private val point: Point<Double>) :
+    SpaceElement<Point<Double>, RealVector, VectorSpace<Double, RealField>>, Point<Double> {
+
+    override val context: VectorSpace<Double, RealField>
+        get() = space(
+            point.size
+        )
+
+    override fun unwrap(): Point<Double> = point
+
+    override fun Point<Double>.wrap(): RealVector =
+        RealVector(this)
+
+    override val size: Int get() = point.size
+
+    override fun get(index: Int): Double = point[index]
+
+    override fun iterator(): Iterator<Double> = point.iterator()
+
+    companion object {
+
+        private val spaceCache = HashMap<Int, BufferVectorSpace<Double, RealField>>()
+
+        inline operator fun invoke(dim: Int, initializer: (Int) -> Double) =
+            RealVector(DoubleBuffer(dim, initializer))
+
+        operator fun invoke(vararg values: Double): RealVector = values.asVector()
+
+        fun space(dim: Int): BufferVectorSpace<Double, RealField> =
+            spaceCache.getOrPut(dim) {
+                BufferVectorSpace(
+                    dim,
+                    RealField
+                ) { size, init -> Buffer.real(size, init) }
+            }
+    }
+}

kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/realBuffer.kt

@@ -0,0 +1,8 @@
+package scientifik.kmath.real
+
+import scientifik.kmath.structures.DoubleBuffer
+
+/**
+ * Simplified [DoubleBuffer] to array comparison
+ */
+fun DoubleBuffer.contentEquals(vararg doubles: Double) = array.contentEquals(doubles)

kmath-for-real/src/commonMain/kotlin/scientifik/kmath/real/realMatrix.kt

@@ -0,0 +1,161 @@
+package scientifik.kmath.real
+
+import scientifik.kmath.linear.MatrixContext
+import scientifik.kmath.linear.RealMatrixContext.elementContext
+import scientifik.kmath.linear.VirtualMatrix
+import scientifik.kmath.operations.sum
+import scientifik.kmath.structures.Buffer
+import scientifik.kmath.structures.DoubleBuffer
+import scientifik.kmath.structures.Matrix
+import scientifik.kmath.structures.asIterable
+import kotlin.math.pow
+
+/*
+ *  Functions for convenient "numpy-like" operations with Double matrices.
+ *
+ *  Initial implementation of these functions is taken from:
+ *    https://github.com/thomasnield/numky/blob/master/src/main/kotlin/org/nield/numky/linear/DoubleOperators.kt
+ *
+ */
+
+/*
+ *  Functions that help create a real (Double) matrix
+ */
+
+typealias RealMatrix = Matrix<Double>
+
+fun realMatrix(rowNum: Int, colNum: Int, initializer: (i: Int, j: Int) -> Double): RealMatrix =
+    MatrixContext.real.produce(rowNum, colNum, initializer)
+
+fun Sequence<DoubleArray>.toMatrix(): RealMatrix = toList().let {
+    MatrixContext.real.produce(it.size, it[0].size) { row, col -> it[row][col] }
+}
+
+fun Matrix<Double>.repeatStackVertical(n: Int): RealMatrix =
+    VirtualMatrix(rowNum * n, colNum) { row, col ->
+        get(if (row == 0) 0 else row % rowNum, col)
+    }
+
+/*
+ *  Operations for matrix and real number
+ */
+
+operator fun Matrix<Double>.times(double: Double): RealMatrix =
+    MatrixContext.real.produce(rowNum, colNum) { row, col ->
+        this[row, col] * double
+    }
+
+operator fun Matrix<Double>.plus(double: Double): RealMatrix =
+    MatrixContext.real.produce(rowNum, colNum) { row, col ->
+        this[row, col] + double
+    }
+
+operator fun Matrix<Double>.minus(double: Double): RealMatrix =
+    MatrixContext.real.produce(rowNum, colNum) { row, col ->
+        this[row, col] - double
+    }
+
+operator fun Matrix<Double>.div(double: Double): RealMatrix =
+    MatrixContext.real.produce(rowNum, colNum) { row, col ->
+        this[row, col] / double
+    }
+
+operator fun Double.times(matrix: Matrix<Double>): RealMatrix =
+    MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
+        this * matrix[row, col]
+    }
+
+operator fun Double.plus(matrix: Matrix<Double>): RealMatrix =
+    MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
+        this + matrix[row, col]
+    }
+
+operator fun Double.minus(matrix: Matrix<Double>): RealMatrix =
+    MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
+        this - matrix[row, col]
+    }
+
+// TODO: does this operation make sense? Should it be 'this/matrix[row, col]'?
+//operator fun Double.div(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) {
+//    row, col -> matrix[row, col] / this
+//}
+
+/*
+ *  Per-element (!) square and power operations
+ */
+
+fun Matrix<Double>.square(): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { row, col ->
+    this[row, col].pow(2)
+}
+
+fun Matrix<Double>.pow(n: Int): RealMatrix = MatrixContext.real.produce(rowNum, colNum) { i, j ->
+    this[i, j].pow(n)
+}
+
+/*
+ * Operations on two matrices (per-element!)
+ */
+
+operator fun Matrix<Double>.times(other: Matrix<Double>): RealMatrix =
+    MatrixContext.real.produce(rowNum, colNum) { row, col ->
+        this[row, col] * other[row, col]
+    }
+
+operator fun Matrix<Double>.plus(other: Matrix<Double>): RealMatrix =
+    MatrixContext.real.add(this, other)
+
+operator fun Matrix<Double>.minus(other: Matrix<Double>): RealMatrix =
+    MatrixContext.real.produce(rowNum, colNum) { row, col ->
+        this[row, col] - other[row, col]
+    }
+
+/*
+ *  Operations on columns
+ */
+
+inline fun Matrix<Double>.appendColumn(crossinline mapper: (Buffer<Double>) -> Double) =
+    MatrixContext.real.produce(rowNum, colNum + 1) { row, col ->
+        if (col < colNum)
+            this[row, col]
+        else
+            mapper(rows[row])
+    }
+
+fun Matrix<Double>.extractColumns(columnRange: IntRange): RealMatrix =
+    MatrixContext.real.produce(rowNum, columnRange.count()) { row, col ->
+        this[row, columnRange.first + col]
+    }
+
+fun Matrix<Double>.extractColumn(columnIndex: Int): RealMatrix =
+    extractColumns(columnIndex..columnIndex)
+
+fun Matrix<Double>.sumByColumn(): DoubleBuffer = DoubleBuffer(colNum) { j ->
+    val column = columns[j]
+    with(elementContext) {
+        sum(column.asIterable())
+    }
+}
+
+fun Matrix<Double>.minByColumn(): DoubleBuffer = DoubleBuffer(colNum) { j ->
+    columns[j].asIterable().min() ?: throw Exception("Cannot produce min on empty column")
+}
+
+fun Matrix<Double>.maxByColumn(): DoubleBuffer = DoubleBuffer(colNum) { j ->
+    columns[j].asIterable().max() ?: throw Exception("Cannot produce min on empty column")
+}
+
+fun Matrix<Double>.averageByColumn(): DoubleBuffer = DoubleBuffer(colNum) { j ->
+    columns[j].asIterable().average()
+}
+
+/*
+ * Operations processing all elements
+ */
+
+fun Matrix<Double>.sum() = elements().map { (_, value) -> value }.sum()
+
+fun Matrix<Double>.min() = elements().map { (_, value) -> value }.min()
+
+fun Matrix<Double>.max() = elements().map { (_, value) -> value }.max()
+
+fun Matrix<Double>.average() = elements().map { (_, value) -> value }.average()

kmath-for-real/src/commonTest/kotlin/scientific.kmath.real/RealMatrixTest.kt

@@ -1,11 +1,12 @@
 package scientific.kmath.real
 
-import scientifik.kmath.real.*
 import scientifik.kmath.linear.VirtualMatrix
 import scientifik.kmath.linear.build
+import scientifik.kmath.real.*
 import scientifik.kmath.structures.Matrix
 import kotlin.test.Test
 import kotlin.test.assertEquals
+import kotlin.test.assertTrue
 
 class RealMatrixTest {
     @Test
@@ -28,11 +29,11 @@ class RealMatrixTest {
 
     @Test
     fun testRepeatStackVertical() {
-        val matrix1 = Matrix.build<Double>(2, 3)(
+        val matrix1 = Matrix.build(2, 3)(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0
         )
-        val matrix2 = Matrix.build<Double>(6, 3)(
+        val matrix2 = Matrix.build(6, 3)(
             1.0, 0.0, 0.0,
             0.0, 1.0, 2.0,
             1.0, 0.0, 0.0,
@@ -45,12 +46,12 @@ class RealMatrixTest {
 
     @Test
     fun testMatrixAndDouble() {
-        val matrix1 = Matrix.build<Double>(2, 3)(
+        val matrix1 = Matrix.build(2, 3)(
             1.0, 0.0, 3.0,
             4.0, 6.0, 2.0
         )
         val matrix2 = (matrix1 * 2.5 + 1.0 - 2.0) / 2.0
-        val expectedResult = Matrix.build<Double>(2, 3)(
+        val expectedResult = Matrix.build(2, 3)(
             0.75, -0.5, 3.25,
             4.5, 7.0, 2.0
         )
@@ -59,13 +60,13 @@ class RealMatrixTest {
 
     @Test
     fun testDoubleAndMatrix() {
-        val matrix1 = Matrix.build<Double>(2, 3)(
+        val matrix1 = Matrix.build(2, 3)(
             1.0, 0.0, 3.0,
             4.0, 6.0, 2.0
         )
         val matrix2 = 20.0 - (10.0 + (5.0 * matrix1))
         //val matrix2 = 10.0 + (5.0 * matrix1)
-        val expectedResult = Matrix.build<Double>(2, 3)(
+        val expectedResult = Matrix.build(2, 3)(
             5.0, 10.0, -5.0,
             -10.0, -20.0, 0.0
         )
@@ -74,15 +75,15 @@ class RealMatrixTest {
 
     @Test
     fun testSquareAndPower() {
-        val matrix1 = Matrix.build<Double>(2, 3)(
+        val matrix1 = Matrix.build(2, 3)(
             -1.0, 0.0, 3.0,
             4.0, -6.0, -2.0
         )
-        val matrix2 = Matrix.build<Double>(2, 3)(
+        val matrix2 = Matrix.build(2, 3)(
             1.0, 0.0, 9.0,
             16.0, 36.0, 4.0
         )
-        val matrix3 = Matrix.build<Double>(2, 3)(
+        val matrix3 = Matrix.build(2, 3)(
             -1.0, 0.0, 27.0,
             64.0, -216.0, -8.0
         )
@@ -92,16 +93,16 @@ class RealMatrixTest {
 
     @Test
     fun testTwoMatrixOperations() {
-        val matrix1 = Matrix.build<Double>(2, 3)(
+        val matrix1 = Matrix.build(2, 3)(
             -1.0, 0.0, 3.0,
             4.0, -6.0, 7.0
         )
-        val matrix2 = Matrix.build<Double>(2, 3)(
+        val matrix2 = Matrix.build(2, 3)(
             1.0, 0.0, 3.0,
             4.0, 6.0, -2.0
         )
         val result = matrix1 * matrix2 + matrix1 - matrix2
-        val expectedResult = Matrix.build<Double>(2, 3)(
+        val expectedResult = Matrix.build(2, 3)(
             -3.0, 0.0, 9.0,
             16.0, -48.0, -5.0
         )
@@ -110,31 +111,41 @@ class RealMatrixTest {
 
     @Test
     fun testColumnOperations() {
-        val matrix1 = Matrix.build<Double>(2, 4)(
+        val matrix1 = Matrix.build(2, 4)(
             -1.0, 0.0, 3.0, 15.0,
             4.0, -6.0, 7.0, -11.0
         )
-        val matrix2 = Matrix.build<Double>(2, 5)(
+        val matrix2 = Matrix.build(2, 5)(
             -1.0, 0.0, 3.0, 15.0, -1.0,
             4.0, -6.0, 7.0, -11.0, 4.0
         )
-        val col1 = Matrix.build<Double>(2, 1)(0.0, -6.0)
+        val col1 = Matrix.build(2, 1)(0.0, -6.0)
-        val cols1to2 = Matrix.build<Double>(2, 2)(
+        val cols1to2 = Matrix.build(2, 2)(
             0.0, 3.0,
             -6.0, 7.0
         )
+
         assertEquals(matrix1.appendColumn { it[0] }, matrix2)
         assertEquals(matrix1.extractColumn(1), col1)
         assertEquals(matrix1.extractColumns(1..2), cols1to2)
-        assertEquals(matrix1.sumByColumn(), Matrix.build<Double>(4, 1)(3.0, -6.0, 10.0, 4.0))
+        //equals should never be called on buffers
-        assertEquals(matrix1.minByColumn(), Matrix.build<Double>(4, 1)(-1.0, -6.0, 3.0, -11.0))
+        assertTrue {
-        assertEquals(matrix1.maxByColumn(), Matrix.build<Double>(4, 1)(4.0, 0.0, 7.0, 15.0))
+            matrix1.sumByColumn().contentEquals(3.0, -6.0, 10.0, 4.0)
-        assertEquals(matrix1.averageByColumn(), Matrix.build<Double>(4, 1)(1.5, -3.0, 5.0, 2.0))
+        } //assertEquals(matrix1.sumByColumn(), DoubleBuffer(3.0, -6.0, 10.0, 4.0))
+        assertTrue {
+            matrix1.minByColumn().contentEquals(-1.0, -6.0, 3.0, -11.0)
+        } //assertEquals(matrix1.minByColumn(), DoubleBuffer(-1.0, -6.0, 3.0, -11.0))
+        assertTrue {
+            matrix1.maxByColumn().contentEquals(4.0, 0.0, 7.0, 15.0)
+        } //assertEquals(matrix1.maxByColumn(), DoubleBuffer(4.0, 0.0, 7.0, 15.0))
+        assertTrue {
+            matrix1.averageByColumn().contentEquals(1.5, -3.0, 5.0, 2.0)
+        } //assertEquals(matrix1.averageByColumn(), DoubleBuffer(1.5, -3.0, 5.0, 2.0))
     }
 
     @Test
     fun testAllElementOperations() {
-        val matrix1 = Matrix.build<Double>(2, 4)(
+        val matrix1 = Matrix.build(2, 4)(
             -1.0, 0.0, 3.0, 15.0,
             4.0, -6.0, 7.0, -11.0
         )

kmath-for-real/src/commonTest/kotlin/scientifik/kmath/linear/VectorTest.kt

@@ -1,5 +1,6 @@
 package scientifik.kmath.linear
 
+import scientifik.kmath.real.RealVector
 import kotlin.test.Test
 import kotlin.test.assertEquals
 

kmath-histograms/src/commonMain/kotlin/scientifik/kmath/histogram/RealHistogram.kt

@@ -1,7 +1,7 @@
 package scientifik.kmath.histogram
 
 import scientifik.kmath.linear.Point
-import scientifik.kmath.linear.asVector
+import scientifik.kmath.real.asVector
 import scientifik.kmath.operations.SpaceOperations
 import scientifik.kmath.structures.*
 import kotlin.math.floor

kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt

@@ -3,7 +3,7 @@ package scietifik.kmath.histogram
 import scientifik.kmath.histogram.RealHistogram
 import scientifik.kmath.histogram.fill
 import scientifik.kmath.histogram.put
-import scientifik.kmath.linear.RealVector
+import scientifik.kmath.real.RealVector
 import kotlin.random.Random
 import kotlin.test.Test
 import kotlin.test.assertEquals

kmath-histograms/src/jvmMain/kotlin/scientifik/kmath/histogram/UnivariateHistogram.kt

@@ -1,7 +1,7 @@
 package scientifik.kmath.histogram
 
-import scientifik.kmath.linear.RealVector
+import scientifik.kmath.real.RealVector
-import scientifik.kmath.linear.asVector
+import scientifik.kmath.real.asVector
 import scientifik.kmath.structures.Buffer
 import java.util.*
 import kotlin.math.floor

kmath-koma/src/commonMain/kotlin/scientifik.kmath.linear/KomaMatrix.kt

@@ -4,6 +4,7 @@ import koma.extensions.fill
 import koma.matrix.MatrixFactory
 import scientifik.kmath.operations.Space
 import scientifik.kmath.structures.Matrix
+import scientifik.kmath.structures.NDStructure
 
 class KomaMatrixContext<T : Any>(
     private val factory: MatrixFactory<koma.matrix.Matrix<T>>,
@@ -85,6 +86,18 @@ class KomaMatrix<T : Any>(val origin: koma.matrix.Matrix<T>, features: Set<Matri
         KomaMatrix(this.origin, this.features + features)
 
     override fun get(i: Int, j: Int): T = origin.getGeneric(i, j)
+
+    override fun equals(other: Any?): Boolean {
+        return NDStructure.equals(this, other as? NDStructure<*> ?: return false)
+    }
+
+    override fun hashCode(): Int {
+        var result = origin.hashCode()
+        result = 31 * result + features.hashCode()
+        return result
+    }
+
+
 }
 
 class KomaVector<T : Any> internal constructor(val origin: koma.matrix.Matrix<T>) : Point<T> {

kmath-memory/src/commonMain/kotlin/scientifik/memory/Memory.kt

@@ -69,3 +69,9 @@ inline fun Memory.write(block: MemoryWriter.() -> Unit) {
  * Allocate the most effective platform-specific memory
  */
 expect fun Memory.Companion.allocate(length: Int): Memory
+
+/**
+ * Wrap a [Memory] around existing [ByteArray]. This operation is unsafe since the array is not copied
+ * and could be mutated independently from the resulting [Memory]
+ */
+expect fun Memory.Companion.wrap(array: ByteArray): Memory

kmath-memory/src/jsMain/kotlin/scientifik/memory/DataViewMemory.kt

@@ -2,14 +2,7 @@ package scientifik.memory
 
 import org.khronos.webgl.ArrayBuffer
 import org.khronos.webgl.DataView
-
+import org.khronos.webgl.Int8Array
-/**
- * Allocate the most effective platform-specific memory
- */
-actual fun Memory.Companion.allocate(length: Int): Memory {
-    val buffer = ArrayBuffer(length)
-    return DataViewMemory(DataView(buffer, 0, length))
-}
 
 class DataViewMemory(val view: DataView) : Memory {
 
@@ -89,3 +82,16 @@ class DataViewMemory(val view: DataView) : Memory {
     override fun writer(): MemoryWriter = writer
 
 }
+
+/**
+ * Allocate the most effective platform-specific memory
+ */
+actual fun Memory.Companion.allocate(length: Int): Memory {
+    val buffer = ArrayBuffer(length)
+    return DataViewMemory(DataView(buffer, 0, length))
+}
+
+actual fun Memory.Companion.wrap(array: ByteArray): Memory {
+    @Suppress("CAST_NEVER_SUCCEEDS") val int8Array = array as Int8Array
+    return DataViewMemory(DataView(int8Array.buffer, int8Array.byteOffset, int8Array.length))
+}

kmath-memory/src/jvmMain/kotlin/scientifik/memory/ByteBufferMemory.kt

@@ -7,14 +7,6 @@ import java.nio.file.Path
 import java.nio.file.StandardOpenOption
 
 
-/**
- * Allocate the most effective platform-specific memory
- */
-actual fun Memory.Companion.allocate(length: Int): Memory {
-    val buffer = ByteBuffer.allocate(length)
-    return ByteBufferMemory(buffer)
-}
-
 private class ByteBufferMemory(
     val buffer: ByteBuffer,
     val startOffset: Int = 0,
@@ -96,6 +88,22 @@ private class ByteBufferMemory(
     override fun writer(): MemoryWriter = writer
 }
 
+/**
+ * Allocate the most effective platform-specific memory
+ */
+actual fun Memory.Companion.allocate(length: Int): Memory {
+    val buffer = ByteBuffer.allocate(length)
+    return ByteBufferMemory(buffer)
+}
+
+actual fun Memory.Companion.wrap(array: ByteArray): Memory {
+    val buffer = ByteBuffer.wrap(array)
+    return ByteBufferMemory(buffer)
+}
+
+fun ByteBuffer.asMemory(startOffset: Int = 0, size: Int = limit()): Memory =
+    ByteBufferMemory(this, startOffset, size)
+
 /**
  * Use direct memory-mapped buffer from file to read something and close it afterwards.
  */

kmath-prob/build.gradle.kts

@@ -8,4 +8,10 @@ kotlin.sourceSets {
             api(project(":kmath-coroutines"))
         }
     }
+    jvmMain{
+        dependencies{
+            api("org.apache.commons:commons-rng-sampling:1.3")
+            api("org.apache.commons:commons-rng-simple:1.3")
+        }
+    }
 }

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

@@ -6,10 +6,16 @@ import kotlin.random.Random
  * A basic generator
  */
 interface RandomGenerator {
+    fun nextBoolean(): Boolean
+
     fun nextDouble(): Double
     fun nextInt(): Int
+    fun nextInt(until: Int): Int
     fun nextLong(): Long
-    fun nextBlock(size: Int): ByteArray
+    fun nextLong(until: Long): Long
+
+    fun fillBytes(array: ByteArray, fromIndex: Int = 0, toIndex: Int = array.size)
+    fun nextBytes(size: Int): ByteArray = ByteArray(size).also { fillBytes(it) }
 
     /**
      * Create a new generator which is independent from current generator (operations on new generator do not affect this one
@@ -21,21 +27,29 @@ interface RandomGenerator {
     fun fork(): RandomGenerator
 
     companion object {
-        val default by lazy { DefaultGenerator(Random.nextLong()) }
+        val default by lazy { DefaultGenerator() }
+
+        fun default(seed: Long) = DefaultGenerator(Random(seed))
     }
 }
 
-class DefaultGenerator(seed: Long?) : RandomGenerator {
+inline class DefaultGenerator(val random: Random = Random) : RandomGenerator {
-    private val random = seed?.let { Random(it) } ?: Random
+    override fun nextBoolean(): Boolean = random.nextBoolean()
 
     override fun nextDouble(): Double = random.nextDouble()
 
     override fun nextInt(): Int = random.nextInt()
+    override fun nextInt(until: Int): Int = random.nextInt(until)
 
     override fun nextLong(): Long = random.nextLong()
 
-    override fun nextBlock(size: Int): ByteArray = random.nextBytes(size)
+    override fun nextLong(until: Long): Long = random.nextLong(until)
 
-    override fun fork(): RandomGenerator = DefaultGenerator(nextLong())
+    override fun fillBytes(array: ByteArray, fromIndex: Int, toIndex: Int) {
+        random.nextBytes(array, fromIndex, toIndex)
+    }
 
+    override fun nextBytes(size: Int): ByteArray = random.nextBytes(size)
+
+    override fun fork(): RandomGenerator = RandomGenerator.default(random.nextLong())
 }

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

@@ -29,3 +29,6 @@ class UniformDistribution(val range: ClosedFloatingPointRange<Double>) : Univari
         }
     }
 }
+
+fun Distribution.Companion.uniform(range: ClosedFloatingPointRange<Double>): UniformDistribution =
+    UniformDistribution(range)

kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/RandomSourceGenerator.kt

@@ -0,0 +1,67 @@
+package scientifik.kmath.prob
+
+import org.apache.commons.rng.UniformRandomProvider
+import org.apache.commons.rng.simple.RandomSource
+
+class RandomSourceGenerator(val source: RandomSource, seed: Long?) : RandomGenerator {
+    internal val random: UniformRandomProvider = seed?.let {
+        RandomSource.create(source, seed)
+    } ?: RandomSource.create(source)
+
+    override fun nextBoolean(): Boolean = random.nextBoolean()
+
+    override fun nextDouble(): Double = random.nextDouble()
+
+    override fun nextInt(): Int = random.nextInt()
+    override fun nextInt(until: Int): Int = random.nextInt(until)
+
+    override fun nextLong(): Long = random.nextLong()
+    override fun nextLong(until: Long): Long = random.nextLong(until)
+
+    override fun fillBytes(array: ByteArray, fromIndex: Int, toIndex: Int) {
+        require(toIndex > fromIndex)
+        random.nextBytes(array, fromIndex, toIndex - fromIndex)
+    }
+
+    override fun fork(): RandomGenerator = RandomSourceGenerator(source, nextLong())
+}
+
+inline class RandomGeneratorProvider(val generator: RandomGenerator) : UniformRandomProvider {
+    override fun nextBoolean(): Boolean = generator.nextBoolean()
+
+    override fun nextFloat(): Float = generator.nextDouble().toFloat()
+
+    override fun nextBytes(bytes: ByteArray) {
+        generator.fillBytes(bytes)
+    }
+
+    override fun nextBytes(bytes: ByteArray, start: Int, len: Int) {
+        generator.fillBytes(bytes, start, start + len)
+    }
+
+    override fun nextInt(): Int = generator.nextInt()
+
+    override fun nextInt(n: Int): Int = generator.nextInt(n)
+
+    override fun nextDouble(): Double = generator.nextDouble()
+
+    override fun nextLong(): Long = generator.nextLong()
+
+    override fun nextLong(n: Long): Long = generator.nextLong(n)
+}
+
+/**
+ * Represent this [RandomGenerator] as commons-rng [UniformRandomProvider] preserving and mirroring its current state.
+ * Getting new value from one of those changes the state of another.
+ */
+fun RandomGenerator.asUniformRandomProvider(): UniformRandomProvider = if (this is RandomSourceGenerator) {
+    random
+} else {
+    RandomGeneratorProvider(this)
+}
+
+fun RandomGenerator.Companion.fromSource(source: RandomSource, seed: Long? = null): RandomSourceGenerator =
+    RandomSourceGenerator(source, seed)
+
+fun RandomGenerator.Companion.mersenneTwister(seed: Long? = null): RandomSourceGenerator =
+    fromSource(RandomSource.MT, seed)

kmath-prob/src/jvmMain/kotlin/scientifik/kmath/prob/distributions.kt

@@ -0,0 +1,109 @@
+package scientifik.kmath.prob
+
+import org.apache.commons.rng.UniformRandomProvider
+import org.apache.commons.rng.sampling.distribution.*
+import scientifik.kmath.chains.BlockingIntChain
+import scientifik.kmath.chains.BlockingRealChain
+import scientifik.kmath.chains.Chain
+import java.util.*
+import kotlin.math.PI
+import kotlin.math.exp
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+abstract class ContinuousSamplerDistribution : Distribution<Double> {
+
+    private inner class ContinuousSamplerChain(val generator: RandomGenerator) : BlockingRealChain() {
+        private val sampler = buildCMSampler(generator)
+
+        override fun nextDouble(): Double = sampler.sample()
+
+        override fun fork(): Chain<Double> = ContinuousSamplerChain(generator.fork())
+    }
+
+    protected abstract fun buildCMSampler(generator: RandomGenerator): ContinuousSampler
+
+    override fun sample(generator: RandomGenerator): BlockingRealChain = ContinuousSamplerChain(generator)
+}
+
+abstract class DiscreteSamplerDistribution : Distribution<Int> {
+
+    private inner class ContinuousSamplerChain(val generator: RandomGenerator) : BlockingIntChain() {
+        private val sampler = buildSampler(generator)
+
+        override fun nextInt(): Int = sampler.sample()
+
+        override fun fork(): Chain<Int> = ContinuousSamplerChain(generator.fork())
+    }
+
+    protected abstract fun buildSampler(generator: RandomGenerator): DiscreteSampler
+
+    override fun sample(generator: RandomGenerator): BlockingIntChain = ContinuousSamplerChain(generator)
+}
+
+enum class NormalSamplerMethod {
+    BoxMuller,
+    Marsaglia,
+    Ziggurat
+}
+
+private fun normalSampler(method: NormalSamplerMethod, provider: UniformRandomProvider): NormalizedGaussianSampler =
+    when (method) {
+        NormalSamplerMethod.BoxMuller -> BoxMullerNormalizedGaussianSampler(provider)
+        NormalSamplerMethod.Marsaglia -> MarsagliaNormalizedGaussianSampler(provider)
+        NormalSamplerMethod.Ziggurat -> ZigguratNormalizedGaussianSampler(provider)
+    }
+
+fun Distribution.Companion.normal(
+    method: NormalSamplerMethod = NormalSamplerMethod.Ziggurat
+): Distribution<Double> = object : ContinuousSamplerDistribution() {
+    override fun buildCMSampler(generator: RandomGenerator): ContinuousSampler {
+        val provider: UniformRandomProvider = generator.asUniformRandomProvider()
+        return normalSampler(method, provider)
+    }
+
+    override fun probability(arg: Double): Double {
+        return exp(-arg.pow(2) / 2) / sqrt(PI * 2)
+    }
+}
+
+fun Distribution.Companion.normal(
+    mean: Double,
+    sigma: Double,
+    method: NormalSamplerMethod = NormalSamplerMethod.Ziggurat
+): ContinuousSamplerDistribution = object : ContinuousSamplerDistribution() {
+    private val sigma2 = sigma.pow(2)
+    private val norm = sigma * sqrt(PI * 2)
+
+    override fun buildCMSampler(generator: RandomGenerator): ContinuousSampler {
+        val provider: UniformRandomProvider = generator.asUniformRandomProvider()
+        val normalizedSampler = normalSampler(method, provider)
+        return GaussianSampler(normalizedSampler, mean, sigma)
+    }
+
+    override fun probability(arg: Double): Double {
+        return exp(-(arg - mean).pow(2) / 2 / sigma2) / norm
+    }
+}
+
+fun Distribution.Companion.poisson(
+    lambda: Double
+): DiscreteSamplerDistribution = object : DiscreteSamplerDistribution() {
+
+    override fun buildSampler(generator: RandomGenerator): DiscreteSampler {
+        return PoissonSampler.of(generator.asUniformRandomProvider(), lambda)
+    }
+
+    private val computedProb: HashMap<Int, Double> = hashMapOf(0 to exp(-lambda))
+
+    override fun probability(arg: Int): Double {
+        require(arg >= 0) { "The argument must be >= 0" }
+        return if (arg > 40) {
+            exp(-(arg - lambda).pow(2) / 2 / lambda) / sqrt(2 * PI * lambda)
+        } else {
+            computedProb.getOrPut(arg) {
+                probability(arg - 1) * lambda / arg
+            }
+        }
+    }
+}

kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/CommonsDistributionsTest.kt

@@ -0,0 +1,28 @@
+package scientifik.kmath.prob
+
+import kotlinx.coroutines.flow.take
+import kotlinx.coroutines.flow.toList
+import kotlinx.coroutines.runBlocking
+import org.junit.jupiter.api.Assertions
+import org.junit.jupiter.api.Test
+
+class CommonsDistributionsTest {
+    @Test
+    fun testNormalDistributionSuspend() {
+        val distribution = Distribution.normal(7.0, 2.0)
+        val generator = RandomGenerator.default(1)
+        val sample = runBlocking {
+            distribution.sample(generator).take(1000).toList()
+        }
+        Assertions.assertEquals(7.0, sample.average(), 0.1)
+    }
+
+    @Test
+    fun testNormalDistributionBlocking() {
+        val distribution = Distribution.normal(7.0, 2.0)
+        val generator = RandomGenerator.default(1)
+        val sample = distribution.sample(generator).nextBlock(1000)
+        Assertions.assertEquals(7.0, sample.average(), 0.1)
+    }
+
+}

kmath-prob/src/jvmTest/kotlin/scientifik/kmath/prob/StatisticTest.kt

@@ -9,7 +9,7 @@ import kotlin.test.Test
 
 class StatisticTest {
     //create a random number generator.
-    val generator = DefaultGenerator(1)
+    val generator = RandomGenerator.default(1)
     //Create a stateless chain from generator.
     val data = generator.chain { nextDouble() }
     //Convert a chaint to Flow and break it into chunks.

settings.gradle.kts

@@ -1,10 +1,12 @@
 pluginManagement {
 
+    val toolsVersion = "0.5.0"
+
     plugins {
-        id("scientifik.mpp") version "0.4.1"
+        id("scientifik.mpp") version toolsVersion
-        id("scientifik.jvm") version "0.4.1"
+        id("scientifik.jvm") version toolsVersion
-        id("scientifik.atomic") version "0.4.1"
+        id("scientifik.atomic") version toolsVersion
-        id("scientifik.publish") version "0.4.1"
+        id("scientifik.publish") version toolsVersion
     }
 
     repositories {
@@ -20,7 +22,7 @@ pluginManagement {
     resolutionStrategy {
         eachPlugin {
             when (requested.id.id) {
-                "scientifik.mpp", "scientifik.jvm", "scientifik.publish" -> useModule("scientifik:gradle-tools:${requested.version}")
+                "scientifik.mpp", "scientifik.jvm", "scientifik.publish" -> useModule("scientifik:gradle-tools:$toolsVersion")
             }
         }
     }