Merge branch 'feature/series' into dev

# Conflicts:
#	build.gradle.kts
#	kmath-commons/src/main/kotlin/space/kscience/kmath/commons/transform/Transformations.kt
#	kmath-core/build.gradle.kts
#	kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/AlgebraND.kt
#	kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferAlgebraND.kt
#	kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/BufferND.kt
#	kmath-core/src/commonMain/kotlin/space/kscience/kmath/nd/ShapeIndices.kt
#	kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/BufferAlgebra.kt
#	kmath-core/src/commonMain/kotlin/space/kscience/kmath/operations/DoubleBufferOps.kt
#	kmath-core/src/commonMain/kotlin/space/kscience/kmath/structures/bufferExtensions.kt
#	kmath-histograms/src/jvmMain/kotlin/space/kscience/kmath/histogram/UnivariateHistogram.kt
#	kmath-stat/src/commonMain/kotlin/space/kscience/kmath/distributions/UniformDistribution.kt
#	kmath-stat/src/commonMain/kotlin/space/kscience/kmath/samplers/Sampler.kt
#	kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BroadcastDoubleTensorAlgebra.kt
#	kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/BufferedTensor.kt
#	kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/DoubleTensorAlgebra.kt
#	kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/broadcastUtils.kt
#	kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/linUtils.kt
#	kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/tensorCastsUtils.kt
#	kmath-tensors/src/commonMain/kotlin/space/kscience/kmath/tensors/core/internal/utils.kt
#	settings.gradle.kts

benchmarks/src/jvmMain/kotlin/space/kscience/kmath/benchmarks/BufferBenchmark.kt

@@ -6,34 +6,75 @@
 package space.kscience.kmath.benchmarks
 
 import kotlinx.benchmark.Benchmark
+import kotlinx.benchmark.Blackhole
 import kotlinx.benchmark.Scope
 import kotlinx.benchmark.State
 import space.kscience.kmath.complex.Complex
