diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2396c6b27..f4d279362 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -41,6 +41,7 @@
 - Refactor histograms. They are marked as prototype
 - `Complex` and related features moved to a separate module `kmath-complex`
 - Refactor AlgebraElement
+- Add `out` projection to `Buffer` generic
 
 ### Deprecated
 
diff --git a/examples/build.gradle.kts b/examples/build.gradle.kts
index 054945eb7..5439d46dc 100644
--- a/examples/build.gradle.kts
+++ b/examples/build.gradle.kts
@@ -76,6 +76,14 @@ benchmark {
         iterationTimeUnit = "ms" // time unity for iterationTime, default is seconds
         include("DotBenchmark")
     }
+
+    configurations.register("expressions") {
+        warmups = 1 // number of warmup iterations
+        iterations = 3 // number of iterations
+        iterationTime = 500 // time in seconds per iteration
+        iterationTimeUnit = "ms" // time unity for iterationTime, default is seconds
+        include("ExpressionsInterpretersBenchmark")
+    }
 }
 
 kotlin.sourceSets.all {
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
index b54cff926..2ebfee2e7 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/LinearAlgebraBenchmark.kt
@@ -10,7 +10,6 @@ import kscience.kmath.linear.Matrix
 import kscience.kmath.linear.MatrixContext
 import kscience.kmath.linear.inverseWithLup
 import kscience.kmath.linear.real
-import kscience.kmath.operations.invoke
 import org.openjdk.jmh.annotations.Scope
 import org.openjdk.jmh.annotations.State
 import kotlin.random.Random
@@ -34,14 +33,14 @@ internal class LinearAlgebraBenchmark {
 
     @Benchmark
     fun cmLUPInversion() {
-        CMMatrixContext {
+        with(CMMatrixContext) {
             inverse(matrix)
         }
     }
 
     @Benchmark
     fun ejmlInverse() {
-        EjmlMatrixContext {
+        with(EjmlMatrixContext) {
             inverse(matrix)
         }
     }
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt
index 5f7559d02..caebf0497 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/NDFieldBenchmark.kt
@@ -2,7 +2,6 @@ package kscience.kmath.benchmarks
 
 import kscience.kmath.nd.*
 import kscience.kmath.operations.RealField
-import kscience.kmath.operations.invoke
 import kscience.kmath.structures.Buffer
 import org.openjdk.jmh.annotations.Benchmark
 import org.openjdk.jmh.annotations.Scope
@@ -12,7 +11,7 @@ import org.openjdk.jmh.annotations.State
 internal class NDFieldBenchmark {
     @Benchmark
     fun autoFieldAdd() {
-        autoField {
+        with(autoField) {
             var res: NDStructure<Double> = one
             repeat(n) { res += one }
         }
@@ -20,7 +19,7 @@ internal class NDFieldBenchmark {
 
     @Benchmark
     fun specializedFieldAdd() {
-        specializedField {
+        with(specializedField) {
             var res: NDStructure<Double> = one
             repeat(n) { res += 1.0 }
         }
@@ -29,7 +28,7 @@ internal class NDFieldBenchmark {
 
     @Benchmark
     fun boxingFieldAdd() {
-        genericField {
+        with(genericField) {
             var res: NDStructure<Double> = one
             repeat(n) { res += 1.0 }
         }
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt
index 5cf6c2ab8..3f08481b3 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorBenchmark.kt
@@ -2,7 +2,6 @@ package kscience.kmath.benchmarks
 
 import kscience.kmath.nd.*
 import kscience.kmath.operations.RealField
-import kscience.kmath.operations.invoke
 import kscience.kmath.viktor.ViktorNDField
 import org.jetbrains.bio.viktor.F64Array
 import org.openjdk.jmh.annotations.Benchmark
@@ -21,7 +20,7 @@ internal class ViktorBenchmark {
 
     @Benchmark
     fun automaticFieldAddition() {
-        autoField {
+        with(autoField) {
             var res: NDStructure<Double> = one
             repeat(n) { res += 1.0 }
         }
@@ -29,7 +28,7 @@ internal class ViktorBenchmark {
 
     @Benchmark
     fun realFieldAddition() {
-        realField {
+        with(realField) {
             var res: NDStructure<Double> = one
             repeat(n) { res += 1.0 }
         }
@@ -37,7 +36,7 @@ internal class ViktorBenchmark {
 
     @Benchmark
     fun viktorFieldAddition() {
-        viktorField {
+        with(viktorField) {
             var res = one
             repeat(n) { res += 1.0 }
         }
diff --git a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorLogBenchmark.kt b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorLogBenchmark.kt
index 914584b8e..57f556e9e 100644
--- a/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorLogBenchmark.kt
+++ b/examples/src/benchmarks/kotlin/kscience/kmath/benchmarks/ViktorLogBenchmark.kt
@@ -2,7 +2,6 @@ package kscience.kmath.benchmarks
 
 import kscience.kmath.nd.*
 import kscience.kmath.operations.RealField
-import kscience.kmath.operations.invoke
 import kscience.kmath.viktor.ViktorNDField
 import org.jetbrains.bio.viktor.F64Array
 import org.openjdk.jmh.annotations.Benchmark
@@ -22,7 +21,7 @@ internal class ViktorLogBenchmark {
 
     @Benchmark
     fun realFieldLog() {
-        realField {
+        with(realField) {
             val fortyTwo = produce { 42.0 }
             var res = one
             repeat(n) { res = ln(fortyTwo) }
@@ -31,7 +30,7 @@ internal class ViktorLogBenchmark {
 
     @Benchmark
     fun viktorFieldLog() {
-        viktorField {
+        with(viktorField) {
             val fortyTwo = produce { 42.0 }
             var res = one
             repeat(n) { res = ln(fortyTwo) }
diff --git a/kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MstExpression.kt b/kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MstExpression.kt
index 03d33aa2b..8da96f770 100644
--- a/kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MstExpression.kt
+++ b/kmath-ast/src/commonMain/kotlin/kscience/kmath/ast/MstExpression.kt
@@ -21,8 +21,17 @@ public class MstExpression<T, out A : Algebra<T>>(public val algebra: A, public
             null
         } ?: arguments.getValue(StringSymbol(value))
 
-        override fun unaryOperationFunction(operation: String): (arg: T) -> T = algebra.unaryOperationFunction(operation)
-        override fun binaryOperationFunction(operation: String): (left: T, right: T) -> T = algebra.binaryOperationFunction(operation)
+        override fun unaryOperation(operation: String, arg: T): T =
+            algebra.unaryOperation(operation, arg)
+
+        override fun binaryOperation(operation: String, left: T, right: T): T =
+            algebra.binaryOperation(operation, left, right)
+
+        override fun unaryOperationFunction(operation: String): (arg: T) -> T =
+            algebra.unaryOperationFunction(operation)
+
+        override fun binaryOperationFunction(operation: String): (left: T, right: T) -> T =
+            algebra.binaryOperationFunction(operation)
 
         @Suppress("UNCHECKED_CAST")
         override fun number(value: Number): T = if (algebra is NumericAlgebra<*>)
diff --git a/kmath-core/api/kmath-core.api b/kmath-core/api/kmath-core.api
index e58e5d45d..1bfc13c3b 100644
--- a/kmath-core/api/kmath-core.api
+++ b/kmath-core/api/kmath-core.api
@@ -4,7 +4,7 @@ public abstract interface class kscience/kmath/domains/Domain {
 }
 
 public final class kscience/kmath/domains/HyperSquareDomain : kscience/kmath/domains/RealDomain {
-	public synthetic fun <init> ([D[DLkotlin/jvm/internal/DefaultConstructorMarker;)V
+	public fun <init> (Lkscience/kmath/structures/Buffer;Lkscience/kmath/structures/Buffer;)V
 	public fun contains (Lkscience/kmath/structures/Buffer;)Z
 	public fun getDimension ()I
 	public fun getLowerBound (I)Ljava/lang/Double;
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/domains/HyperSquareDomain.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/domains/HyperSquareDomain.kt
index b45cf6bf5..20f19864f 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/domains/HyperSquareDomain.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/domains/HyperSquareDomain.kt
@@ -16,6 +16,7 @@
 package kscience.kmath.domains
 
 import kscience.kmath.linear.Point
+import kscience.kmath.structures.Buffer
 import kscience.kmath.structures.RealBuffer
 import kscience.kmath.structures.indices
 
@@ -25,20 +26,20 @@ import kscience.kmath.structures.indices
  *
  * @author Alexander Nozik
  */
-public class HyperSquareDomain(private val lower: RealBuffer, private val upper: RealBuffer) : RealDomain {
+public class HyperSquareDomain(private val lower: Buffer<Double>, private val upper: Buffer<Double>) : RealDomain {
     public override val dimension: Int get() = lower.size
 
     public override operator fun contains(point: Point<Double>): Boolean = point.indices.all { i ->
         point[i] in lower[i]..upper[i]
     }
 
-    public override fun getLowerBound(num: Int, point: Point<Double>): Double? = lower[num]
+    public override fun getLowerBound(num: Int, point: Point<Double>): Double = lower[num]
 
-    public override fun getLowerBound(num: Int): Double? = lower[num]
+    public override fun getLowerBound(num: Int): Double = lower[num]
 
-    public override fun getUpperBound(num: Int, point: Point<Double>): Double? = upper[num]
+    public override fun getUpperBound(num: Int, point: Point<Double>): Double = upper[num]
 
-    public override fun getUpperBound(num: Int): Double? = upper[num]
+    public override fun getUpperBound(num: Int): Double = upper[num]
 
     public override fun nearestInDomain(point: Point<Double>): Point<Double> {
         val res = DoubleArray(point.size) { i ->
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDStructure.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDStructure.kt
index 4aa1c7d52..c3de02925 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDStructure.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/nd/NDStructure.kt
@@ -78,8 +78,7 @@ public interface NDStructure<T> {
             strides: Strides,
             bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
             initializer: (IntArray) -> T,
-        ): NDBuffer<T> =
-            NDBuffer(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
+        ): NDBuffer<T> = NDBuffer(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
 
         /**
          * Inline create NDStructure with non-boxing buffer implementation if it is possible
@@ -87,15 +86,13 @@ public interface NDStructure<T> {
         public inline fun <reified T : Any> auto(
             strides: Strides,
             crossinline initializer: (IntArray) -> T,
-        ): NDBuffer<T> =
-            NDBuffer(strides, Buffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) })
+        ): NDBuffer<T> = NDBuffer(strides, Buffer.auto(strides.linearSize) { i -> initializer(strides.index(i)) })
 
         public inline fun <T : Any> auto(
             type: KClass<T>,
             strides: Strides,
             crossinline initializer: (IntArray) -> T,
-        ): NDBuffer<T> =
-            NDBuffer(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) })
+        ): NDBuffer<T> = NDBuffer(strides, Buffer.auto(type, strides.linearSize) { i -> initializer(strides.index(i)) })
 
         public fun <T> build(
             shape: IntArray,
@@ -106,8 +103,7 @@ public interface NDStructure<T> {
         public inline fun <reified T : Any> auto(
             shape: IntArray,
             crossinline initializer: (IntArray) -> T,
-        ): NDBuffer<T> =
-            auto(DefaultStrides(shape), initializer)
+        ): NDBuffer<T> = auto(DefaultStrides(shape), initializer)
 
         @JvmName("autoVarArg")
         public inline fun <reified T : Any> auto(
@@ -120,8 +116,7 @@ public interface NDStructure<T> {
             type: KClass<T>,
             vararg shape: Int,
             crossinline initializer: (IntArray) -> T,
-        ): NDBuffer<T> =
-            auto(type, DefaultStrides(shape), initializer)
+        ): NDBuffer<T> = auto(type, DefaultStrides(shape), initializer)
     }
 }
 
diff --git a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffer.kt b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffer.kt
index be48055c7..d214eeba1 100644
--- a/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffer.kt
+++ b/kmath-core/src/commonMain/kotlin/kscience/kmath/structures/Buffer.kt
@@ -21,7 +21,7 @@ public typealias MutableBufferFactory<T> = (Int, (Int) -> T) -> MutableBuffer<T>
  *
  * @param T the type of elements contained in the buffer.
  */
-public interface Buffer<T> {
+public interface Buffer<out T> {
     /**
      * The size of this buffer.
      */
diff --git a/kmath-histograms/build.gradle.kts b/kmath-histograms/build.gradle.kts
index 40196416e..6031b3e85 100644
--- a/kmath-histograms/build.gradle.kts
+++ b/kmath-histograms/build.gradle.kts
@@ -1,4 +1,10 @@
-plugins { id("ru.mipt.npm.mpp") }
+plugins {
+    id("ru.mipt.npm.mpp")
+}
+
+kscience {
+    useAtomic()
+}
 
 kotlin.sourceSets {
     commonMain {
@@ -6,8 +12,8 @@ kotlin.sourceSets {
             api(project(":kmath-core"))
         }
     }
-    commonTest{
-        dependencies{
+    commonTest {
+        dependencies {
             implementation(project(":kmath-for-real"))
         }
     }
diff --git a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/Counters.kt b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/Counters.kt
index 7a263a9fc..bdc03cc1d 100644
--- a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/Counters.kt
+++ b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/Counters.kt
@@ -1,20 +1,49 @@
 package kscience.kmath.histogram
 
-/*
+import kotlinx.atomicfu.atomic
+import kotlinx.atomicfu.getAndUpdate
+import kscience.kmath.operations.RealField
+import kscience.kmath.operations.Space
+
+/**
  * Common representation for atomic counters
- * TODO replace with atomics
  */
-
-public expect class LongCounter() {
-    public fun decrement()
-    public fun increment()
-    public fun reset()
-    public fun sum(): Long
-    public fun add(l: Long)
+public interface Counter<T : Any> {
+    public fun add(delta: T)
+    public val value: T
+    public companion object{
+        public fun real(): ObjectCounter<Double> = ObjectCounter(RealField)
+    }
 }
 
-public expect class DoubleCounter() {
-    public fun reset()
-    public fun sum(): Double
-    public fun add(d: Double)
+public class IntCounter : Counter<Int> {
+    private val innerValue = atomic(0)
+
+    override fun add(delta: Int) {
+        innerValue += delta
+    }
+
+    override val value: Int get() = innerValue.value
 }
+
+public class LongCounter : Counter<Long> {
+    private val innerValue = atomic(0L)
+
+    override fun add(delta: Long) {
+        innerValue += delta
+    }
+
+    override val value: Long get() = innerValue.value
+}
+
+public class ObjectCounter<T : Any>(public val space: Space<T>) : Counter<T> {
+    private val innerValue = atomic(space.zero)
+
+    override fun add(delta: T) {
+        innerValue.getAndUpdate { space.run { it + delta } }
+    }
+
+    override val value: T get() = innerValue.value
+}
+
+
diff --git a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/Histogram.kt b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/Histogram.kt
index 370a01215..91032fc31 100644
--- a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/Histogram.kt
+++ b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/Histogram.kt
@@ -6,18 +6,16 @@ import kscience.kmath.structures.ArrayBuffer
 import kscience.kmath.structures.RealBuffer
 
 /**
- * The bin in the histogram. The histogram is by definition always done in the real space
+ * The binned data element. Could be a histogram bin with a number of counts or an artificial construct
  */
 public interface Bin<T : Any> : Domain<T> {
     /**
      * The value of this bin.
      */
     public val value: Number
-
-    public val center: Point<T>
 }
 
-public interface Histogram<T : Any, out B : Bin<T>> : Iterable<B> {
+public interface Histogram<T : Any, out B : Bin<T>> {
     /**
      * Find existing bin, corresponding to given coordinates
      */
@@ -27,28 +25,31 @@ public interface Histogram<T : Any, out B : Bin<T>> : Iterable<B> {
      * Dimension of the histogram
      */
     public val dimension: Int
+
+    public val bins: Iterable<B>
 }
 
-public interface MutableHistogram<T : Any, out B : Bin<T>> : Histogram<T, B> {
+public fun interface HistogramBuilder<T : Any>  {
 
     /**
      * Increment appropriate bin
      */
-    public fun putWithWeight(point: Point<out T>, weight: Double)
+    public fun putValue(point: Point<out T>, value: Number)
 
-    public fun put(point: Point<out T>): Unit = putWithWeight(point, 1.0)
 }
 
-public fun <T : Any> MutableHistogram<T, *>.put(vararg point: T): Unit = put(ArrayBuffer(point))
+public fun <T : Any, B : Bin<T>> HistogramBuilder<T>.put(point: Point<out T>): Unit = putValue(point, 1.0)
 
-public fun MutableHistogram<Double, *>.put(vararg point: Number): Unit =
+public fun <T : Any> HistogramBuilder<T>.put(vararg point: T): Unit = put(ArrayBuffer(point))
+
+public fun HistogramBuilder<Double>.put(vararg point: Number): Unit =
     put(RealBuffer(point.map { it.toDouble() }.toDoubleArray()))
 
-public fun MutableHistogram<Double, *>.put(vararg point: Double): Unit = put(RealBuffer(point))
-public fun <T : Any> MutableHistogram<T, *>.fill(sequence: Iterable<Point<T>>): Unit = sequence.forEach { put(it) }
+public fun HistogramBuilder<Double>.put(vararg point: Double): Unit = put(RealBuffer(point))
+public fun <T : Any> HistogramBuilder<T>.fill(sequence: Iterable<Point<T>>): Unit = sequence.forEach { put(it) }
 
 /**
  * Pass a sequence builder into histogram
  */
-public fun <T : Any> MutableHistogram<T, *>.fill(block: suspend SequenceScope<Point<T>>.() -> Unit): Unit =
+public fun <T : Any> HistogramBuilder<T>.fill(block: suspend SequenceScope<Point<T>>.() -> Unit): Unit =
     fill(sequence(block).asIterable())
diff --git a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/IndexedHistogramSpace.kt b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/IndexedHistogramSpace.kt
new file mode 100644
index 000000000..53761452f
--- /dev/null
+++ b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/IndexedHistogramSpace.kt
@@ -0,0 +1,76 @@
+package kscience.kmath.histogram
+
+import kscience.kmath.domains.Domain
+import kscience.kmath.linear.Point
+import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.nd.NDSpace
+import kscience.kmath.nd.NDStructure
+import kscience.kmath.nd.Strides
+import kscience.kmath.operations.Space
+import kscience.kmath.operations.SpaceElement
+import kscience.kmath.operations.invoke
+
+/**
+ * A simple histogram bin based on domain
+ */
+public data class DomainBin<T : Comparable<T>>(
+    public val domain: Domain<T>,
+    public override val value: Number,
+) : Bin<T>, Domain<T> by domain
+
+@OptIn(UnstableKMathAPI::class)
+public class IndexedHistogram<T : Comparable<T>, V : Any>(
+    override val context: IndexedHistogramSpace<T, V>,
+    public val values: NDStructure<V>,
+) : Histogram<T, Bin<T>>, SpaceElement<IndexedHistogram<T, V>, IndexedHistogramSpace<T, V>> {
+
+    override fun get(point: Point<T>): Bin<T>? {
+        val index = context.getIndex(point) ?: return null
+        return context.produceBin(index, values[index])
+    }
+
+    override val dimension: Int get() = context.strides.shape.size
+
+    override val bins: Iterable<Bin<T>>
+        get() = context.strides.indices().map {
+            context.produceBin(it, values[it])
+        }.asIterable()
+
+}
+
+/**
+ * A space for producing histograms with values in a NDStructure
+ */
+public interface IndexedHistogramSpace<T : Comparable<T>, V : Any> : Space<IndexedHistogram<T, V>> {
+    //public val valueSpace: Space<V>
+    public val strides: Strides
+    public val histogramValueSpace: NDSpace<V, *> //= NDAlgebra.space(valueSpace, Buffer.Companion::boxing, *shape),
+
+    /**
+     * Resolve index of the bin including given [point]
+     */
+    public fun getIndex(point: Point<T>): IntArray?
+
+    /**
+     * Get a bin domain represented by given index
+     */
+    public fun getDomain(index: IntArray): Domain<T>?
+
+    public fun produceBin(index: IntArray, value: V): Bin<T>
+
+    public fun produce(builder: HistogramBuilder<T>.() -> Unit): IndexedHistogram<T, V>
+
+    override fun add(a: IndexedHistogram<T, V>, b: IndexedHistogram<T, V>): IndexedHistogram<T, V> {
+        require(a.context == this) { "Can't operate on a histogram produced by external space" }
+        require(b.context == this) { "Can't operate on a histogram produced by external space" }
+        return IndexedHistogram(this, histogramValueSpace.invoke { a.values + b.values })
+    }
+
+    override fun multiply(a: IndexedHistogram<T, V>, k: Number): IndexedHistogram<T, V> {
+        require(a.context == this) { "Can't operate on a histogram produced by external space" }
+        return IndexedHistogram(this, histogramValueSpace.invoke { a.values * k })
+    }
+
+    override val zero: IndexedHistogram<T, V> get() = produce { }
+}
+
diff --git a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt
deleted file mode 100644
index 21d873806..000000000
--- a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogram.kt
+++ /dev/null
@@ -1,165 +0,0 @@
-package kscience.kmath.histogram
-
-import kscience.kmath.linear.Point
-import kscience.kmath.nd.DefaultStrides
-import kscience.kmath.nd.NDStructure
-import kscience.kmath.operations.SpaceOperations
-import kscience.kmath.operations.invoke
-import kscience.kmath.structures.*
-import kotlin.math.floor
-
-public data class MultivariateBinDefinition<T : Comparable<T>>(
-    public val space: SpaceOperations<Point<T>>,
-    public val center: Point<T>,
-    public val sizes: Point<T>,
-) {
-    public fun contains(vector: Point<out T>): Boolean {
-        require(vector.size == center.size) { "Dimension mismatch for input vector. Expected ${center.size}, but found ${vector.size}" }
-        val upper = space { center + sizes / 2.0 }
-        val lower = space { center - sizes / 2.0 }
-        return vector.asSequence().mapIndexed { i, value -> value in lower[i]..upper[i] }.all { it }
-    }
-}
-
-
-public class MultivariateBin<T : Comparable<T>>(
-    public val definition: MultivariateBinDefinition<T>,
-    public val count: Long,
-    public override val value: Double,
-) : Bin<T> {
-    public override val dimension: Int
-        get() = definition.center.size
-
-    public override val center: Point<T>
-        get() = definition.center
-
-    public override operator fun contains(point: Point<T>): Boolean = definition.contains(point)
-}
-
-/**
- * Uniform multivariate histogram with fixed borders. Based on NDStructure implementation with complexity of m for bin search, where m is the number of dimensions.
- */
-public class RealHistogram(
-    private val lower: Buffer<Double>,
-    private val upper: Buffer<Double>,
-    private val binNums: IntArray = IntArray(lower.size) { 20 },
-) : MutableHistogram<Double, MultivariateBin<Double>> {
-    private val strides = DefaultStrides(IntArray(binNums.size) { binNums[it] + 2 })
-    private val counts: NDStructure<LongCounter> = NDStructure.auto(strides) { LongCounter() }
-    private val values: NDStructure<DoubleCounter> = NDStructure.auto(strides) { DoubleCounter() }
-    public override val dimension: Int get() = lower.size
-    private val binSize = RealBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] }
-
-    init {
-        // argument checks
-        require(lower.size == upper.size) { "Dimension mismatch in histogram lower and upper limits." }
-        require(lower.size == binNums.size) { "Dimension mismatch in bin count." }
-        require(!(0 until dimension).any { upper[it] - lower[it] < 0 }) { "Range for one of axis is not strictly positive" }
-    }
-
-    /**
-     * Get internal [NDStructure] bin index for given axis
-     */
-    private fun getIndex(axis: Int, value: Double): Int = when {
-        value >= upper[axis] -> binNums[axis] + 1 // overflow
-        value < lower[axis] -> 0 // underflow
-        else -> floor((value - lower[axis]) / binSize[axis]).toInt() + 1
-    }
-
-    private fun getIndex(point: Buffer<out Double>): IntArray = IntArray(dimension) { getIndex(it, point[it]) }
-
-    private fun getCount(index: IntArray): Long = counts[index].sum()
-
-    public fun getCount(point: Buffer<out Double>): Long = getCount(getIndex(point))
-
-    private fun getValue(index: IntArray): Double = values[index].sum()
-
-    public fun getValue(point: Buffer<out Double>): Double = getValue(getIndex(point))
-
-    private fun getBinDefinition(index: IntArray): MultivariateBinDefinition<Double> {
-        val center = index.mapIndexed { axis, i ->
-            when (i) {
-                0 -> Double.NEGATIVE_INFINITY
-                strides.shape[axis] - 1 -> Double.POSITIVE_INFINITY
-                else -> lower[axis] + (i.toDouble() - 0.5) * binSize[axis]
-            }
-        }.asBuffer()
-
-        return MultivariateBinDefinition(RealBufferFieldOperations, center, binSize)
-    }
-
-    public fun getBinDefinition(point: Buffer<out Double>): MultivariateBinDefinition<Double> = getBinDefinition(getIndex(point))
-
-    public override operator fun get(point: Buffer<out Double>): MultivariateBin<Double>? {
-        val index = getIndex(point)
-        return MultivariateBin(getBinDefinition(index), getCount(index),getValue(index))
-    }
-
-//    fun put(point: Point<out Double>){
-//        val index = getIndex(point)
-//        values[index].increment()
-//    }
-
-    public override fun putWithWeight(point: Buffer<out Double>, weight: Double) {
-        val index = getIndex(point)
-        counts[index].increment()
-        values[index].add(weight)
-    }
-
-    public override operator fun iterator(): Iterator<MultivariateBin<Double>> =
-        strides.indices().map { index->
-            MultivariateBin(getBinDefinition(index), counts[index].sum(), values[index].sum())
-        }.iterator()
-
-    /**
-     * NDStructure containing number of events in bins without weights
-     */
-    public fun counts(): NDStructure<Long> = NDStructure.auto(counts.shape) { counts[it].sum() }
-
-    /**
-     * NDStructure containing values of bins including weights
-     */
-    public fun values(): NDStructure<Double> = NDStructure.auto(values.shape) { values[it].sum() }
-
-    public companion object {
-        /**
-         * Use it like
-         * ```
-         *FastHistogram.fromRanges(
-         *  (-1.0..1.0),
-         *  (-1.0..1.0)
-         *)
-         *```
-         */
-        public fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): RealHistogram = RealHistogram(
-            ranges.map(ClosedFloatingPointRange<Double>::start).asBuffer(),
-            ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asBuffer()
-        )
-
-        /**
-         * Use it like
-         * ```
-         *FastHistogram.fromRanges(
-         *  (-1.0..1.0) to 50,
-         *  (-1.0..1.0) to 32
-         *)
-         *```
-         */
-        public fun fromRanges(vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>): RealHistogram =
-            RealHistogram(
-                ListBuffer(
-                    ranges
-                        .map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
-                        .map(ClosedFloatingPointRange<Double>::start)
-                ),
-
-                ListBuffer(
-                    ranges
-                        .map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
-                        .map(ClosedFloatingPointRange<Double>::endInclusive)
-                ),
-
-                ranges.map(Pair<ClosedFloatingPointRange<Double>, Int>::second).toIntArray()
-            )
-    }
-}
diff --git a/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogramSpace.kt b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogramSpace.kt
new file mode 100644
index 000000000..c9beb3a73
--- /dev/null
+++ b/kmath-histograms/src/commonMain/kotlin/kscience/kmath/histogram/RealHistogramSpace.kt
@@ -0,0 +1,121 @@
+package kscience.kmath.histogram
+
+import kscience.kmath.domains.Domain
+import kscience.kmath.domains.HyperSquareDomain
+import kscience.kmath.nd.*
+import kscience.kmath.structures.*
+import kotlin.math.floor
+
+public class RealHistogramSpace(
+    private val lower: Buffer<Double>,
+    private val upper: Buffer<Double>,
+    private val binNums: IntArray = IntArray(lower.size) { 20 },
+) : IndexedHistogramSpace<Double, Double> {
+
+    init {
+        // argument checks
+        require(lower.size == upper.size) { "Dimension mismatch in histogram lower and upper limits." }
+        require(lower.size == binNums.size) { "Dimension mismatch in bin count." }
+        require(!lower.indices.any { upper[it] - lower[it] < 0 }) { "Range for one of axis is not strictly positive" }
+    }
+
+    public val dimension: Int get() = lower.size
+
+    private val shape = IntArray(binNums.size) { binNums[it] + 2 }
+    override val histogramValueSpace: RealNDField = NDAlgebra.real(*shape)
+
+    override val strides: Strides get() = histogramValueSpace.strides
+    private val binSize = RealBuffer(dimension) { (upper[it] - lower[it]) / binNums[it] }
+
+    /**
+     * Get internal [NDStructure] bin index for given axis
+     */
+    private fun getIndex(axis: Int, value: Double): Int = when {
+        value >= upper[axis] -> binNums[axis] + 1 // overflow
+        value < lower[axis] -> 0 // underflow
+        else -> floor((value - lower[axis]) / binSize[axis]).toInt()
+    }
+
+    override fun getIndex(point: Buffer<Double>): IntArray = IntArray(dimension) {
+        getIndex(it, point[it])
+    }
+
+    override fun getDomain(index: IntArray): Domain<Double> {
+        val lowerBoundary = index.mapIndexed { axis, i ->
+            when (i) {
+                0 -> Double.NEGATIVE_INFINITY
+                strides.shape[axis] - 1 -> upper[axis]
+                else -> lower[axis] + (i.toDouble()) * binSize[axis]
+            }
+        }.asBuffer()
+
+        val upperBoundary = index.mapIndexed { axis, i ->
+            when (i) {
+                0 -> lower[axis]
+                strides.shape[axis] - 1 -> Double.POSITIVE_INFINITY
+                else -> lower[axis] + (i.toDouble() + 1) * binSize[axis]
+            }
+        }.asBuffer()
+
+        return HyperSquareDomain(lowerBoundary, upperBoundary)
+    }
+
+
+    override fun produceBin(index: IntArray, value: Double): Bin<Double> {
+        val domain = getDomain(index)
+        return DomainBin(domain, value)
+    }
+
+    override fun produce(builder: HistogramBuilder<Double>.() -> Unit): IndexedHistogram<Double, Double> {
+        val ndCounter = NDStructure.auto(strides) { Counter.real() }
+        val hBuilder = HistogramBuilder<Double> { point, value ->
+            val index = getIndex(point)
+            ndCounter[index].add(1.0)
+        }
+        hBuilder.apply(builder)
+        val values: NDBuffer<Double> = ndCounter.mapToBuffer { it.value }
+        return IndexedHistogram(this, values)
+    }
+
+    public companion object {
+        /**
+         * Use it like
+         * ```
+         *FastHistogram.fromRanges(
+         *  (-1.0..1.0),
+         *  (-1.0..1.0)
+         *)
+         *```
+         */
+        public fun fromRanges(vararg ranges: ClosedFloatingPointRange<Double>): RealHistogramSpace = RealHistogramSpace(
+            ranges.map(ClosedFloatingPointRange<Double>::start).asBuffer(),
+            ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asBuffer()
+        )
+
+        /**
+         * Use it like
+         * ```
+         *FastHistogram.fromRanges(
+         *  (-1.0..1.0) to 50,
+         *  (-1.0..1.0) to 32
+         *)
+         *```
+         */
+        public fun fromRanges(vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>): RealHistogramSpace =
+            RealHistogramSpace(
+                ListBuffer(
+                    ranges
+                        .map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
+                        .map(ClosedFloatingPointRange<Double>::start)
+                ),
+
+                ListBuffer(
+                    ranges
+                        .map(Pair<ClosedFloatingPointRange<Double>, Int>::first)
+                        .map(ClosedFloatingPointRange<Double>::endInclusive)
+                ),
+
+                ranges.map(Pair<ClosedFloatingPointRange<Double>, Int>::second).toIntArray()
+            )
+    }
+}
\ No newline at end of file
diff --git a/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt b/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt
index 87a2b3e68..44b9410d2 100644
--- a/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt
+++ b/kmath-histograms/src/commonTest/kotlin/scietifik/kmath/histogram/MultivariateHistogramTest.kt
@@ -1,8 +1,8 @@
 package scietifik.kmath.histogram
 
-import kscience.kmath.histogram.RealHistogram
-import kscience.kmath.histogram.fill
+import kscience.kmath.histogram.RealHistogramSpace
 import kscience.kmath.histogram.put
+import kscience.kmath.operations.invoke
 import kscience.kmath.real.RealVector
 import kscience.kmath.real.invoke
 import kotlin.random.Random
@@ -11,12 +11,14 @@ import kotlin.test.*
 internal class MultivariateHistogramTest {
     @Test
     fun testSinglePutHistogram() {
-        val histogram = RealHistogram.fromRanges(
+        val hSpace = RealHistogramSpace.fromRanges(
             (-1.0..1.0),
             (-1.0..1.0)
         )
-        histogram.put(0.55, 0.55)
-        val bin = histogram.find { it.value.toInt() > 0 } ?: fail()
+        val histogram = hSpace.produce {
+            put(0.55, 0.55)
+        }
+        val bin = histogram.bins.find { it.value.toInt() > 0 } ?: fail()
         assertTrue { bin.contains(RealVector(0.55, 0.55)) }
         assertTrue { bin.contains(RealVector(0.6, 0.5)) }
         assertFalse { bin.contains(RealVector(-0.55, 0.55)) }
@@ -24,7 +26,7 @@ internal class MultivariateHistogramTest {
 
     @Test
     fun testSequentialPut() {
-        val histogram = RealHistogram.fromRanges(
+        val hSpace = RealHistogramSpace.fromRanges(
             (-1.0..1.0),
             (-1.0..1.0),
             (-1.0..1.0)
@@ -34,12 +36,45 @@ internal class MultivariateHistogramTest {
         fun nextDouble() = random.nextDouble(-1.0, 1.0)
 
         val n = 10000
-
-        histogram.fill {
+        val histogram = hSpace.produce {
             repeat(n) {
-                yield(RealVector(nextDouble(), nextDouble(), nextDouble()))
+                put(nextDouble(), nextDouble(), nextDouble())
             }
         }
-        assertEquals(n, histogram.sumBy { it.value.toInt() })
+        assertEquals(n, histogram.bins.sumBy { it.value.toInt() })
+    }
+
+    @Test
+    fun testHistogramAlgebra() {
+        val hSpace = RealHistogramSpace.fromRanges(
+            (-1.0..1.0),
+            (-1.0..1.0),
+            (-1.0..1.0)
+        ).invoke {
+            val random = Random(1234)
+
+            fun nextDouble() = random.nextDouble(-1.0, 1.0)
+            val n = 10000
+            val histogram1 = produce {
+                repeat(n) {
+                    put(nextDouble(), nextDouble(), nextDouble())
+                }
+            }
+            val histogram2 = produce {
+                repeat(n) {
+                    put(nextDouble(), nextDouble(), nextDouble())
+                }
+            }
+            val res = histogram1 - histogram2
+            assertTrue {
+                strides.indices().all { index ->
+                    res.values[index] <= histogram1.values[index]
+                }
+            }
+            assertTrue {
+                res.bins.count() >= histogram1.bins.count()
+            }
+            assertEquals(0.0, res.bins.sumByDouble { it.value.toDouble() })
+        }
     }
 }
\ No newline at end of file
diff --git a/kmath-histograms/src/jsMain/kotlin/kscience/kmath/histogram/Counters.kt b/kmath-histograms/src/jsMain/kotlin/kscience/kmath/histogram/Counters.kt
deleted file mode 100644
index d0fa1f4c2..000000000
--- a/kmath-histograms/src/jsMain/kotlin/kscience/kmath/histogram/Counters.kt
+++ /dev/null
@@ -1,37 +0,0 @@
-package kscience.kmath.histogram
-
-public actual class LongCounter {
-    private var sum: Long = 0L
-
-    public actual fun decrement() {
-        sum--
-    }
-
-    public actual fun increment() {
-        sum++
-    }
-
-    public actual fun reset() {
-        sum = 0
-    }
-
-    public actual fun sum(): Long = sum
-
-    public actual fun add(l: Long) {
-        sum += l
-    }
-}
-
-public actual class DoubleCounter {
-    private var sum: Double = 0.0
-
-    public actual fun reset() {
-        sum = 0.0
-    }
-
-    public actual fun sum(): Double = sum
-
-    public actual fun add(d: Double) {
-        sum += d
-    }
-}
diff --git a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/Counters.kt b/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/Counters.kt
deleted file mode 100644
index efbd185ef..000000000
--- a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/Counters.kt
+++ /dev/null
@@ -1,7 +0,0 @@
-package kscience.kmath.histogram
-
-import java.util.concurrent.atomic.DoubleAdder
-import java.util.concurrent.atomic.LongAdder
-
-public actual typealias LongCounter = LongAdder
-public actual typealias DoubleCounter = DoubleAdder
diff --git a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/TreeHistogramSpace.kt b/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/TreeHistogramSpace.kt
new file mode 100644
index 000000000..67964fbeb
--- /dev/null
+++ b/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/TreeHistogramSpace.kt
@@ -0,0 +1,155 @@
+package kscience.kmath.histogram
+
+import kscience.kmath.domains.UnivariateDomain
+import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.operations.Space
+import kscience.kmath.structures.Buffer
+import java.util.*
+import kotlin.math.abs
+import kotlin.math.sqrt
+
+private fun <B : ClosedFloatingPointRange<Double>> TreeMap<Double, B>.getBin(value: Double): B? {
+    // check ceiling entry and return it if it is what needed
+    val ceil = ceilingEntry(value)?.value
+    if (ceil != null && value in ceil) return ceil
+    //check floor entry
+    val floor = floorEntry(value)?.value
+    if (floor != null && value in floor) return floor
+    //neither is valid, not found
+    return null
+}
+
+@UnstableKMathAPI
+public class TreeHistogram(
+    override val context: TreeHistogramSpace,
+    private val binMap: TreeMap<Double, out UnivariateBin>,
+) : UnivariateHistogram {
+    override fun get(value: Double): UnivariateBin? = binMap.getBin(value)
+    override val dimension: Int get() = 1
+    override val bins: Collection<UnivariateBin> get() = binMap.values
+}
+
+/**
+ * A space for univariate histograms with variable bin borders based on a tree map
+ */
+@UnstableKMathAPI
+public class TreeHistogramSpace(
+    public val binFactory: (Double) -> UnivariateDomain,
+) : Space<UnivariateHistogram> {
+
+    private class BinCounter(val domain: UnivariateDomain, val counter: Counter<Double> = Counter.real()) :
+        ClosedFloatingPointRange<Double> by domain.range
+
+    public fun produce(builder: UnivariateHistogramBuilder.() -> Unit): UnivariateHistogram {
+        val bins: TreeMap<Double, BinCounter> = TreeMap()
+        val hBuilder = object : UnivariateHistogramBuilder {
+
+            fun get(value: Double): BinCounter? = bins.getBin(value)
+
+            fun createBin(value: Double): BinCounter {
+                val binDefinition = binFactory(value)
+                val newBin = BinCounter(binDefinition)
+                synchronized(this) { bins[binDefinition.center] = newBin }
+                return newBin
+            }
+
+            /**
+             * Thread safe put operation
+             */
+            override fun putValue(at: Double, value: Double) {
+                (get(at) ?: createBin(at)).apply {
+                    counter.add(value)
+                }
+            }
+
+            override fun putValue(point: Buffer<Double>, value: Number) {
+                put(point[0], value.toDouble())
+            }
+        }
+        hBuilder.apply(builder)
+        val resBins = TreeMap<Double, UnivariateBin>()
+        bins.forEach { key, binCounter ->
+            val count = binCounter.counter.value
+            resBins[key] = UnivariateBin(binCounter.domain, count, sqrt(count))
+        }
+        return TreeHistogram(this, resBins)
+    }
+
+    override fun add(
+        a: UnivariateHistogram,
+        b: UnivariateHistogram,
+    ): UnivariateHistogram {
+        require(a.context == this) { "Histogram $a does not belong to this context" }
+        require(b.context == this) { "Histogram $b does not belong to this context" }
+        val bins = TreeMap<Double, UnivariateBin>().apply {
+            (a.bins.map { it.domain } union b.bins.map { it.domain }).forEach { def ->
+                put(def.center,
+                    UnivariateBin(
+                        def,
+                        value = (a[def.center]?.value ?: 0.0) + (b[def.center]?.value ?: 0.0),
+                        standardDeviation = (a[def.center]?.standardDeviation
+                            ?: 0.0) + (b[def.center]?.standardDeviation ?: 0.0)
+                    )
+                )
+            }
+        }
+        return TreeHistogram(this, bins)
+    }
+
+    override fun multiply(a: UnivariateHistogram, k: Number): UnivariateHistogram {
+        val bins = TreeMap<Double, UnivariateBin>().apply {
+            a.bins.forEach { bin ->
+                put(bin.domain.center,
+                    UnivariateBin(
+                        bin.domain,
+                        value = bin.value * k.toDouble(),
+                        standardDeviation = abs(bin.standardDeviation * k.toDouble())
+                    )
+                )
+            }
+        }
+
+        return TreeHistogram(this, bins)
+    }
+
+    override val zero: UnivariateHistogram = produce { }
+
+    public companion object {
+        /**
+         * Build and fill a [UnivariateHistogram]. Returns a read-only histogram.
+         */
+        public fun uniform(
+            binSize: Double,
+            start: Double = 0.0,
+        ): TreeHistogramSpace = TreeHistogramSpace { value ->
+            val center = start + binSize * Math.floor((value - start) / binSize + 0.5)
+            UnivariateDomain((center - binSize / 2)..(center + binSize / 2))
+        }
+
+        /**
+         * Create a histogram with custom cell borders
+         */
+        public fun custom(borders: DoubleArray): TreeHistogramSpace {
+            val sorted = borders.sortedArray()
+
+            return TreeHistogramSpace { value ->
+                when {
+                    value < sorted.first() -> UnivariateDomain(
+                        Double.NEGATIVE_INFINITY..sorted.first()
+                    )
+
+                    value > sorted.last() -> UnivariateDomain(
+                        sorted.last()..Double.POSITIVE_INFINITY
+                    )
+
+                    else -> {
+                        val index = sorted.indices.first { value > sorted[it] }
+                        val left = sorted[index]
+                        val right = sorted[index + 1]
+                        UnivariateDomain(left..right)
+                    }
+                }
+            }
+        }
+    }
+}
\ No newline at end of file
diff --git a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt b/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt
index 10aa9f8ca..68b9019d0 100644
--- a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt
+++ b/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogram.kt
@@ -1,76 +1,38 @@
 package kscience.kmath.histogram
 
-import kscience.kmath.linear.Point
+import kscience.kmath.domains.UnivariateDomain
 import kscience.kmath.misc.UnstableKMathAPI
+import kscience.kmath.operations.Space
 import kscience.kmath.operations.SpaceElement
 import kscience.kmath.structures.Buffer
-import kscience.kmath.structures.asBuffer
 import kscience.kmath.structures.asSequence
-import java.util.*
-import kotlin.math.floor
 
-//TODO move to common
 
-public class UnivariateBin(
-    public val position: Double,
-    public val size: Double,
-) : Bin<Double> {
-    //internal mutation operations
-    internal val counter: LongCounter = LongCounter()
-    internal val weightCounter: DoubleCounter = DoubleCounter()
-
-    /**
-     * The precise number of events ignoring weighting
-     */
-    public val count: Long get() = counter.sum()
-
-    /**
-     * The value of histogram including weighting
-     */
-    public override val value: Double get() = weightCounter.sum()
-
-    public override val center: Point<Double> get() = doubleArrayOf(position).asBuffer()
-    public override val dimension: Int get() = 1
-
-    public operator fun contains(value: Double): Boolean = value in (position - size / 2)..(position + size / 2)
-    public override fun contains(point: Buffer<Double>): Boolean = contains(point[0])
-}
+public val UnivariateDomain.center: Double get() = (range.endInclusive - range.start) / 2
 
 /**
- * Univariate histogram with log(n) bin search speed
+ * A univariate bin based an a range
+ * @param value The value of histogram including weighting
+ * @param standardDeviation Standard deviation of the bin value. Zero or negative if not applicable
  */
-@OptIn(UnstableKMathAPI::class)
-public abstract class UnivariateHistogram protected constructor(
-    protected val bins: TreeMap<Double, UnivariateBin> = TreeMap(),
-) : Histogram<Double, UnivariateBin>, SpaceElement<UnivariateHistogram, UnivariateHistogramSpace> {
-
-    public operator fun get(value: Double): UnivariateBin? {
-        // check ceiling entry and return it if it is what needed
-        val ceil = bins.ceilingEntry(value)?.value
-        if (ceil != null && value in ceil) return ceil
-        //check floor entry
-        val floor = bins.floorEntry(value)?.value
-        if (floor != null && value in floor) return floor
-        //neither is valid, not found
-        return null
-    }
-
-    public override operator fun get(point: Buffer<out Double>): UnivariateBin? = get(point[0])
+public class UnivariateBin(
+    public val domain: UnivariateDomain,
+    override val value: Double,
+    public val standardDeviation: Double,
+) : Bin<Double>, ClosedFloatingPointRange<Double> by domain.range {
 
     public override val dimension: Int get() = 1
 
-    public override operator fun iterator(): Iterator<UnivariateBin> = bins.values.iterator()
+    public override fun contains(point: Buffer<Double>): Boolean = point.size == 1 && contains(point[0])
+}
+
+@OptIn(UnstableKMathAPI::class)
+public interface UnivariateHistogram : Histogram<Double, UnivariateBin>,
+    SpaceElement<UnivariateHistogram, Space<UnivariateHistogram>> {
+    public operator fun get(value: Double): UnivariateBin?
+    public override operator fun get(point: Buffer<Double>): UnivariateBin? = get(point[0])
 
     public companion object {
-        /**
-         * Build a histogram with a uniform binning with a start at [start] and a bin size of [binSize]
-         */
-        public fun uniformBuilder(binSize: Double, start: Double = 0.0): UnivariateHistogramBuilder =
-            UnivariateHistogramSpace { value ->
-                val center = start + binSize * floor((value - start) / binSize + 0.5)
-                UnivariateBin(center, binSize)
-            }.builder()
-
         /**
          * Build and fill a [UnivariateHistogram]. Returns a read-only histogram.
          */
@@ -78,35 +40,7 @@ public abstract class UnivariateHistogram protected constructor(
             binSize: Double,
             start: Double = 0.0,
             builder: UnivariateHistogramBuilder.() -> Unit,
-        ): UnivariateHistogram = uniformBuilder(binSize, start).apply(builder)
-
-        /**
-         * Create a histogram with custom cell borders
-         */
-        public fun customBuilder(borders: DoubleArray): UnivariateHistogramBuilder {
-            val sorted = borders.sortedArray()
-
-            return UnivariateHistogramSpace { value ->
-                when {
-                    value < sorted.first() -> UnivariateBin(
-                        Double.NEGATIVE_INFINITY,
-                        Double.MAX_VALUE
-                    )
-
-                    value > sorted.last() -> UnivariateBin(
-                        Double.POSITIVE_INFINITY,
-                        Double.MAX_VALUE
-                    )
-
-                    else -> {
-                        val index = sorted.indices.first { value > sorted[it] }
-                        val left = sorted[index]
-                        val right = sorted[index + 1]
-                        UnivariateBin((left + right) / 2, (right - left))
-                    }
-                }
-            }.builder()
-        }
+        ): UnivariateHistogram = TreeHistogramSpace.uniform(binSize, start).produce(builder)
 
         /**
          * Build and fill a histogram with custom borders. Returns a read-only histogram.
@@ -114,48 +48,26 @@ public abstract class UnivariateHistogram protected constructor(
         public fun custom(
             borders: DoubleArray,
             builder: UnivariateHistogramBuilder.() -> Unit,
-        ): UnivariateHistogram = customBuilder(borders).apply(builder)
+        ): UnivariateHistogram = TreeHistogramSpace.custom(borders).produce(builder)
+
     }
 }
 
-public class UnivariateHistogramBuilder internal constructor(
-    override val context: UnivariateHistogramSpace,
-) : UnivariateHistogram(), MutableHistogram<Double, UnivariateBin> {
-
-    private fun createBin(value: Double): UnivariateBin = context.binFactory(value).also {
-        synchronized(this) { bins[it.position] = it }
-    }
-
+@UnstableKMathAPI
+public interface UnivariateHistogramBuilder : HistogramBuilder<Double> {
     /**
      * Thread safe put operation
      */
-    public fun put(value: Double, weight: Double = 1.0) {
-        (get(value) ?: createBin(value)).apply {
-            counter.increment()
-            weightCounter.add(weight)
-        }
-    }
+    public fun putValue(at: Double, value: Double = 1.0)
 
-    override fun putWithWeight(point: Buffer<out Double>, weight: Double) {
-        put(point[0], weight)
-    }
-
-    /**
-     * Put several items into a single bin
-     */
-    public fun putMany(value: Double, count: Int, weight: Double = count.toDouble()) {
-        (get(value) ?: createBin(value)).apply {
-            counter.add(count.toLong())
-            weightCounter.add(weight)
-        }
-    }
+    override fun putValue(point: Buffer<Double>, value: Number)
 }
 
 @UnstableKMathAPI
-public fun UnivariateHistogramBuilder.fill(items: Iterable<Double>): Unit = items.forEach(::put)
+public fun UnivariateHistogramBuilder.fill(items: Iterable<Double>): Unit = items.forEach(this::putValue)
 
 @UnstableKMathAPI
-public fun UnivariateHistogramBuilder.fill(array: DoubleArray): Unit = array.forEach(::put)
+public fun UnivariateHistogramBuilder.fill(array: DoubleArray): Unit = array.forEach(this::putValue)
 
 @UnstableKMathAPI
-public fun UnivariateHistogramBuilder.fill(buffer: Buffer<Double>): Unit = buffer.asSequence().forEach(::put)
\ No newline at end of file
+public fun UnivariateHistogramBuilder.fill(buffer: Buffer<Double>): Unit = buffer.asSequence().forEach(this::putValue)
\ No newline at end of file
diff --git a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogramSpace.kt b/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogramSpace.kt
deleted file mode 100644
index 0deeb0a97..000000000
--- a/kmath-histograms/src/jvmMain/kotlin/kscience/kmath/histogram/UnivariateHistogramSpace.kt
+++ /dev/null
@@ -1,25 +0,0 @@
-package kscience.kmath.histogram
-
-import kscience.kmath.operations.Space
-
-public class UnivariateHistogramSpace(public val binFactory: (Double) -> UnivariateBin) : Space<UnivariateHistogram> {
-
-    public fun builder(): UnivariateHistogramBuilder = UnivariateHistogramBuilder(this)
-
-    public fun produce(builder: UnivariateHistogramBuilder.() -> Unit): UnivariateHistogram = builder().apply(builder)
-
-    override fun add(
-        a: UnivariateHistogram,
-        b: UnivariateHistogram,
-    ): UnivariateHistogram {
-        require(a.context == this){"Histogram $a does not belong to this context"}
-        require(b.context == this){"Histogram $b does not belong to this context"}
-        TODO()
-    }
-
-    override fun multiply(a: UnivariateHistogram, k: Number): UnivariateHistogram {
-        TODO("Not yet implemented")
-    }
-
-    override val zero: UnivariateHistogram = produce {  }
-}
\ No newline at end of file