+import space.kscience.kmath.complex.ComplexField
 import space.kscience.kmath.complex.complex
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.MutableBuffer
+import space.kscience.kmath.structures.getDouble
+import space.kscience.kmath.structures.permute
 
 @State(Scope.Benchmark)
 internal class BufferBenchmark {
-    @Benchmark
-    fun genericDoubleBufferReadWrite() {
-        val buffer = DoubleBuffer(size) { it.toDouble() }
 
+    @Benchmark
+    fun doubleArrayReadWrite(blackhole: Blackhole) {
+        val buffer = DoubleArray(size) { it.toDouble() }
+        var res = 0.0
         (0 until size).forEach {
-            buffer[it]
+            res += buffer[it]
         }
+        blackhole.consume(res)
     }
 
     @Benchmark
-    fun complexBufferReadWrite() {
+    fun doubleBufferReadWrite(blackhole: Blackhole) {
-        val buffer = MutableBuffer.complex(size / 2) { Complex(it.toDouble(), -it.toDouble()) }
+        val buffer = DoubleBuffer(size) { it.toDouble() }
-
+        var res = 0.0
-        (0 until size / 2).forEach {
+        (0 until size).forEach {
-            buffer[it]
+            res += buffer[it]
         }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun bufferViewReadWrite(blackhole: Blackhole) {
+        val buffer = DoubleBuffer(size) { it.toDouble() }.permute(reversedIndices)
+        var res = 0.0
+        (0 until size).forEach {
+            res += buffer[it]
+        }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun bufferViewReadWriteSpecialized(blackhole: Blackhole) {
+        val buffer = DoubleBuffer(size) { it.toDouble() }.permute(reversedIndices)
+        var res = 0.0
+        (0 until size).forEach {
+            res += buffer.getDouble(it)
+        }
+        blackhole.consume(res)
+    }
+
+    @Benchmark
+    fun complexBufferReadWrite(blackhole: Blackhole) = ComplexField {
+        val buffer = Buffer.complex(size / 2) { Complex(it.toDouble(), -it.toDouble()) }
+
+        var res = zero
+        (0 until size / 2).forEach {
+            res += buffer[it]
+        }
+
+        blackhole.consume(res)
     }
 
     private companion object {
         private const val size = 100
+        private val reversedIndices = IntArray(size){it}.apply { reverse() }
     }
 }

examples/src/main/kotlin/space/kscience/kmath/fit/chiSquared.kt

@@ -12,16 +12,16 @@ import space.kscience.kmath.commons.optimization.CMOptimizer
 import space.kscience.kmath.distributions.NormalDistribution
 import space.kscience.kmath.expressions.chiSquaredExpression
 import space.kscience.kmath.expressions.symbol
-import space.kscience.kmath.operations.asIterable
-import space.kscience.kmath.operations.toList
 import space.kscience.kmath.optimization.FunctionOptimizationTarget
 import space.kscience.kmath.optimization.optimizeWith
 import space.kscience.kmath.optimization.resultPoint
 import space.kscience.kmath.optimization.resultValue
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.real.DoubleVector
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.structures.asIterable
+import space.kscience.kmath.structures.toList
 import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
 import space.kscience.plotly.models.TraceValues

examples/src/main/kotlin/space/kscience/kmath/fit/qowFit.kt

@@ -13,15 +13,15 @@ import space.kscience.kmath.distributions.NormalDistribution
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.expressions.binding
 import space.kscience.kmath.expressions.symbol
-import space.kscience.kmath.operations.asIterable
-import space.kscience.kmath.operations.toList
 import space.kscience.kmath.optimization.QowOptimizer
 import space.kscience.kmath.optimization.chiSquaredOrNull
 import space.kscience.kmath.optimization.fitWith
 import space.kscience.kmath.optimization.resultPoint
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.real.map
 import space.kscience.kmath.real.step
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.structures.asIterable
+import space.kscience.kmath.structures.toList
 import space.kscience.plotly.*
 import space.kscience.plotly.models.ScatterMode
 import kotlin.math.abs

examples/src/main/kotlin/space/kscience/kmath/series/analyzeDif.kt

@@ -0,0 +1,47 @@
+package space.kscience.kmath.series
+
+
+import kotlinx.html.FlowContent
+import kotlinx.html.h1
+import space.kscience.kmath.operations.algebra
+import space.kscience.kmath.operations.bufferAlgebra
+import space.kscience.kmath.stat.ksComparisonStatistic
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.slice
+import space.kscience.kmath.structures.toList
+import space.kscience.plotly.*
+import kotlin.math.PI
+
+fun main() = with(Double.algebra.bufferAlgebra.seriesAlgebra()) {
+    fun FlowContent.plotSeries(buffer: Buffer<Double>) {
+        val ls = buffer.labels
+        plot {
+            scatter {
+                x.numbers = ls
+                y.numbers = buffer.toList()
+            }
+            layout {
+                xaxis {
+                    range(0.0..100.0)
+                }
+            }
+        }
+    }
+
+
+    val s1 = series(100) { sin(2 * PI * it / 100) + 1.0 }
+    val s2 = s1.slice(20U..50U).moveTo(40)
+
+    val s3: Buffer<Double> = s1.zip(s2) { l, r -> l + r } //s1 + s2
+    val s4 = ln(s3)
+
+    val kmTest = ksComparisonStatistic(s1, s2)
+
+    Plotly.page {
+        h1 { +"This is my plot" }
+        plotSeries(s1)
+        plotSeries(s2)
+        plotSeries(s4)
+    }.makeFile()
+
+}

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

@@ -10,6 +10,7 @@ import kotlinx.coroutines.async
 import kotlinx.coroutines.runBlocking
 import org.apache.commons.rng.sampling.distribution.BoxMullerNormalizedGaussianSampler
 import org.apache.commons.rng.simple.RandomSource
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.samplers.GaussianSampler
 import java.time.Duration
 import java.time.Instant

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

@@ -9,6 +9,7 @@ import kotlinx.coroutines.runBlocking
 import space.kscience.kmath.chains.Chain
 import space.kscience.kmath.chains.combineWithState
 import space.kscience.kmath.distributions.NormalDistribution
+import space.kscience.kmath.random.RandomGenerator
 
 private data class AveragingChainState(var num: Int = 0, var value: Double = 0.0)
 

examples/src/main/kotlin/space/kscience/kmath/structures/StreamDoubleFieldND.kt

@@ -37,7 +37,7 @@ class StreamDoubleFieldND(override val shape: IntArray) : FieldND<Double, Double
                 this@StreamDoubleFieldND.shape,
                 shape
             )
-            this is BufferND && this.indices == this@StreamDoubleFieldND.strides -> this.buffer as DoubleBuffer
+            this is BufferND && this.shapeIndices == this@StreamDoubleFieldND.strides -> this.buffer as DoubleBuffer
             else -> DoubleBuffer(strides.linearSize) { offset -> get(strides.index(offset)) }
         }
 

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

@@ -7,10 +7,11 @@ package space.kscience.kmath.commons.random
 
 import kotlinx.coroutines.runBlocking
 import space.kscience.kmath.misc.PerformancePitfall
-import space.kscience.kmath.samplers.GaussianSampler
 import space.kscience.kmath.misc.toIntExact
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
-import space.kscience.kmath.stat.next
+import space.kscience.kmath.samplers.GaussianSampler
+import space.kscience.kmath.samplers.next
+
 
 public class CMRandomGeneratorWrapper(
     public val factory: (IntArray) -> RandomGenerator,

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

@@ -13,11 +13,11 @@ import space.kscience.kmath.expressions.Symbol.Companion.x
 import space.kscience.kmath.expressions.Symbol.Companion.y
 import space.kscience.kmath.expressions.chiSquaredExpression
 import space.kscience.kmath.expressions.symbol
-import space.kscience.kmath.operations.map
 import space.kscience.kmath.optimization.*
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.structures.DoubleBuffer
 import space.kscience.kmath.structures.asBuffer
+import space.kscience.kmath.structures.map
 import kotlin.math.pow
 import kotlin.test.Test
 

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

@@ -6,8 +6,8 @@
 package space.kscience.kmath.expressions
 
 import space.kscience.kmath.operations.ExtendedField
-import space.kscience.kmath.operations.asIterable
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.asIterable
 import space.kscience.kmath.structures.indices
 import kotlin.jvm.JvmName
 

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

@@ -15,7 +15,7 @@ import kotlin.math.pow
 import kotlin.math.pow as kpow
 
 public class DoubleBufferND(
-    indexes: ShapeIndexer,
+    indexes: ShapeIndices,
     override val buffer: DoubleBuffer,
 ) : MutableBufferND<Double>(indexes, buffer)
 
@@ -35,7 +35,7 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
         arg: DoubleBufferND,
         transform: (Double) -> Double,
     ): DoubleBufferND {
-        val indexes = arg.indices
+        val indexes = arg.shapeIndices
         val array = arg.buffer.array
         return DoubleBufferND(indexes, DoubleBuffer(indexes.linearSize) { transform(array[it]) })
     }
@@ -45,8 +45,8 @@ public sealed class DoubleFieldOpsND : BufferedFieldOpsND<Double, DoubleField>(D
         r: DoubleBufferND,
         block: (l: Double, r: Double) -> Double,
     ): DoubleBufferND {
-        require(l.indices == r.indices) { "Zip requires the same shapes, but found ${l.shape} on the left and ${r.shape} on the right" }
+        require(l.shapeIndices == r.shapeIndices) { "Zip requires the same shapes, but found ${l.shape} on the left and ${r.shape} on the right" }
-        val indexes = l.indices
+        val indexes = l.shapeIndices
         val lArray = l.buffer.array
         val rArray = r.buffer.array
         return DoubleBufferND(indexes, DoubleBuffer(indexes.linearSize) { block(lArray[it], rArray[it]) })

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

@@ -6,10 +6,10 @@
 package space.kscience.kmath.nd
 
 import space.kscience.kmath.misc.PerformancePitfall
-import space.kscience.kmath.operations.asSequence
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.MutableBuffer
 import space.kscience.kmath.structures.asMutableBuffer
+import space.kscience.kmath.structures.asSequence
 import kotlin.jvm.JvmInline
 
 /**

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

@@ -54,7 +54,7 @@ public interface StructureND<out T> : Featured<StructureFeature>, WithShape {
      * @return the lazy sequence of pairs of indices to values.
      */
     @PerformancePitfall
-    public fun elements(): Sequence<Pair<IntArray, T>> = indices.asSequence().map { it to get(it) }
+    public fun elements(): Sequence<Pair<IntArray, T>> = shapeIndices.asSequence().map { it to get(it) }
 
     /**
      * Feature is some additional structure information that allows to access it special properties or hints.
@@ -71,7 +71,7 @@ public interface StructureND<out T> : Featured<StructureFeature>, WithShape {
             if (st1 === st2) return true
 
             // fast comparison of buffers if possible
-            if (st1 is BufferND && st2 is BufferND && st1.indices == st2.indices)
+            if (st1 is BufferND && st2 is BufferND && st1.shapeIndices == st2.shapeIndices)
                 return Buffer.contentEquals(st1.buffer, st2.buffer)
 
             //element by element comparison if it could not be avoided
@@ -87,7 +87,7 @@ public interface StructureND<out T> : Featured<StructureFeature>, WithShape {
             if (st1 === st2) return true
 
             // fast comparison of buffers if possible
-            if (st1 is BufferND && st2 is BufferND && st1.indices == st2.indices)
+            if (st1 is BufferND && st2 is BufferND && st1.shapeIndices == st2.shapeIndices)
                 return Buffer.contentEquals(st1.buffer, st2.buffer)
 
             //element by element comparison if it could not be avoided

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

@@ -5,7 +5,7 @@
 
 package space.kscience.kmath.structures
 
-import space.kscience.kmath.operations.asSequence
+import space.kscience.kmath.operations.WithSize
 import kotlin.jvm.JvmInline
 import kotlin.reflect.KClass
 
@@ -50,11 +50,11 @@ public fun interface MutableBufferFactory<T> : BufferFactory<T> {
  *
  * @param T the type of elements contained in the buffer.
  */
-public interface Buffer<out T> {
+public interface Buffer<out T> : WithSize {
     /**
      * The size of this buffer.
      */
-    public val size: Int
+    override val size: Int
 
     /**
      * Gets element at given index.
@@ -64,7 +64,7 @@ public interface Buffer<out T> {
     /**
      * Iterates over all elements.
      */
-    public operator fun iterator(): Iterator<T>
+    public operator fun iterator(): Iterator<T> = indices.asSequence().map(::get).iterator()
 
     override fun toString(): String
 
@@ -122,7 +122,14 @@ public interface Buffer<out T> {
 /**
  * Returns an [IntRange] of the valid indices for this [Buffer].
  */
-public val Buffer<*>.indices: IntRange get() = 0 until size
+public val <T> Buffer<T>.indices: IntRange get() = 0 until size
+
+public operator fun <T> Buffer<T>.get(index: UInt): T = get(index.toInt())
+
+/**
+ * if index is in range of buffer, return the value. Otherwise, return null.
+ */
+public fun <T> Buffer<T>.getOrNull(index: Int): T? = if (index in indices) get(index) else null
 
 public fun <T> Buffer<T>.first(): T {
     require(size > 0) { "Can't get the first element of empty buffer" }

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

@@ -0,0 +1,148 @@
+package space.kscience.kmath.structures
+
+import space.kscience.kmath.misc.UnstableKMathAPI
+
+/**
+ * A buffer that wraps an original buffer
+ */
+public interface BufferView<T> : Buffer<T> {
+    public val origin: Buffer<T>
+
+    /**
+     * Get the index in [origin] buffer from index in this buffer.
+     * Return -1 if element not present in the original buffer
+     * This method should be used internally to optimize non-boxing access.
+     */
+    @UnstableKMathAPI
+    public fun originIndex(index: Int): Int
+}
+
+/**
+ * A zero-copy buffer that "sees" only part of original buffer. Slice can't go beyond original buffer borders.
+ */
+public class BufferSlice<T>(
+    override val origin: Buffer<T>,
+    public val offset: UInt = 0U,
+    override val size: Int,
+) : BufferView<T> {
+
+    init {
+        require(size > 0) { "Size must be positive" }
+        require(offset + size.toUInt() <= origin.size.toUInt()) {
+            "End of buffer ${offset + size.toUInt()} is beyond the end of origin buffer size ${origin.size}"
+        }
+    }
+
+    override fun get(index: Int): T = if (index >= size) {
+        throw IndexOutOfBoundsException("$index is out of ${0 until size} rage")
+    } else {
+        origin[index.toUInt() + offset]
+    }
+
+    override fun iterator(): Iterator<T> =
+        (offset until (offset + size.toUInt())).asSequence().map { origin[it] }.iterator()
+
+    @UnstableKMathAPI
+    override fun originIndex(index: Int): Int = if (index >= size) -1 else index - offset.toInt()
+
+    override fun toString(): String = "$origin[$offset..${offset + size.toUInt()}"
+}
+
+/**
+ * An expanded buffer that could include the whole initial buffer ot its part and fills all space beyond it borders with [defaultValue].
+ *
+ * The [offset] parameter shows the shift of expanded buffer start relative to origin start and could be both positive and negative.
+ */
+public class BufferExpanded<T>(
+    override val origin: Buffer<T>,
+    public val defaultValue: T,
+    public val offset: Int = 0,
+    override val size: Int = origin.size,
+) : BufferView<T> {
+
+    init {
+        require(size > 0) { "Size must be positive" }
+    }
+
+    override fun get(index: Int): T = when (index) {
+        !in 0 until size -> throw IndexOutOfBoundsException("Index $index is not in $indices")
+        in -offset until origin.size - offset -> origin[index + offset]
+        else -> defaultValue
+    }
+
+    @UnstableKMathAPI
+    override fun originIndex(index: Int): Int = if (index in -offset until origin.size - offset) index + offset else -1
+
+    override fun toString(): String = "$origin[$offset..${offset + size}]"
+}
+
+/**
+ * Zero-copy select a slice inside the original buffer
+ */
+public fun <T> Buffer<T>.slice(range: UIntRange): BufferView<T> = if (this is BufferSlice) {
+    BufferSlice(
+        origin,
+        this.offset + range.first,
+        (range.last - range.first).toInt() + 1
+    )
+} else {
+    BufferSlice(
+        this,
+        range.first,
+        (range.last - range.first).toInt() + 1
+    )
+}
+
+/**
+ * Resize original buffer to a given range using given [range], filling additional segments with [defaultValue].
+ * Range left border could be negative to designate adding new blank segment to the beginning of the buffer
+ */
+public fun <T> Buffer<T>.expand(
+    range: IntRange,
+    defaultValue: T,
+): BufferView<T> = if (range.first >= 0 && range.last < size) {
+    BufferSlice(
+        this,
+        range.first.toUInt(),
+        (range.last - range.first) + 1
+    )
+} else {
+    BufferExpanded(
+        this,
+        defaultValue,
+        range.first,
+        (range.last - range.first) + 1
+    )
+}
+
+/**
+ * A [BufferView] that overrides indexing of the original buffer
+ */
+public class PermutatedBuffer<T>(
+    override val origin: Buffer<T>,
+    private val permutations: IntArray,
+) : BufferView<T> {
+    init {
+        permutations.forEach { index ->
+            if (index !in origin.indices) {
+                throw IndexOutOfBoundsException("Index $index is not in ${origin.indices}")
+            }
+        }
+    }
+
+    override val size: Int get() = permutations.size
+
+    override fun get(index: Int): T = origin[permutations[index]]
+
+    override fun iterator(): Iterator<T> = permutations.asSequence().map { origin[it] }.iterator()
+
+    @UnstableKMathAPI
+    override fun originIndex(index: Int): Int = if (index in permutations.indices) permutations[index] else -1
+
+    override fun toString(): String = Buffer.toString(this)
+}
+
+/**
+ * Created a permuted view of given buffer using provided [indices]
+ */
+public fun <T> Buffer<T>.permute(indices: IntArray): PermutatedBuffer<T> = PermutatedBuffer(this, indices)

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

@@ -14,7 +14,7 @@ import kotlin.jvm.JvmInline
  * @property array the underlying array.
  */
 @JvmInline
-public value class DoubleBuffer(public val array: DoubleArray) : MutableBuffer<Double> {
+public value class DoubleBuffer(public val array: DoubleArray) : PrimitiveBuffer<Double> {
     override val size: Int get() = array.size
 
     override operator fun get(index: Int): Double = array[index]

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

@@ -14,7 +14,7 @@ import kotlin.jvm.JvmInline
  * @author Iaroslav Postovalov
  */
 @JvmInline
-public value class FloatBuffer(public val array: FloatArray) : MutableBuffer<Float> {
+public value class FloatBuffer(public val array: FloatArray) : PrimitiveBuffer<Float> {
     override val size: Int get() = array.size
 
     override operator fun get(index: Int): Float = array[index]

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

@@ -13,7 +13,7 @@ import kotlin.jvm.JvmInline
  * @property array the underlying array.
  */
 @JvmInline
-public value class IntBuffer(public val array: IntArray) : MutableBuffer<Int> {
+public value class IntBuffer(public val array: IntArray) : PrimitiveBuffer<Int> {
     override val size: Int get() = array.size
 
     override operator fun get(index: Int): Int = array[index]

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

@@ -13,7 +13,7 @@ import kotlin.jvm.JvmInline
  * @property array the underlying array.
  */
 @JvmInline
-public value class LongBuffer(public val array: LongArray) : MutableBuffer<Long> {
+public value class LongBuffer(public val array: LongArray) : PrimitiveBuffer<Long> {
     override val size: Int get() = array.size
 
     override operator fun get(index: Int): Long = array[index]

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

@@ -94,4 +94,7 @@ public interface MutableBuffer<T> : Buffer<T> {
         public inline fun <reified T : Any> auto(size: Int, initializer: (Int) -> T): MutableBuffer<T> =
             auto(T::class, size, initializer)
     }
-}
+}
+
+
+public sealed interface PrimitiveBuffer<T>: MutableBuffer<T>

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

@@ -0,0 +1,38 @@
+package space.kscience.kmath.structures
+
+import space.kscience.kmath.misc.UnstableKMathAPI
+
+/**
+ * Non-boxing access to primitive [Double]
+ */
+
+@UnstableKMathAPI
+public fun Buffer<Double>.getDouble(index: Int): Double = if (this is BufferView) {
+    val originIndex = originIndex(index)
+    if( originIndex>=0) {
+        origin.getDouble(originIndex)
+    } else {
+        get(index)
+    }
+} else if (this is DoubleBuffer) {
+    array[index]
+} else {
+    get(index)
+}
+
+/**
+ * Non-boxing access to primitive [Int]
+ */
+@UnstableKMathAPI
+public fun Buffer<Int>.getInt(index: Int): Int = if (this is BufferView) {
+    val originIndex = originIndex(index)
+    if( originIndex>=0) {
+        origin.getInt(originIndex)
+    } else {
+        get(index)
+    }
+} else if (this is IntBuffer) {
+    array[index]
+} else {
+    get(index)
+}

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

@@ -0,0 +1,27 @@
+package space.kscience.kmath.structures
+
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertFails
+
+internal class BufferExpandedTest {
+    private val buffer = (0..100).toList().asBuffer()
+
+    @Test
+    fun shrink(){
+        val view = buffer.slice(20U..30U)
+        assertEquals(20, view[0])
+        assertEquals(30, view[10])
+        assertFails { view[11] }
+    }
+
+    @Test
+    fun expandNegative(){
+        val view: BufferView<Int> = buffer.expand(-20..113,0)
+        assertEquals(0,view[4])
+        assertEquals(0,view[123])
+        assertEquals(100, view[120])
+        assertFails { view[-2] }
+        assertFails { view[134] }
+    }
+}

kmath-coroutines/src/jvmTest/kotlin/space/kscience/kmath/streaming/RingBufferTest.kt

@@ -7,7 +7,7 @@ package space.kscience.kmath.streaming
 
 import kotlinx.coroutines.flow.*
 import kotlinx.coroutines.runBlocking
-import space.kscience.kmath.operations.asSequence
+import space.kscience.kmath.structures.asSequence
 import kotlin.test.Test
 import kotlin.test.assertEquals
 

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

@@ -13,9 +13,9 @@ import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.algebra
-import space.kscience.kmath.operations.asIterable
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.asIterable
 import kotlin.math.pow
 
 /*

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

@@ -13,8 +13,8 @@ import space.kscience.kmath.structures.DoubleBuffer
  * Map one [BufferND] using function without indices.
  */
 public inline fun BufferND<Double>.mapInline(crossinline transform: DoubleField.(Double) -> Double): BufferND<Double> {
-    val array = DoubleArray(indices.linearSize) { offset -> DoubleField.transform(buffer[offset]) }
+    val array = DoubleArray(shapeIndices.linearSize) { offset -> DoubleField.transform(buffer[offset]) }
-    return BufferND(indices, DoubleBuffer(array))
+    return BufferND(shapeIndices, DoubleBuffer(array))
 }
 
 /**

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

@@ -5,10 +5,10 @@
 
 package space.kscience.kmath.integration
 
-import space.kscience.kmath.operations.map
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
 import space.kscience.kmath.structures.asBuffer
+import space.kscience.kmath.structures.map
 import kotlin.jvm.Synchronized
 import kotlin.math.ulp
 import kotlin.native.concurrent.ThreadLocal

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

@@ -12,10 +12,14 @@ import space.kscience.kmath.interpolation.SplineInterpolator
 import space.kscience.kmath.interpolation.interpolatePolynomials
 import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
-import space.kscience.kmath.operations.*
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.operations.Field
+import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.operations.sum
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.DoubleBuffer
 import space.kscience.kmath.structures.MutableBufferFactory
+import space.kscience.kmath.structures.map
 
 /**
  * Compute analytical indefinite integral of this [PiecewisePolynomial], keeping all intervals intact

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector2DTest.kt

@@ -5,7 +5,8 @@
 
 package space.kscience.kmath.geometry
 
-import space.kscience.kmath.operations.toList
+
+import space.kscience.kmath.structures.toList
 import kotlin.test.Test
 import kotlin.test.assertEquals
 

kmath-geometry/src/commonTest/kotlin/space/kscience/kmath/geometry/Vector3DTest.kt

@@ -5,7 +5,7 @@
 
 package space.kscience.kmath.geometry
 
-import space.kscience.kmath.operations.toList
+import space.kscience.kmath.structures.toList
 import kotlin.test.Test
 import kotlin.test.assertEquals
 

kmath-jupyter/src/main/kotlin/space/kscience/kmath/jupyter/KMathJupyter.kt

@@ -21,9 +21,9 @@ import space.kscience.kmath.expressions.MST
 import space.kscience.kmath.expressions.MstRing
 import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.nd.Structure2D
-import space.kscience.kmath.operations.asSequence
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.asSequence
 
 /**
  * A function for conversion of number to MST for pretty print

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

@@ -6,8 +6,8 @@
 package space.kscience.kmath.distributions
 
 import space.kscience.kmath.chains.Chain
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
-import space.kscience.kmath.stat.Sampler
+import space.kscience.kmath.samplers.Sampler
 
 /**
  * A distribution of typed objects.

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

@@ -7,7 +7,7 @@ package space.kscience.kmath.distributions
 
 import space.kscience.kmath.chains.Chain
 import space.kscience.kmath.chains.SimpleChain
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
 
 /**
  * A multivariate distribution that takes a map of parameters.

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

@@ -6,11 +6,11 @@
 package space.kscience.kmath.distributions
 
 import space.kscience.kmath.chains.Chain
-import space.kscience.kmath.internal.InternalErf
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.samplers.GaussianSampler
+import space.kscience.kmath.samplers.InternalErf
 import space.kscience.kmath.samplers.NormalizedGaussianSampler
 import space.kscience.kmath.samplers.ZigguratNormalizedGaussianSampler
-import space.kscience.kmath.stat.RandomGenerator
 import kotlin.math.*
 
 /**

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/random/MCScope.kt

@@ -3,7 +3,7 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
  */
 
-package space.kscience.kmath.stat
+package space.kscience.kmath.random
 
 import kotlinx.coroutines.*
 import kotlin.contracts.InvocationKind
@@ -22,6 +22,10 @@ public class MCScope(
     public val random: RandomGenerator,
 )
 
+public fun MCScope.asCoroutineScope(): CoroutineScope = object : CoroutineScope {
+    override val coroutineContext: CoroutineContext get() = this@asCoroutineScope.coroutineContext
+}
+
 /**
  * Launches a supervised Monte-Carlo scope
  */

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

@@ -3,7 +3,7 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
  */
 
-package space.kscience.kmath.stat
+package space.kscience.kmath.random
 
 import space.kscience.kmath.chains.BlockingDoubleChain
 import space.kscience.kmath.chains.Chain

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

@@ -3,7 +3,7 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
  */
 
-package space.kscience.kmath.stat
+package space.kscience.kmath.random
 
 import space.kscience.kmath.structures.DoubleBuffer
 import kotlin.random.Random

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

@@ -6,8 +6,7 @@
 package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.BlockingDoubleChain
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
-import space.kscience.kmath.stat.Sampler
 import space.kscience.kmath.structures.DoubleBuffer
 import kotlin.math.ln
 import kotlin.math.pow
@@ -67,7 +66,7 @@ public class AhrensDieterExponentialSampler(public val mean: Double) : Sampler<D
             var qi = 0.0
 
             DoubleArray(16) { i ->
-                qi += ln2.pow(i + 1.0) / space.kscience.kmath.internal.InternalUtils.factorial(i + 1)
+                qi += ln2.pow(i + 1.0) / InternalUtils.factorial(i + 1)
                 qi
             }
         }

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

@@ -6,10 +6,8 @@
 package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.Chain
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
-import space.kscience.kmath.stat.Sampler
+import space.kscience.kmath.random.chain
-import space.kscience.kmath.stat.chain
-import space.kscience.kmath.stat.next
 import kotlin.math.*
 
 /**
@@ -25,14 +23,14 @@ import kotlin.math.*
  */
 public class AhrensDieterMarsagliaTsangGammaSampler private constructor(
     alpha: Double,
-    theta: Double
+    theta: Double,
 ) : Sampler<Double> {
     private val delegate: BaseGammaSampler =
         if (alpha < 1) AhrensDieterGammaSampler(alpha, theta) else MarsagliaTsangGammaSampler(alpha, theta)
 
     private abstract class BaseGammaSampler internal constructor(
         protected val alpha: Double,
-        protected val theta: Double
+        protected val theta: Double,
     ) : Sampler<Double> {
         init {
             require(alpha > 0) { "alpha is not strictly positive: $alpha" }
@@ -119,7 +117,7 @@ public class AhrensDieterMarsagliaTsangGammaSampler private constructor(
     public companion object {
         public fun of(
             alpha: Double,
-            theta: Double
+            theta: Double,
         ): Sampler<Double> = AhrensDieterMarsagliaTsangGammaSampler(alpha, theta)
     }
 }

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

@@ -6,10 +6,8 @@
 package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.Chain
-import space.kscience.kmath.internal.InternalUtils
+import space.kscience.kmath.random.RandomGenerator
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.chain
-import space.kscience.kmath.stat.Sampler
-import space.kscience.kmath.stat.chain
 import kotlin.math.ceil
 import kotlin.math.max
 import kotlin.math.min

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

@@ -6,7 +6,7 @@
 package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.BlockingDoubleChain
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.structures.DoubleBuffer
 import kotlin.math.*
 

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

@@ -7,7 +7,7 @@ package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.BlockingDoubleChain
 import space.kscience.kmath.chains.map
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
 
 /**
  * Sampling from a Gaussian distribution with given mean and standard deviation.
@@ -21,7 +21,7 @@ import space.kscience.kmath.stat.RandomGenerator
 public class GaussianSampler(
     public val mean: Double,
     public val standardDeviation: Double,
-    private val normalized: NormalizedGaussianSampler = BoxMullerSampler
+    private val normalized: NormalizedGaussianSampler = BoxMullerSampler,
 ) : BlockingDoubleSampler {
 
     init {

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

@@ -3,7 +3,7 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
  */
 
-package space.kscience.kmath.internal
+package space.kscience.kmath.samplers
 
 import kotlin.math.abs
 

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

@@ -3,7 +3,7 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
  */
 
-package space.kscience.kmath.internal
+package space.kscience.kmath.samplers
 
 import kotlin.math.*
 

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

@@ -3,7 +3,7 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
  */
 
-package space.kscience.kmath.internal
+package space.kscience.kmath.samplers
 
 import kotlin.math.ln
 import kotlin.math.min

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

@@ -6,8 +6,7 @@
 package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.BlockingIntChain
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
-import space.kscience.kmath.stat.Sampler
 import space.kscience.kmath.structures.IntBuffer
 import kotlin.math.exp
 

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

@@ -6,7 +6,7 @@
 package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.BlockingDoubleChain
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.structures.DoubleBuffer
 import kotlin.math.ln
 import kotlin.math.sqrt

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

@@ -6,10 +6,9 @@
 package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.BlockingDoubleChain
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
-import space.kscience.kmath.stat.Sampler
 
-public interface BlockingDoubleSampler: Sampler<Double>{
+public interface BlockingDoubleSampler : Sampler<Double> {
     override fun sample(generator: RandomGenerator): BlockingDoubleChain
 }
 
@@ -18,6 +17,6 @@ public interface BlockingDoubleSampler: Sampler<Double>{
  * Marker interface for a sampler that generates values from an N(0,1)
  * [Gaussian distribution](https://en.wikipedia.org/wiki/Normal_distribution).
  */
-public fun interface NormalizedGaussianSampler : BlockingDoubleSampler{
+public fun interface NormalizedGaussianSampler : BlockingDoubleSampler {
     public companion object
 }

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

@@ -6,10 +6,8 @@
 package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.BlockingIntChain
-import space.kscience.kmath.internal.InternalUtils
 import space.kscience.kmath.misc.toIntExact
-import space.kscience.kmath.stat.RandomGenerator
+import space.kscience.kmath.random.RandomGenerator
-import space.kscience.kmath.stat.Sampler
 import space.kscience.kmath.structures.IntBuffer
 import kotlin.math.*
 

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

@@ -3,7 +3,7 @@
  * Use of this source code is governed by the Apache 2.0 license that can be found in the license/LICENSE.txt file.
  */
 
-package space.kscience.kmath.stat
+package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.Chain
 import space.kscience.kmath.chains.ConstantChain
@@ -12,6 +12,7 @@ import space.kscience.kmath.chains.zip
 import space.kscience.kmath.operations.Group
 import space.kscience.kmath.operations.ScaleOperations
 import space.kscience.kmath.operations.invoke
+import space.kscience.kmath.random.RandomGenerator
 
 /**
  * Implements [Sampler] by sampling only certain [value].

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

@@ -6,8 +6,7 @@
 package space.kscience.kmath.samplers
 
 import space.kscience.kmath.chains.BlockingDoubleChain
-import space.kscience.kmath.misc.toIntExact
+import space.kscience.kmath.random.RandomGenerator
-import space.kscience.kmath.stat.RandomGenerator
 import space.kscience.kmath.structures.DoubleBuffer
 import kotlin.math.*
 

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/SeriesAlgebra.kt

@@ -0,0 +1,193 @@
+package space.kscience.kmath.series
+
+import space.kscience.kmath.operations.BufferAlgebra
+import space.kscience.kmath.operations.Ring
+import space.kscience.kmath.operations.RingOps
+import space.kscience.kmath.stat.StatisticalAlgebra
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.BufferView
+import kotlin.math.max
+import kotlin.math.min
+
+@PublishedApi
+internal fun IntRange.intersect(other: IntRange): IntRange =
+    max(first, other.first)..min(last, other.last)
+
+@PublishedApi
+internal val IntRange.size: Int
+    get() = last - first + 1
+
+@PublishedApi
+internal operator fun IntRange.contains(other: IntRange): Boolean = (other.first in this) && (other.last in this)
+
+//TODO add permutation sort
+//TODO check rank statistics
+
+
+public interface Series<T> : Buffer<T> {
+    public val origin: Buffer<T>
+
+    /**
+     * Absolute position of start of this [Series] in [SeriesAlgebra]
+     */
+    public val position: Int
+}
+
+public val <T> Series<T>.absoluteIndices: IntRange get() = position until position + size
+
+/**
+ * A [BufferView] with index offset (both positive and negative) and possible size change
+ */
+private class OffsetBufer<T>(
+    override val origin: Buffer<T>,
+    override val position: Int,
+    override val size: Int = origin.size,
+) : Series<T>, Buffer<T> by origin {
+
+    init {
+        require(size > 0) { "Size must be positive" }
+        require(size <= origin.size) { "Slice size is larger than the original buffer" }
+    }
+
+    override fun toString(): String = "$origin-->${position}"
+}
+
+/**
+ * A scope to operation on series
+ */
+public class SeriesAlgebra<T, out A : Ring<T>, out BA : BufferAlgebra<T, A>, L>(
+    override val bufferAlgebra: BA,
+    private val labelResolver: (Int) -> L,
+) : RingOps<Buffer<T>>, StatisticalAlgebra<T, A, BA> {
+
+    public val Buffer<T>.indices: IntRange
+        get() = if (this is Series) {
+            absoluteIndices
+        } else {
+            0 until size
+        }
+
+    /**
+     * Get the value by absolute index in the series algebra or return null if index is out of range
+     */
+    public fun Buffer<T>.getAbsoluteOrNull(index: Int): T? = when {
+        index !in indices -> null
+        this is Series -> origin[index - position]
+        else -> get(index)
+    }
+
+    /**
+     * Get the value by absolute index in the series algebra or throw [IndexOutOfBoundsException] if index is out of range
+     */
+    public fun Buffer<T>.getAbsolute(index: Int): T =
+        getAbsoluteOrNull(index) ?: throw IndexOutOfBoundsException("Index $index is not in $indices")
+
+    /**
+     * Create an offset series with index starting point at [index]
+     */
+    public fun Buffer<T>.moveTo(index: Int): Series<T> = if (this is Series) {
+        OffsetBufer(origin, index, size)
+    } else {
+        OffsetBufer(this, index, size)
+    }
+
+    public val Buffer<T>.offset: Int get() = if (this is Series) position else 0
+
+    /**
+     * Build a new series
+     */
+    public fun series(size: Int, fromIndex: Int = 0, block: A.(label: L) -> T): Series<T> {
+        return bufferFactory(size) {
+            val index = it + fromIndex
+            elementAlgebra.block(labelResolver(index))
+        }.moveTo(fromIndex)
+    }
+
+    /**
+     * Get a label buffer for given buffer.
+     */
+    public val Buffer<T>.labels: List<L> get() = indices.map(labelResolver)
+
+
+    /**
+     * Try to resolve element by label and return null if element with a given label is not found
+     */
+    public operator fun Buffer<T>.get(label: L): T? {
+        val index = labels.indexOf(label)
+        if (index == -1) return null
+        return getAbsolute(index + offset)
+    }
+
+    /**
+     * Map a series to another series of the same size
+     */
+    public inline fun Buffer<T>.map(crossinline transform: A.(T) -> T): Series<T> {
+        val buf = bufferFactory(size) {
+            elementAlgebra.transform(getAbsolute(it))
+        }
+        return buf.moveTo(indices.first)
+    }
+
+    /**
+     * Map series to another series of the same size with label
+     */
+    public inline fun Buffer<T>.mapWithLabel(crossinline transform: A.(arg: T, label: L) -> T): Series<T> {
+        val labels = labels
+        val buf = bufferFactory(size) {
+            elementAlgebra.transform(getAbsolute(it), labels[it])
+        }
+        return buf.moveTo(indices.first)
+    }
+
+    public inline fun <R> Buffer<T>.fold(initial: R, operation: A.(acc: R, T) -> R): R {
+        var accumulator = initial
+        for (index in this.indices) accumulator = elementAlgebra.operation(accumulator, getAbsolute(index))
+        return accumulator
+    }
+
+    public inline fun <R> Buffer<T>.foldWithLabel(initial: R, operation: A.(acc: R, arg: T, label: L) -> R): R {
+        val labels = labels
+        var accumulator = initial
+        for (index in this.indices) accumulator =
+            elementAlgebra.operation(accumulator, getAbsolute(index), labels[index])
+        return accumulator
+    }
+
+    /**
+     * Zip two buffers in the range whe they overlap
+     */
+    public inline fun Buffer<T>.zip(
+        other: Buffer<T>,
+        crossinline operation: A.(left: T, right: T) -> T,
+    ): Series<T> {
+        val newRange = indices.intersect(other.indices)
+        return bufferFactory(newRange.size) {
+            elementAlgebra.operation(
+                getAbsolute(it),
+                other.getAbsolute(it)
+            )
+        }.moveTo(newRange.first)
+    }
+
+    override fun Buffer<T>.unaryMinus(): Buffer<T> = map { -it }
+
+    override fun add(left: Buffer<T>, right: Buffer<T>): Series<T> = left.zip(right) { l, r -> l + r }
+
+    override fun multiply(left: Buffer<T>, right: Buffer<T>): Buffer<T> = left.zip(right) { l, r -> l * r }
+}
+
+public fun <T, A : Ring<T>, BA : BufferAlgebra<T, A>, L> BA.seriesAlgebra(labels: Iterable<L>): SeriesAlgebra<T, A, BA, L> {
+    val l = labels.toList()
+    return SeriesAlgebra(this) {
+        if (it in l.indices) l[it] else error("Index $it is outside of labels range ${l.indices}")
+    }
+}
+
+public fun <T, A : Ring<T>, BA : BufferAlgebra<T, A>, L> BA.seriesAlgebra(labelGenerator: (Int) -> L): SeriesAlgebra<T, A, BA, L> =
+    SeriesAlgebra(this, labelGenerator)
+
+/**
+ * Create a series algebra using offset as a label
+ */
+public fun <T, A : Ring<T>, BA : BufferAlgebra<T, A>> BA.seriesAlgebra(): SeriesAlgebra<T, A, BA, Int> =
+    SeriesAlgebra(this) { it }

kmath-stat/src/commonMain/kotlin/space/kscience/kmath/series/seriesExtensions.kt

@@ -0,0 +1,97 @@
+package space.kscience.kmath.series
+
+import space.kscience.kmath.operations.BufferAlgebra
+import space.kscience.kmath.operations.ExponentialOperations
+import space.kscience.kmath.operations.PowerOperations
+import space.kscience.kmath.operations.TrigonometricOperations
+import space.kscience.kmath.structures.Buffer
+
+
+//trigonometric
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.sin(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : TrigonometricOperations<Buffer<T>> =
+    bufferAlgebra.sin(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.cos(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : TrigonometricOperations<Buffer<T>> =
+    bufferAlgebra.cos(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.tan(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : TrigonometricOperations<Buffer<T>> =
+    bufferAlgebra.tan(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.asin(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : TrigonometricOperations<Buffer<T>> =
+    bufferAlgebra.asin(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.acos(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : TrigonometricOperations<Buffer<T>> =
+    bufferAlgebra.acos(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.atan(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : TrigonometricOperations<Buffer<T>> =
+    bufferAlgebra.atan(arg).moveTo(arg.offset)
+
+
+//exponential
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.exp(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : ExponentialOperations<Buffer<T>> =
+    bufferAlgebra.exp(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.ln(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : ExponentialOperations<Buffer<T>> =
+    bufferAlgebra.ln(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.sinh(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : ExponentialOperations<Buffer<T>> =
+    bufferAlgebra.sinh(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.cosh(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : ExponentialOperations<Buffer<T>> =
+    bufferAlgebra.cosh(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.tanh(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : ExponentialOperations<Buffer<T>> =
+    bufferAlgebra.tanh(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.asinh(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : ExponentialOperations<Buffer<T>> =
+    bufferAlgebra.asinh(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.acosh(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : ExponentialOperations<Buffer<T>> =
+    bufferAlgebra.acosh(arg).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.atanh(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : ExponentialOperations<Buffer<T>> =
+    bufferAlgebra.atanh(arg).moveTo(arg.offset)
+
+
+//power
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.power(
+    arg: Buffer<T>,
+    pow: Number,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : PowerOperations<Buffer<T>> =
+    bufferAlgebra.power(arg, pow).moveTo(arg.offset)
+
+public fun <T, BA> SeriesAlgebra<T, *, BA, *>.sqrt(
+    arg: Buffer<T>,
+): Series<T> where BA : BufferAlgebra<T, *>, BA : PowerOperations<Buffer<T>> =
+    bufferAlgebra.sqrt(arg).moveTo(arg.offset)

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

@@ -45,8 +45,10 @@ public class Mean<T>(
     public companion object {
         @Deprecated("Use Double.mean instead")
         public val double: Mean<Double> = Mean(DoubleField) { sum, count -> sum / count }
+
         @Deprecated("Use Int.mean instead")
         public val int: Mean<Int> = Mean(IntRing) { sum, count -> sum / count }
+
         @Deprecated("Use Long.mean instead")
         public val long: Mean<Long> = Mean(LongRing) { sum, count -> sum / count }
 
@@ -60,6 +62,6 @@ public class Mean<T>(
 //TODO replace with optimized version which respects overflow
 public val Double.Companion.mean: Mean<Double> get() = Mean(DoubleField) { sum, count -> sum / count }
 public val Int.Companion.mean: Mean<Int> get() = Mean(IntRing) { sum, count -> sum / count }
-public val Long.Companion.mean: Mean<Long> get() =  Mean(LongRing) { sum, count -> sum / count }
+public val Long.Companion.mean: Mean<Long> get() = Mean(LongRing) { sum, count -> sum / count }
 
 

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

@@ -5,8 +5,8 @@
 
 package space.kscience.kmath.stat
 
-import space.kscience.kmath.operations.asSequence
 import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.asSequence
 
 /**
  * Non-composable median

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

@@ -0,0 +1,35 @@
+package space.kscience.kmath.stat
+
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.asIterable
+
+/**
+ * Rank statistics
+ */
+public class Rank<T : Comparable<T>> : BlockingStatistic<T, IntArray> {
+    override fun evaluateBlocking(data: Buffer<T>): IntArray = Companion.evaluate(data)
+
+    public companion object {
+        public fun <T : Comparable<T>> evaluate(data: Buffer<T>): IntArray {
+            // https://www.geeksforgeeks.org/rank-elements-array/
+            val permutations = ArrayList<Pair<T, Int>>(data.size)
+            data.asIterable().mapIndexedTo(permutations) { i, v -> v to i }
+            permutations.sortBy { it.first }
+            var rank = 1
+            var i = 0
+            val r = IntArray(data.size)
+            while (i < data.size) {
+                var j = i
+                while (j < data.size - 1 && permutations[j].first == permutations[j + 1]) ++j
+                val n = j - i + 1
+                (0 until n).map { k ->
+                    val idx = permutations[i + k].second
+                    r[idx] = rank + ((n - 1) * 0.5f).toInt()
+                }
+                rank += n
+                i += n
+            }
+            return r
+        }
+    }
+}

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

@@ -0,0 +1,68 @@
+package space.kscience.kmath.stat
+
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.BufferFactory
+import space.kscience.kmath.structures.asIterable
+import space.kscience.kmath.structures.sorted
+
+public interface StatisticalAlgebra<T, out A : Algebra<T>, out BA : BufferAlgebra<T, A>> : Algebra<Buffer<T>> {
+    public val bufferAlgebra: BA
+    public val elementAlgebra: A get() = bufferAlgebra.elementAlgebra
+    public val bufferFactory: BufferFactory<T> get() = bufferAlgebra.bufferFactory
+}
+
+/**
+ * Compute [empirical CDF function](https://en.wikipedia.org/wiki/Empirical_distribution_function)
+ */
+public fun <T : Comparable<T>> StatisticalAlgebra<T, *, *>.ecdf(buffer: Buffer<T>): (T) -> Double = { arg ->
+    buffer.asIterable().count { it < arg }.toDouble() / buffer.size
+}
+
+/**
+ * Resulting value of kolmogorov-smirnov two-sample statistic
+ */
+@UnstableKMathAPI
+public data class KMComparisonResult<T : Comparable<T>>(val n: Int, val m: Int, val value: T)
+
+/**
+ * Kolmogorov-Smirnov sample comparison test
+ * Implementation copied from https://commons.apache.org/proper/commons-math/javadocs/api-3.6.1/index.html?org/apache/commons/math3/stat/inference/KolmogorovSmirnovTest.html
+ */
+@UnstableKMathAPI
+public fun <T : Comparable<T>, A, BA : BufferAlgebra<T, A>> StatisticalAlgebra<T, A, BA>.ksComparisonStatistic(
+    x: Buffer<T>,
+    y: Buffer<T>,
+): KMComparisonResult<T> where A : Group<T>, A : NumericAlgebra<T> = elementAlgebra.invoke {
+    // Copy and sort the sample arrays
+    val sx = x.sorted()
+    val sy = y.sorted()
+    val n = sx.size
+    val m = sy.size
+
+    var rankX = 0
+    var rankY = 0
+    var curD: T = zero
+
+    // Find the max difference between cdf_x and cdf_y
+    var supD: T = zero
+    do {
+        val z = if (sx[rankX] <= sy[rankY]) sx[rankX] else sy[rankY]
+        while (rankX < n && sx[rankX].compareTo(z) == 0) {
+            rankX += 1
+            curD += number(m)
+        }
+
+        while (rankY < m && sy[rankY].compareTo(z) == 0) {
+            rankY += 1
+            curD -= number(n)
+        }
+
+        when {
+            curD > supD -> supD = curD
+            -curD > supD -> supD = -curD
+        }
+    } while (rankX < n && rankY < m)
+    return KMComparisonResult(n, m, supD)
+}

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

@@ -7,13 +7,17 @@ package space.kscience.kmath.stat
 
 import org.apache.commons.rng.UniformRandomProvider
 import org.apache.commons.rng.simple.RandomSource
+import space.kscience.kmath.random.RandomGenerator
 
 /**
  * Implements [RandomGenerator] by delegating all operations to [RandomSource].
  *
  * @property source the underlying [RandomSource] object.
  */
-public class RandomSourceGenerator internal constructor(public val source: RandomSource, seed: Long?) : RandomGenerator {
+public class RandomSourceGenerator internal constructor(
+    public val source: RandomSource,
+    seed: Long?,
+) : RandomGenerator {
     internal val random: UniformRandomProvider = seed?.let { RandomSource.create(source, seed) }
         ?: RandomSource.create(source)
 

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

@@ -8,6 +8,7 @@ package space.kscience.kmath.stat
 import kotlinx.coroutines.runBlocking
 import org.junit.jupiter.api.Assertions
 import org.junit.jupiter.api.Test
+import space.kscience.kmath.random.RandomGenerator
 import space.kscience.kmath.samplers.GaussianSampler
 
 internal class CommonsDistributionsTest {

kmath-stat/src/jvmTest/kotlin/space/kscience/kmath/stat/MCScopeTest.kt

@@ -6,20 +6,22 @@
 package space.kscience.kmath.stat
 
 import kotlinx.coroutines.*
+import space.kscience.kmath.random.launch
+import space.kscience.kmath.random.mcScope
 import java.util.*
 import kotlin.test.Test
 import kotlin.test.assertEquals
 
 data class RandomResult(val branch: String, val order: Int, val value: Int)
 
-typealias ATest = suspend CoroutineScope.() -> Set<RandomResult>
+internal typealias ATest = suspend () -> Set<RandomResult>
 
-class MCScopeTest {
+internal class MCScopeTest {
     val simpleTest: ATest = {
         mcScope(1111) {
             val res = Collections.synchronizedSet(HashSet<RandomResult>())
 
-            launch {
+            launch{
                 //println(random)
                 repeat(10) {
                     delay(10)

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

@@ -6,6 +6,10 @@
 package space.kscience.kmath.stat
 
 import kotlinx.coroutines.runBlocking
+import space.kscience.kmath.random.RandomGenerator
+import space.kscience.kmath.random.chain
+import space.kscience.kmath.samplers.Sampler
+import space.kscience.kmath.samplers.sampleBuffer
 import kotlin.test.Test
 
 class SamplerTest {

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

@@ -9,6 +9,8 @@ import kotlinx.coroutines.flow.first
 import kotlinx.coroutines.flow.last
 import kotlinx.coroutines.flow.take
 import kotlinx.coroutines.runBlocking
+import space.kscience.kmath.random.RandomGenerator
+import space.kscience.kmath.random.chain
 import space.kscience.kmath.streaming.chunked
 import kotlin.test.Test
 import kotlin.test.assertEquals