CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ## [Unreleased]
 ### Added
+- Autodiff for generic algebra elements in core!
+- Algebra now has an obligatory `bufferFactory` (#477).
 
 ### Changed
 - Kotlin 1.7

README.md

@@ -44,7 +44,7 @@ module definitions below. The module stability could have the following levels:
 * **PROTOTYPE**. On this level there are no compatibility guarantees. All methods and classes form those modules could
   break any moment. You can still use it, but be sure to fix the specific version.
 * **EXPERIMENTAL**. The general API is decided, but some changes could be made. Volatile API is marked
-  with `@UnstableKmathAPI` or other stability warning annotations.
+  with `@UnstableKMathAPI` or other stability warning annotations.
 * **DEVELOPMENT**. API breaking generally follows semantic versioning ideology. There could be changes in minor
   versions, but not in patch versions. API is protected
   with [binary-compatibility-validator](https://github.com/Kotlin/binary-compatibility-validator) tool.

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

@@ -16,7 +16,6 @@ import space.kscience.kmath.linear.linearSpace
 import space.kscience.kmath.multik.multikAlgebra
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.invoke
-import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.tensorflow.produceWithTF
 import space.kscience.kmath.tensors.core.DoubleTensorAlgebra
 import space.kscience.kmath.tensors.core.tensorAlgebra
@@ -84,7 +83,7 @@ internal class DotBenchmark {
     }
 
     @Benchmark
-    fun bufferedDot(blackhole: Blackhole) = with(DoubleField.linearSpace(Buffer.Companion::auto)) {
+    fun bufferedDot(blackhole: Blackhole) = with(DoubleField.linearSpace) {
         blackhole.consume(matrix1 dot matrix2)
     }
 

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

@@ -20,7 +20,6 @@ import space.kscience.kmath.nd.ndAlgebra
 import space.kscience.kmath.nd.one
 import space.kscience.kmath.nd4j.nd4j
 import space.kscience.kmath.operations.DoubleField
-import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.tensors.core.DoubleTensor
 import space.kscience.kmath.tensors.core.one
 import space.kscience.kmath.tensors.core.tensorAlgebra
@@ -28,12 +27,6 @@ import space.kscience.kmath.viktor.viktorAlgebra
 
 @State(Scope.Benchmark)
 internal class NDFieldBenchmark {
-    @Benchmark
-    fun autoFieldAdd(blackhole: Blackhole) = with(autoField) {
-        var res: StructureND<Double> = one(shape)
-        repeat(n) { res += 1.0 }
-        blackhole.consume(res)
-    }
 
     @Benchmark
     fun specializedFieldAdd(blackhole: Blackhole) = with(specializedField) {
@@ -95,9 +88,8 @@ internal class NDFieldBenchmark {
         private const val dim = 1000
         private const val n = 100
         private val shape = intArrayOf(dim, dim)
-        private val autoField = BufferedFieldOpsND(DoubleField, Buffer.Companion::auto)
         private val specializedField = DoubleField.ndAlgebra
-        private val genericField = BufferedFieldOpsND(DoubleField, Buffer.Companion::boxing)
+        private val genericField = BufferedFieldOpsND(DoubleField)
         private val nd4jField = DoubleField.nd4j
         private val multikField = DoubleField.multikAlgebra
         private val viktorField = DoubleField.viktorAlgebra

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

@@ -10,25 +10,19 @@ import kotlinx.benchmark.Blackhole
 import kotlinx.benchmark.Scope
 import kotlinx.benchmark.State
 import org.jetbrains.bio.viktor.F64Array
-import space.kscience.kmath.nd.*
+import space.kscience.kmath.nd.Shape
+import space.kscience.kmath.nd.StructureND
+import space.kscience.kmath.nd.ndAlgebra
+import space.kscience.kmath.nd.one
 import space.kscience.kmath.operations.DoubleField
-import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.viktor.ViktorFieldND
 
 @State(Scope.Benchmark)
 internal class ViktorBenchmark {
-    @Benchmark
-    fun automaticFieldAddition(blackhole: Blackhole) {
-        with(autoField) {
-            var res: StructureND<Double> = one(shape)
-            repeat(n) { res += 1.0 }
-            blackhole.consume(res)
-        }
-    }
 
     @Benchmark
-    fun realFieldAddition(blackhole: Blackhole) {
+    fun doubleFieldAddition(blackhole: Blackhole) {
-        with(realField) {
+        with(doubleField) {
             var res: StructureND<Double> = one(shape)
             repeat(n) { res += 1.0 }
             blackhole.consume(res)
@@ -58,8 +52,7 @@ internal class ViktorBenchmark {
         private val shape = Shape(dim, dim)
 
         // automatically build context most suited for given type.
-        private val autoField = BufferedFieldOpsND(DoubleField, Buffer.Companion::auto)
+        private val doubleField = DoubleField.ndAlgebra
-        private val realField = DoubleField.ndAlgebra
         private val viktorField = ViktorFieldND(dim, dim)
     }
 }

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

@@ -10,19 +10,17 @@ import kotlinx.benchmark.Blackhole
 import kotlinx.benchmark.Scope
 import kotlinx.benchmark.State
 import org.jetbrains.bio.viktor.F64Array
-import space.kscience.kmath.nd.BufferedFieldOpsND
 import space.kscience.kmath.nd.Shape
 import space.kscience.kmath.nd.ndAlgebra
 import space.kscience.kmath.nd.one
 import space.kscience.kmath.operations.DoubleField
-import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.viktor.ViktorFieldND
 
 @State(Scope.Benchmark)
 internal class ViktorLogBenchmark {
     @Benchmark
     fun realFieldLog(blackhole: Blackhole) {
-        with(realField) {
+        with(doubleField) {
             val fortyTwo = structureND(shape) { 42.0 }
             var res = one(shape)
             repeat(n) { res = ln(fortyTwo) }
@@ -54,8 +52,7 @@ internal class ViktorLogBenchmark {
         private val shape = Shape(dim, dim)
 
         // automatically build context most suited for given type.
-        private val autoField = BufferedFieldOpsND(DoubleField, Buffer.Companion::auto)
+        private val doubleField = DoubleField.ndAlgebra
-        private val realField = DoubleField.ndAlgebra
         private val viktorField = ViktorFieldND(dim, dim)
     }
 }

docs/templates/README-TEMPLATE.md

@@ -44,7 +44,7 @@ module definitions below. The module stability could have the following levels:
 * **PROTOTYPE**. On this level there are no compatibility guarantees. All methods and classes form those modules could
   break any moment. You can still use it, but be sure to fix the specific version.
 * **EXPERIMENTAL**. The general API is decided, but some changes could be made. Volatile API is marked
-  with `@UnstableKmathAPI` or other stability warning annotations.
+  with `@UnstableKMathAPI` or other stability warning annotations.
 * **DEVELOPMENT**. API breaking generally follows semantic versioning ideology. There could be changes in minor
   versions, but not in patch versions. API is protected
   with [binary-compatibility-validator](https://github.com/Kotlin/binary-compatibility-validator) tool.

examples/src/main/kotlin/space/kscience/kmath/operations/complexDemo.kt

@@ -7,7 +7,6 @@ package space.kscience.kmath.operations
 
 import space.kscience.kmath.complex.Complex
 import space.kscience.kmath.complex.algebra
-import space.kscience.kmath.complex.bufferAlgebra
 import space.kscience.kmath.complex.ndAlgebra
 import space.kscience.kmath.nd.BufferND
 import space.kscience.kmath.nd.StructureND
@@ -18,7 +17,7 @@ fun main() = Complex.algebra {
     println(complex * 8 - 5 * i)
 
     //flat buffer
-    val buffer = with(bufferAlgebra){
+    val buffer = with(bufferAlgebra) {
         buffer(8) { Complex(it, -it) }.map { Complex(it.im, it.re) }
     }
     println(buffer)
@@ -30,7 +29,7 @@ fun main() = Complex.algebra {
     println(element)
 
     // 1d element operation
-    val result: StructureND<Complex> = ndAlgebra{
+    val result: StructureND<Complex> = ndAlgebra {
         val a = structureND(8) { (it) -> i * it - it.toDouble() }
         val b = 3
         val c = Complex(1.0, 1.0)

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

@@ -32,12 +32,10 @@ fun main() {
     val shape = Shape(dim, dim)
 
 
-    // automatically build context most suited for given type.
-    val autoField = BufferedFieldOpsND(DoubleField, Buffer.Companion::auto)
     // specialized nd-field for Double. It works as generic Double field as well.
-    val realField = DoubleField.ndAlgebra
+    val doubleField = DoubleField.ndAlgebra
-    //A generic boxing field. It should be used for objects, not primitives.
+    //A generic field. It should be used for objects, not primitives.
-    val boxingField = BufferedFieldOpsND(DoubleField, Buffer.Companion::boxing)
+    val genericField = BufferedFieldOpsND(DoubleField)
     // Nd4j specialized field.
     val nd4jField = DoubleField.nd4j
     //viktor field
@@ -46,14 +44,14 @@ fun main() {
     val parallelField = DoubleField.ndStreaming(dim, dim)
 
     measureAndPrint("Boxing addition") {
-        boxingField {
+        genericField {
             var res: StructureND<Double> = one(shape)
             repeat(n) { res += 1.0 }
         }
     }
 
     measureAndPrint("Specialized addition") {
-        realField {
+        doubleField {
             var res: StructureND<Double> = one(shape)
             repeat(n) { res += 1.0 }
         }
@@ -80,15 +78,8 @@ fun main() {
         }
     }
 
-    measureAndPrint("Automatic field addition") {
-        autoField {
-            var res: StructureND<Double> = one(shape)
-            repeat(n) { res += 1.0 }
-        }
-    }
-
     measureAndPrint("Lazy addition") {
-        val res = realField.one(shape).mapAsync(GlobalScope) {
+        val res = doubleField.one(shape).mapAsync(GlobalScope) {
             var c = 0.0
             repeat(n) {
                 c += 1.0

gradle.properties

@@ -12,4 +12,4 @@ org.gradle.configureondemand=true
 org.gradle.parallel=true
 org.gradle.jvmargs=-Xmx4096m
 
-toolsVersion=0.11.7-kotlin-1.7.0
+toolsVersion=0.11.8-kotlin-1.7.10

kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/Complex.kt

@@ -10,10 +10,7 @@ import space.kscience.kmath.memory.MemorySpec
 import space.kscience.kmath.memory.MemoryWriter
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.*
-import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.*
-import space.kscience.kmath.structures.MemoryBuffer
-import space.kscience.kmath.structures.MutableBuffer
-import space.kscience.kmath.structures.MutableMemoryBuffer
 import kotlin.math.*
 
 /**
@@ -54,6 +51,9 @@ public object ComplexField :
     Norm<Complex, Complex>,
     NumbersAddOps<Complex>,
     ScaleOperations<Complex> {
+    override val bufferFactory: MutableBufferFactory<Complex> = MutableBufferFactory { size, init ->
+        MutableMemoryBuffer.create(Complex, size, init)
+    }
 
     override val zero: Complex = 0.0.toComplex()
     override val one: Complex = 1.0.toComplex()

kmath-complex/src/commonMain/kotlin/space/kscience/kmath/complex/ComplexFieldND.kt

@@ -56,11 +56,6 @@ public sealed class ComplexFieldOpsND : BufferedFieldOpsND<Complex, ComplexField
     public companion object : ComplexFieldOpsND()
 }
 
-@UnstableKMathAPI
-public val ComplexField.bufferAlgebra: BufferFieldOps<Complex, ComplexField>
-    get() = bufferAlgebra(Buffer.Companion::complex)
-
-
 @OptIn(UnstableKMathAPI::class)
 public class ComplexFieldND(override val shape: Shape) :
     ComplexFieldOpsND(), FieldND<Complex, ComplexField>,

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

@@ -0,0 +1,456 @@
+/*
+ * Copyright 2018-2021 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
+ */
+
+package space.kscience.kmath.expressions
+
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.*
+import space.kscience.kmath.structures.Buffer
+import space.kscience.kmath.structures.MutableBuffer
+import space.kscience.kmath.structures.MutableBufferFactory
+import space.kscience.kmath.structures.asBuffer
+import kotlin.math.max
+import kotlin.math.min
+
+/**
+ * Class representing both the value and the differentials of a function.
+ *
+ * This class is the workhorse of the differentiation package.
+ *
+ * This class is an implementation of the extension to Rall's numbers described in Dan Kalman's paper
+ * [Doubly Recursive Multivariate Automatic Differentiation](http://www1.american.edu/cas/mathstat/People/kalman/pdffiles/mmgautodiff.pdf),
+ * Mathematics Magazine, vol. 75, no. 3, June 2002. Rall's numbers are an extension to the real numbers used
+ * throughout mathematical expressions; they hold the derivative together with the value of a function. Dan Kalman's
+ * derivative structures hold all partial derivatives up to any specified order, with respect to any number of free
+ * parameters. Rall's numbers therefore can be seen as derivative structures for order one derivative and one free
+ * parameter, and real numbers can be seen as derivative structures with zero order derivative and no free parameters.
+ *
+ * Derived from
+ * [Commons Math's `DerivativeStructure`](https://github.com/apache/commons-math/blob/924f6c357465b39beb50e3c916d5eb6662194175/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/differentiation/DerivativeStructure.java).
+ */
+@UnstableKMathAPI
+public interface DS<T, A : Ring<T>> {
+    public val derivativeAlgebra: DSAlgebra<T, A>
+    public val data: Buffer<T>
+}
+
+/**
+ * Get a partial derivative.
+ *
+ * @param orders derivation orders with respect to each variable (if all orders are 0, the value is returned).
+ * @return partial derivative.
+ * @see value
+ */
+@UnstableKMathAPI
+private fun <T, A : Ring<T>> DS<T, A>.getPartialDerivative(vararg orders: Int): T =
+    data[derivativeAlgebra.compiler.getPartialDerivativeIndex(*orders)]
+
+/**
+ * Provide a partial derivative with given symbols. On symbol could me mentioned multiple times
+ */
+@UnstableKMathAPI
+public fun <T, A : Ring<T>> DS<T, A>.derivative(symbols: List<Symbol>): T {
+    require(symbols.size <= derivativeAlgebra.order) { "The order of derivative ${symbols.size} exceeds computed order ${derivativeAlgebra.order}" }
+    val ordersCount: Map<String, Int> = symbols.map { it.identity }.groupBy { it }.mapValues { it.value.size }
+    return getPartialDerivative(*symbols.map { ordersCount[it] ?: 0 }.toIntArray())
+}
+
+/**
+ * Provide a partial derivative with given symbols. On symbol could me mentioned multiple times
+ */
+@UnstableKMathAPI
+public fun <T, A : Ring<T>> DS<T, A>.derivative(vararg symbols: Symbol): T {
+    require(symbols.size <= derivativeAlgebra.order) { "The order of derivative ${symbols.size} exceeds computed order ${derivativeAlgebra.order}" }
+    val ordersCount: Map<String, Int> = symbols.map { it.identity }.groupBy { it }.mapValues { it.value.size }
+    return getPartialDerivative(*symbols.map { ordersCount[it] ?: 0 }.toIntArray())
+}
+
+/**
+ * The value part of the derivative structure.
+ *
+ * @see getPartialDerivative
+ */
+@UnstableKMathAPI
+public val <T, A : Ring<T>> DS<T, A>.value: T get() = data[0]
+
+@UnstableKMathAPI
+public abstract class DSAlgebra<T, A : Ring<T>>(
+    public val algebra: A,
+    public val order: Int,
+    bindings: Map<Symbol, T>,
+    public val valueBufferFactory: MutableBufferFactory<T> = algebra.bufferFactory,
+) : ExpressionAlgebra<T, DS<T, A>>, SymbolIndexer {
+
+    /**
+     * Get the compiler for number of free parameters and order.
+     *
+     * @return cached rules set.
+     */
+    @PublishedApi
+    internal val compiler: DSCompiler<T, A> by lazy {
+        val numberOfVariables = bindings.size
+        // get the cached compilers
+        val cache: Array<Array<DSCompiler<T, A>?>>? = null
+
+        // we need to create more compilers
+        val maxParameters: Int = max(numberOfVariables, cache?.size ?: 0)
+        val maxOrder: Int = max(order, if (cache == null) 0 else cache[0].size)
+        val newCache: Array<Array<DSCompiler<T, A>?>> = Array(maxParameters + 1) { arrayOfNulls(maxOrder + 1) }
+
+        if (cache != null) {
+            // preserve the already created compilers
+            for (i in cache.indices) {
+                cache[i].copyInto(newCache[i], endIndex = cache[i].size)
+            }
+        }
+
+        // create the array in increasing diagonal order
+        for (diag in 0..numberOfVariables + order) {
+            for (o in max(0, diag - numberOfVariables)..min(order, diag)) {
+                val p: Int = diag - o
+                if (newCache[p][o] == null) {
+                    val valueCompiler: DSCompiler<T, A>? = if (p == 0) null else newCache[p - 1][o]!!
+                    val derivativeCompiler: DSCompiler<T, A>? = if (o == 0) null else newCache[p][o - 1]!!
+
+                    newCache[p][o] = DSCompiler(
+                        algebra,
+                        valueBufferFactory,
+                        p,
+                        o,
+                        valueCompiler,
+                        derivativeCompiler,
+                    )
+                }
+            }
+        }
+
+        return@lazy newCache[numberOfVariables][order]!!
+    }
+
+    private val variables: Map<Symbol, DSSymbol> by lazy {
+        bindings.entries.mapIndexed { index, (key, value) ->
+            key to DSSymbol(
+                index,
+                key,
+                value,
+            )
+        }.toMap()
+    }
+    override val symbols: List<Symbol> = bindings.map { it.key }
+
+    private fun bufferForVariable(index: Int, value: T): Buffer<T> {
+        val buffer = valueBufferFactory(compiler.size) { algebra.zero }
+        buffer[0] = value
+        if (compiler.order > 0) {
+            // the derivative of the variable with respect to itself is 1.
+
+            val indexOfDerivative = compiler.getPartialDerivativeIndex(*IntArray(symbols.size).apply {
+                set(index, 1)
+            })
+
+            buffer[indexOfDerivative] = algebra.one
+        }
+        return buffer
+    }
+
+    @UnstableKMathAPI
+    private inner class DSImpl(
+        override val data: Buffer<T>,
+    ) : DS<T, A> {
+        override val derivativeAlgebra: DSAlgebra<T, A> get() = this@DSAlgebra
+    }
+
+    protected fun DS(data: Buffer<T>): DS<T, A> = DSImpl(data)
+
+
+    /**
+     * Build an instance representing a variable.
+     *
+     * Instances built using this constructor are considered to be the free variables with respect to which
+     * differentials are computed. As such, their differential with respect to themselves is +1.
+     */
+    public fun variable(
+        index: Int,
+        value: T,
+    ): DS<T, A> {
+        require(index < compiler.freeParameters) { "number is too large: $index >= ${compiler.freeParameters}" }
+        return DS(bufferForVariable(index, value))
+    }
+
+    /**
+     * Build an instance from all its derivatives.
+     *
+     * @param derivatives derivatives sorted according to [DSCompiler.getPartialDerivativeIndex].
+     */
+    public fun ofDerivatives(
+        vararg derivatives: T,
+    ): DS<T, A> {
+        require(derivatives.size == compiler.size) { "dimension mismatch: ${derivatives.size} and ${compiler.size}" }
+        val data = derivatives.asBuffer()
+
+        return DS(data)
+    }
+
+    /**
+     * A class implementing both [DS] and [Symbol].
+     */
+    @UnstableKMathAPI
+    public inner class DSSymbol internal constructor(
+        index: Int,
+        symbol: Symbol,
+        value: T,
+    ) : Symbol by symbol, DS<T, A> {
+        override val derivativeAlgebra: DSAlgebra<T, A> get() = this@DSAlgebra
+        override val data: Buffer<T> = bufferForVariable(index, value)
+    }
+
+    public override fun const(value: T): DS<T, A> {
+        val buffer = valueBufferFactory(compiler.size) { algebra.zero }
+        buffer[0] = value
+
+        return DS(buffer)
+    }
+
+    override fun bindSymbolOrNull(value: String): DSSymbol? = variables[StringSymbol(value)]
+
+    override fun bindSymbol(value: String): DSSymbol =
+        bindSymbolOrNull(value) ?: error("Symbol '$value' is not supported in $this")
+
+    public fun bindSymbolOrNull(symbol: Symbol): DSSymbol? = variables[symbol.identity]
+
+    public fun bindSymbol(symbol: Symbol): DSSymbol =
+        bindSymbolOrNull(symbol.identity) ?: error("Symbol '${symbol}' is not supported in $this")
+
+    public fun DS<T, A>.derivative(symbols: List<Symbol>): T {
+        require(symbols.size <= order) { "The order of derivative ${symbols.size} exceeds computed order $order" }
+        val ordersCount = symbols.groupBy { it }.mapValues { it.value.size }
+        return getPartialDerivative(*variables.keys.map { ordersCount[it] ?: 0 }.toIntArray())
+    }
+
+    public fun DS<T, A>.derivative(vararg symbols: Symbol): T = derivative(symbols.toList())
+
+}
+
+
+/**
+ * A ring over [DS].
+ *
+ * @property order The derivation order.
+ * @param bindings The map of bindings values. All bindings are considered free parameters.
+ */
+@UnstableKMathAPI
+public open class DSRing<T, A>(
+    algebra: A,
+    order: Int,
+    bindings: Map<Symbol, T>,
+    valueBufferFactory: MutableBufferFactory<T>,
+) : DSAlgebra<T, A>(algebra, order, bindings, valueBufferFactory),
+    Ring<DS<T, A>>, ScaleOperations<DS<T, A>>,
+    NumericAlgebra<DS<T, A>>,
+    NumbersAddOps<DS<T, A>> where A : Ring<T>, A : NumericAlgebra<T>, A : ScaleOperations<T> {
+
+    override fun bindSymbolOrNull(value: String): DSSymbol? =
+        super<DSAlgebra>.bindSymbolOrNull(value)
+
+    override fun DS<T, A>.unaryMinus(): DS<T, A> = mapData { -it }
+
+    /**
+     * Create a copy of given [Buffer] and modify it according to [block]
+     */
+    protected inline fun DS<T, A>.transformDataBuffer(block: A.(MutableBuffer<T>) -> Unit): DS<T, A> {
+        require(derivativeAlgebra == this@DSRing) { "All derivative operations should be done in the same algebra" }
+        val newData = valueBufferFactory(compiler.size) { data[it] }
+        algebra.block(newData)
+        return DS(newData)
+    }
+
+    protected fun DS<T, A>.mapData(block: A.(T) -> T): DS<T, A> {
+        require(derivativeAlgebra == this@DSRing) { "All derivative operations should be done in the same algebra" }
+        val newData: Buffer<T> = data.map(valueBufferFactory) {
+            algebra.block(it)
+        }
+        return DS(newData)
+    }
+
+    protected fun DS<T, A>.mapDataIndexed(block: (Int, T) -> T): DS<T, A> {
+        require(derivativeAlgebra == this@DSRing) { "All derivative operations should be done in the same algebra" }
+        val newData: Buffer<T> = data.mapIndexed(valueBufferFactory, block)
+        return DS(newData)
+    }
+
+    override val zero: DS<T, A> by lazy {
+        const(algebra.zero)
+    }
+
+    override val one: DS<T, A> by lazy {
+        const(algebra.one)
+    }
+
+    override fun number(value: Number): DS<T, A> = const(algebra.number(value))
+
+    override fun add(left: DS<T, A>, right: DS<T, A>): DS<T, A> = left.transformDataBuffer { result ->
+        require(right.derivativeAlgebra == this@DSRing) { "All derivative operations should be done in the same algebra" }
+        compiler.add(left.data, 0, right.data, 0, result, 0)
+    }
+
+    override fun scale(a: DS<T, A>, value: Double): DS<T, A> = a.mapData {
+        it.times(value)
+    }
+
+    override fun multiply(
+        left: DS<T, A>,
+        right: DS<T, A>,
+    ): DS<T, A> = left.transformDataBuffer { result ->
+        compiler.multiply(left.data, 0, right.data, 0, result, 0)
+    }
+//
+//    override fun DS<T, A>.minus(arg: DS): DS<T, A> = transformDataBuffer { result ->
+//        subtract(data, 0, arg.data, 0, result, 0)
+//    }
+
+    override operator fun DS<T, A>.plus(other: Number): DS<T, A> = transformDataBuffer {
+        it[0] += number(other)
+    }
+
+//
+//    override operator fun DS<T, A>.minus(other: Number): DS<T, A> =
+//        this + (-other.toDouble())
+
+    override operator fun Number.plus(other: DS<T, A>): DS<T, A> = other + this
+    override operator fun Number.minus(other: DS<T, A>): DS<T, A> = other - this
+}
+
+@UnstableKMathAPI
+public class DerivativeStructureRingExpression<T, A>(
+    public val algebra: A,
+    public val elementBufferFactory: MutableBufferFactory<T> = algebra.bufferFactory,
+    public val function: DSRing<T, A>.() -> DS<T, A>,
+) : DifferentiableExpression<T> where A : Ring<T>, A : ScaleOperations<T>, A : NumericAlgebra<T> {
+    override operator fun invoke(arguments: Map<Symbol, T>): T =
+        DSRing(algebra, 0, arguments, elementBufferFactory).function().value
+
+    override fun derivativeOrNull(symbols: List<Symbol>): Expression<T> = Expression { arguments ->
+        with(
+            DSRing(
+                algebra,
+                symbols.size,
+                arguments,
+                elementBufferFactory
+            )
+        ) { function().derivative(symbols) }
+    }
+}
+
+/**
+ * A field over commons-math [DerivativeStructure].
+ *
+ * @property order The derivation order.
+ * @param bindings The map of bindings values. All bindings are considered free parameters.
+ */
+@UnstableKMathAPI
+public class DSField<T, A : ExtendedField<T>>(
+    algebra: A,
+    order: Int,
+    bindings: Map<Symbol, T>,
+    valueBufferFactory: MutableBufferFactory<T>,
+) : DSRing<T, A>(algebra, order, bindings, valueBufferFactory), ExtendedField<DS<T, A>> {
+    override fun number(value: Number): DS<T, A> = const(algebra.number(value))
+
+    override fun divide(left: DS<T, A>, right: DS<T, A>): DS<T, A> = left.transformDataBuffer { result ->
+        compiler.divide(left.data, 0, right.data, 0, result, 0)
+    }
+
+    override fun sin(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.sin(arg.data, 0, result, 0)
+    }
+
+    override fun cos(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.cos(arg.data, 0, result, 0)
+    }
+
+    override fun tan(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.tan(arg.data, 0, result, 0)
+    }
+
+    override fun asin(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.asin(arg.data, 0, result, 0)
+    }
+
+    override fun acos(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.acos(arg.data, 0, result, 0)
+    }
+
+    override fun atan(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.atan(arg.data, 0, result, 0)
+    }
+
+    override fun sinh(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.sinh(arg.data, 0, result, 0)
+    }
+
+    override fun cosh(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.cosh(arg.data, 0, result, 0)
+    }
+
+    override fun tanh(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.tanh(arg.data, 0, result, 0)
+    }
+
+    override fun asinh(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.asinh(arg.data, 0, result, 0)
+    }
+
+    override fun acosh(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.acosh(arg.data, 0, result, 0)
+    }
+
+    override fun atanh(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.atanh(arg.data, 0, result, 0)
+    }
+
+    override fun power(arg: DS<T, A>, pow: Number): DS<T, A> = when (pow) {
+        is Int -> arg.transformDataBuffer { result ->
+            compiler.pow(arg.data, 0, pow, result, 0)
+        }
+        else -> arg.transformDataBuffer { result ->
+            compiler.pow(arg.data, 0, pow.toDouble(), result, 0)
+        }
+    }
+
+    override fun sqrt(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.sqrt(arg.data, 0, result, 0)
+    }
+
+    public fun power(arg: DS<T, A>, pow: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.pow(arg.data, 0, pow.data, 0, result, 0)
+    }
+
+    override fun exp(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.exp(arg.data, 0, result, 0)
+    }
+
+    override fun ln(arg: DS<T, A>): DS<T, A> = arg.transformDataBuffer { result ->
+        compiler.ln(arg.data, 0, result, 0)
+    }
+}
+
+@UnstableKMathAPI
+public class DSFieldExpression<T, A : ExtendedField<T>>(
+    public val algebra: A,
+    private val valueBufferFactory: MutableBufferFactory<T> = algebra.bufferFactory,
+    public val function: DSField<T, A>.() -> DS<T, A>,
+) : DifferentiableExpression<T> {
+    override operator fun invoke(arguments: Map<Symbol, T>): T =
+        DSField(algebra, 0, arguments, valueBufferFactory).function().value
+
+    override fun derivativeOrNull(symbols: List<Symbol>): Expression<T> = Expression { arguments ->
+        DSField(
+            algebra,
+            symbols.size,
+            arguments,
+            valueBufferFactory,
+        ).run { function().derivative(symbols) }
+    }
+}

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

@@ -11,13 +11,12 @@ import space.kscience.kmath.nd.as2D
 import space.kscience.kmath.nd.asND
 import space.kscience.kmath.operations.*
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.BufferFactory
 import space.kscience.kmath.structures.VirtualBuffer
 import space.kscience.kmath.structures.indices
 
 
 public class BufferedLinearSpace<T, out A : Ring<T>>(
-    private val bufferAlgebra: BufferAlgebra<T, A>
+    private val bufferAlgebra: BufferAlgebra<T, A>,
 ) : LinearSpace<T, A> {
     override val elementAlgebra: A get() = bufferAlgebra.elementAlgebra
 
@@ -91,5 +90,5 @@ public class BufferedLinearSpace<T, out A : Ring<T>>(
 }
 
 
-public fun <T, A : Ring<T>> A.linearSpace(bufferFactory: BufferFactory<T>): BufferedLinearSpace<T, A> =
+public val <T, A : Ring<T>> A.linearSpace: BufferedLinearSpace<T, A>
-    BufferedLinearSpace(BufferRingOps(this, bufferFactory))
+    get() = BufferedLinearSpace(BufferRingOps(this))

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

@@ -11,12 +11,9 @@ import space.kscience.kmath.nd.Structure2D
 import space.kscience.kmath.nd.StructureFeature
 import space.kscience.kmath.nd.as1D
 import space.kscience.kmath.operations.BufferRingOps
-import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.Ring
 import space.kscience.kmath.operations.invoke
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.BufferFactory
-import space.kscience.kmath.structures.DoubleBuffer
 import kotlin.reflect.KClass
 
 /**
@@ -187,18 +184,9 @@ public interface LinearSpace<T, out A : Ring<T>> {
          * A structured matrix with custom buffer
          */
         public fun <T : Any, A : Ring<T>> buffered(
-            algebra: A,
+            algebra: A
-            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
+        ): LinearSpace<T, A> = BufferedLinearSpace(BufferRingOps(algebra))
-        ): LinearSpace<T, A> = BufferedLinearSpace(BufferRingOps(algebra, bufferFactory))
 
-        @Deprecated("use DoubleField.linearSpace")
-        public val double: LinearSpace<Double, DoubleField> = buffered(DoubleField, ::DoubleBuffer)
-
-        /**
-         * Automatic buffered matrix, unboxed if it is possible
-         */
-        public inline fun <reified T : Any, A : Ring<T>> auto(ring: A): LinearSpace<T, A> =
-            buffered(ring, Buffer.Companion::auto)
     }
 }
 

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

@@ -27,5 +27,5 @@ public annotation class UnstableKMathAPI
     RequiresOptIn.Level.WARNING,
 )
 public annotation class PerformancePitfall(
-    val message: String = "Potential performance problem"
+    val message: String = "Potential performance problem",
 )

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

@@ -10,7 +10,6 @@ package space.kscience.kmath.nd
 import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.*
-import space.kscience.kmath.structures.BufferFactory
 
 public interface BufferAlgebraND<T, out A : Algebra<T>> : AlgebraND<T, A> {
     public val indexerBuilder: (IntArray) -> ShapeIndexer
@@ -60,7 +59,7 @@ public inline fun <T, A : Algebra<T>> BufferAlgebraND<T, A>.mapInline(
     return BufferND(
         indexes,
         bufferAlgebra.run {
-            bufferFactory(buffer.size) { elementAlgebra.transform(buffer[it]) }
+            elementBufferFactory(buffer.size) { elementAlgebra.transform(buffer[it]) }
         }
     )
 }
@@ -74,7 +73,7 @@ internal inline fun <T, A : Algebra<T>> BufferAlgebraND<T, A>.mapIndexedInline(
     return BufferND(
         indexes,
         bufferAlgebra.run {
-            bufferFactory(buffer.size) { elementAlgebra.transform(indexes.index(it), buffer[it]) }
+            elementBufferFactory(buffer.size) { elementAlgebra.transform(indexes.index(it), buffer[it]) }
         }
     )
 }
@@ -91,7 +90,7 @@ internal inline fun <T, A : Algebra<T>> BufferAlgebraND<T, A>.zipInline(
     return BufferND(
         indexes,
         bufferAlgebra.run {
-            bufferFactory(lbuffer.size) { elementAlgebra.block(lbuffer[it], rbuffer[it]) }
+            elementBufferFactory(lbuffer.size) { elementAlgebra.block(lbuffer[it], rbuffer[it]) }
         }
     )
 }
@@ -116,9 +115,8 @@ public open class BufferedFieldOpsND<T, out A : Field<T>>(
 
     public constructor(
         elementAlgebra: A,
-        bufferFactory: BufferFactory<T>,
         indexerBuilder: (IntArray) -> ShapeIndexer = BufferAlgebraND.defaultIndexerBuilder,
-    ) : this(BufferFieldOps(elementAlgebra, bufferFactory), indexerBuilder)
+    ) : this(BufferFieldOps(elementAlgebra), indexerBuilder)
 
     @OptIn(PerformancePitfall::class)
     override fun scale(a: StructureND<T>, value: Double): StructureND<T> = a.map { it * value }

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

@@ -69,7 +69,7 @@ public class MutableBufferND<T>(
  * Transform structure to a new structure using provided [MutableBufferFactory] and optimizing if argument is [MutableBufferND]
  */
 public inline fun <T, reified R : Any> MutableStructureND<T>.mapToMutableBuffer(
-    factory: MutableBufferFactory<R> = MutableBuffer.Companion::auto,
+    factory: MutableBufferFactory<R> = MutableBufferFactory(MutableBuffer.Companion::auto),
     crossinline transform: (T) -> R,
 ): MutableBufferND<R> {
     return if (this is MutableBufferND<T>)

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

@@ -120,7 +120,7 @@ public interface StructureND<out T> : Featured<StructureFeature>, WithShape {
          */
         public fun <T> buffered(
             strides: Strides,
-            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
+            bufferFactory: BufferFactory<T> = BufferFactory.boxing(),
             initializer: (IntArray) -> T,
         ): BufferND<T> = BufferND(strides, bufferFactory(strides.linearSize) { i -> initializer(strides.index(i)) })
 
@@ -140,7 +140,7 @@ public interface StructureND<out T> : Featured<StructureFeature>, WithShape {
 
         public fun <T> buffered(
             shape: IntArray,
-            bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
+            bufferFactory: BufferFactory<T> = BufferFactory.boxing(),
             initializer: (IntArray) -> T,
         ): BufferND<T> = buffered(DefaultStrides(shape), bufferFactory, initializer)
 

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

@@ -8,12 +8,7 @@ package space.kscience.kmath.operations
 import space.kscience.kmath.expressions.Symbol
 import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.operations.Ring.Companion.optimizedPower
-
+import space.kscience.kmath.structures.MutableBufferFactory
-/**
- * Stub for DSL the [Algebra] is.
- */
-@DslMarker
-public annotation class KMathContext
 
 /**
  * Represents an algebraic structure.
@@ -21,6 +16,12 @@ public annotation class KMathContext
  * @param T the type of element of this structure.
  */
 public interface Algebra<T> {
+
+    /**
+     * Provide a factory for buffers, associated with this [Algebra]
+     */
+    public val bufferFactory: MutableBufferFactory<T> get() = MutableBufferFactory.boxing()
+
     /**
      * Wraps a raw string to [T] object. This method is designed for three purposes:
      *

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

@@ -10,7 +10,6 @@ import space.kscience.kmath.nd.BufferedRingOpsND
 import space.kscience.kmath.operations.BigInt.Companion.BASE
 import space.kscience.kmath.operations.BigInt.Companion.BASE_SIZE
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.MutableBuffer
 import kotlin.math.log2
 import kotlin.math.max
 import kotlin.math.min
@@ -528,19 +527,11 @@ public fun String.parseBigInteger(): BigInt? {
 
 public val BigInt.algebra: BigIntField get() = BigIntField
 
-@Deprecated("Use BigInt::buffer")
-public inline fun Buffer.Companion.bigInt(size: Int, initializer: (Int) -> BigInt): Buffer<BigInt> =
-    boxing(size, initializer)
-
 public inline fun BigInt.Companion.buffer(size: Int, initializer: (Int) -> BigInt): Buffer<BigInt> =
     Buffer.boxing(size, initializer)
 
-@Deprecated("Use BigInt::mutableBuffer")
+public inline fun BigInt.Companion.mutableBuffer(size: Int, initializer: (Int) -> BigInt): Buffer<BigInt> =
-public inline fun MutableBuffer.Companion.bigInt(size: Int, initializer: (Int) -> BigInt): MutableBuffer<BigInt> =
-    boxing(size, initializer)
-
-public inline fun BigInt.mutableBuffer(size: Int, initializer: (Int) -> BigInt): Buffer<BigInt> =
     Buffer.boxing(size, initializer)
 
 public val BigIntField.nd: BufferedRingOpsND<BigInt, BigIntField>
-    get() = BufferedRingOpsND(BufferRingOps(BigIntField, BigInt::buffer))
+    get() = BufferedRingOpsND(BufferRingOps(BigIntField))

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

@@ -7,8 +7,6 @@ package space.kscience.kmath.operations
 
 import space.kscience.kmath.structures.Buffer
 import space.kscience.kmath.structures.BufferFactory
-import space.kscience.kmath.structures.DoubleBuffer
-import space.kscience.kmath.structures.ShortBuffer
 
 public interface WithSize {
     public val size: Int
@@ -19,11 +17,11 @@ public interface WithSize {
  */
 public interface BufferAlgebra<T, out A : Algebra<T>> : Algebra<Buffer<T>> {
     public val elementAlgebra: A
-    public val bufferFactory: BufferFactory<T>
+    public val elementBufferFactory: BufferFactory<T> get() = elementAlgebra.bufferFactory
 
     public fun buffer(size: Int, vararg elements: T): Buffer<T> {
         require(elements.size == size) { "Expected $size elements but found ${elements.size}" }
-        return bufferFactory(size) { elements[it] }
+        return elementBufferFactory(size) { elements[it] }
     }
 
     //TODO move to multi-receiver inline extension
@@ -36,13 +34,13 @@ public interface BufferAlgebra<T, out A : Algebra<T>> : Algebra<Buffer<T>> {
 
     override fun unaryOperationFunction(operation: String): (arg: Buffer<T>) -> Buffer<T> {
         val operationFunction = elementAlgebra.unaryOperationFunction(operation)
-        return { arg -> bufferFactory(arg.size) { operationFunction(arg[it]) } }
+        return { arg -> elementBufferFactory(arg.size) { operationFunction(arg[it]) } }
     }
 
     override fun binaryOperationFunction(operation: String): (left: Buffer<T>, right: Buffer<T>) -> Buffer<T> {
         val operationFunction = elementAlgebra.binaryOperationFunction(operation)
         return { left, right ->
-            bufferFactory(left.size) { operationFunction(left[it], right[it]) }
+            elementBufferFactory(left.size) { operationFunction(left[it], right[it]) }
         }
     }
 }
@@ -53,7 +51,7 @@ public interface BufferAlgebra<T, out A : Algebra<T>> : Algebra<Buffer<T>> {
 private inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.mapInline(
     buffer: Buffer<T>,
     crossinline block: A.(T) -> T,
-): Buffer<T> = bufferFactory(buffer.size) { elementAlgebra.block(buffer[it]) }
+): Buffer<T> = elementBufferFactory(buffer.size) { elementAlgebra.block(buffer[it]) }
 
 /**
  * Inline map
@@ -61,7 +59,7 @@ private inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.mapInline(
 private inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.mapIndexedInline(
     buffer: Buffer<T>,
     crossinline block: A.(index: Int, arg: T) -> T,
-): Buffer<T> = bufferFactory(buffer.size) { elementAlgebra.block(it, buffer[it]) }
+): Buffer<T> = elementBufferFactory(buffer.size) { elementAlgebra.block(it, buffer[it]) }
 
 /**
  * Inline zip
@@ -72,15 +70,15 @@ private inline fun <T, A : Algebra<T>> BufferAlgebra<T, A>.zipInline(
     crossinline block: A.(l: T, r: T) -> T,
 ): Buffer<T> {
     require(l.size == r.size) { "Incompatible buffer sizes. left: ${l.size}, right: ${r.size}" }
-    return bufferFactory(l.size) { elementAlgebra.block(l[it], r[it]) }
+    return elementBufferFactory(l.size) { elementAlgebra.block(l[it], r[it]) }
 }
 
 public fun <T> BufferAlgebra<T, *>.buffer(size: Int, initializer: (Int) -> T): Buffer<T> {
-    return bufferFactory(size, initializer)
+    return elementBufferFactory(size, initializer)
 }
 
 public fun <T, A> A.buffer(initializer: (Int) -> T): Buffer<T> where A : BufferAlgebra<T, *>, A : WithSize {
-    return bufferFactory(size, initializer)
+    return elementBufferFactory(size, initializer)
 }
 
 public fun <T, A : TrigonometricOperations<T>> BufferAlgebra<T, A>.sin(arg: Buffer<T>): Buffer<T> =
@@ -131,7 +129,6 @@ public fun <T, A : PowerOperations<T>> BufferAlgebra<T, A>.pow(arg: Buffer<T>, p
 
 public open class BufferRingOps<T, A : Ring<T>>(
     override val elementAlgebra: A,
-    override val bufferFactory: BufferFactory<T>,
 ) : BufferAlgebra<T, A>, RingOps<Buffer<T>> {
 
     override fun add(left: Buffer<T>, right: Buffer<T>): Buffer<T> = zipInline(left, right) { l, r -> l + r }
@@ -146,15 +143,13 @@ public open class BufferRingOps<T, A : Ring<T>>(
 }
 
 public val ShortRing.bufferAlgebra: BufferRingOps<Short, ShortRing>
-    get() = BufferRingOps(ShortRing, ::ShortBuffer)
+    get() = BufferRingOps(ShortRing)
 
 public open class BufferFieldOps<T, A : Field<T>>(
     elementAlgebra: A,
-    bufferFactory: BufferFactory<T>,
+) : BufferRingOps<T, A>(elementAlgebra), BufferAlgebra<T, A>, FieldOps<Buffer<T>>, ScaleOperations<Buffer<T>> {
-) : BufferRingOps<T, A>(elementAlgebra, bufferFactory), BufferAlgebra<T, A>, FieldOps<Buffer<T>>,
-    ScaleOperations<Buffer<T>> {
 
-//    override fun add(left: Buffer<T>, right: Buffer<T>): Buffer<T> = zipInline(left, right) { l, r -> l + r }
+    //    override fun add(left: Buffer<T>, right: Buffer<T>): Buffer<T> = zipInline(left, right) { l, r -> l + r }
 //    override fun multiply(left: Buffer<T>, right: Buffer<T>): Buffer<T> = zipInline(left, right) { l, r -> l * r }
     override fun divide(left: Buffer<T>, right: Buffer<T>): Buffer<T> = zipInline(left, right) { l, r -> l / r }
 
@@ -167,30 +162,26 @@ public open class BufferFieldOps<T, A : Field<T>>(
 
 public class BufferField<T, A : Field<T>>(
     elementAlgebra: A,
-    bufferFactory: BufferFactory<T>,
     override val size: Int,
-) : BufferFieldOps<T, A>(elementAlgebra, bufferFactory), Field<Buffer<T>>, WithSize {
+) : BufferFieldOps<T, A>(elementAlgebra), Field<Buffer<T>>, WithSize {
 
-    override val zero: Buffer<T> = bufferFactory(size) { elementAlgebra.zero }
+    override val zero: Buffer<T> = elementAlgebra.bufferFactory(size) { elementAlgebra.zero }
-    override val one: Buffer<T> = bufferFactory(size) { elementAlgebra.one }
+    override val one: Buffer<T> = elementAlgebra.bufferFactory(size) { elementAlgebra.one }
 }
 
 /**
  * Generate full buffer field from given buffer operations
  */
 public fun <T, A : Field<T>> BufferFieldOps<T, A>.withSize(size: Int): BufferField<T, A> =
-    BufferField(elementAlgebra, bufferFactory, size)
+    BufferField(elementAlgebra, size)
 
 //Double buffer specialization
 
 public fun BufferField<Double, *>.buffer(vararg elements: Number): Buffer<Double> {
     require(elements.size == size) { "Expected $size elements but found ${elements.size}" }
-    return bufferFactory(size) { elements[it].toDouble() }
+    return elementBufferFactory(size) { elements[it].toDouble() }
 }
 
-public fun <T, A : Field<T>> A.bufferAlgebra(bufferFactory: BufferFactory<T>): BufferFieldOps<T, A> =
+public val <T, A : Field<T>> A.bufferAlgebra: BufferFieldOps<T, A>
-    BufferFieldOps(this, bufferFactory)
+    get() = BufferFieldOps(this)
-
-public val DoubleField.bufferAlgebra: BufferFieldOps<Double, DoubleField>
-    get() = BufferFieldOps(DoubleField, ::DoubleBuffer)
 

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

@@ -6,12 +6,10 @@
 package space.kscience.kmath.operations
 
 import space.kscience.kmath.linear.Point
-import space.kscience.kmath.misc.UnstableKMathAPI
 import space.kscience.kmath.structures.Buffer
-import space.kscience.kmath.structures.BufferFactory
 import space.kscience.kmath.structures.DoubleBuffer
+import space.kscience.kmath.structures.MutableBufferFactory
 import space.kscience.kmath.structures.asBuffer
-
 import kotlin.math.*
 
 /**
@@ -21,7 +19,7 @@ public abstract class DoubleBufferOps : BufferAlgebra<Double, DoubleField>, Exte
     Norm<Buffer<Double>, Double> {
 
     override val elementAlgebra: DoubleField get() = DoubleField
-    override val bufferFactory: BufferFactory<Double> get() = ::DoubleBuffer
+    override val elementBufferFactory: MutableBufferFactory<Double> get() = elementAlgebra.bufferFactory
 
     override fun Buffer<Double>.map(block: DoubleField.(Double) -> Double): DoubleBuffer =
         mapInline { DoubleField.block(it) }

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

@@ -61,31 +61,39 @@ public inline fun <reified T> Buffer<T>.toTypedArray(): Array<T> = Array(size, :
 /**
  * Create a new buffer from this one with the given mapping function and using [Buffer.Companion.auto] buffer factory.
  */
-public inline fun <T : Any, reified R : Any> Buffer<T>.map(block: (T) -> R): Buffer<R> =
+public inline fun <T, reified R : Any> Buffer<T>.map(block: (T) -> R): Buffer<R> =
     Buffer.auto(size) { block(get(it)) }
 
 /**
  * Create a new buffer from this one with the given mapping function.
  * Provided [bufferFactory] is used to construct the new buffer.
  */
-public inline fun <T : Any, R : Any> Buffer<T>.map(
+public inline fun <T, R> Buffer<T>.map(
     bufferFactory: BufferFactory<R>,
     crossinline block: (T) -> R,
 ): Buffer<R> = bufferFactory(size) { block(get(it)) }
 
 /**
- * Create a new buffer from this one with the given indexed mapping function.
+ * Create a new buffer from this one with the given mapping (indexed) function.
- * Provided [BufferFactory] is used to construct the new buffer.
+ * Provided [bufferFactory] is used to construct the new buffer.
  */
-public inline fun <T : Any, reified R : Any> Buffer<T>.mapIndexed(
+public inline fun <T, R> Buffer<T>.mapIndexed(
-    bufferFactory: BufferFactory<R> = Buffer.Companion::auto,
+    bufferFactory: BufferFactory<R>,
     crossinline block: (index: Int, value: T) -> R,
 ): Buffer<R> = bufferFactory(size) { block(it, get(it)) }
 
+/**
+ * Create a new buffer from this one with the given indexed mapping function.
+ * Provided [BufferFactory] is used to construct the new buffer.
+ */
+public inline fun <T, reified R : Any> Buffer<T>.mapIndexed(
+    crossinline block: (index: Int, value: T) -> R,
+): Buffer<R> = Buffer.auto(size) { block(it, get(it)) }
+
 /**
  * Fold given buffer according to [operation]
  */
-public inline fun <T : Any, R> Buffer<T>.fold(initial: R, operation: (acc: R, T) -> R): R {
+public inline fun <T, R> Buffer<T>.fold(initial: R, operation: (acc: R, T) -> R): R {
     var accumulator = initial
     for (index in this.indices) accumulator = operation(accumulator, get(index))
     return accumulator
@@ -104,9 +112,9 @@ public inline fun <T : Any, R> Buffer<T>.foldIndexed(initial: R, operation: (ind
  * Zip two buffers using given [transform].
  */
 @UnstableKMathAPI
-public inline fun <T1 : Any, T2 : Any, reified R : Any> Buffer<T1>.zip(
+public inline fun <T1, T2 : Any, reified R : Any> Buffer<T1>.zip(
     other: Buffer<T2>,
-    bufferFactory: BufferFactory<R> = Buffer.Companion::auto,
+    bufferFactory: BufferFactory<R> = BufferFactory.auto(),
     crossinline transform: (T1, T2) -> R,
 ): Buffer<R> {
     require(size == other.size) { "Buffer size mismatch in zip: expected $size but found ${other.size}" }

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

@@ -5,6 +5,7 @@
 
 package space.kscience.kmath.operations
 
+import space.kscience.kmath.structures.*
 import kotlin.math.pow as kpow
 
 /**
@@ -65,6 +66,8 @@ public interface ExtendedField<T> : ExtendedFieldOps<T>, Field<T>, PowerOperatio
  */
 @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
 public object DoubleField : ExtendedField<Double>, Norm<Double, Double>, ScaleOperations<Double> {
+    override val bufferFactory: MutableBufferFactory<Double> = MutableBufferFactory(::DoubleBuffer)
+
     override inline val zero: Double get() = 0.0
     override inline val one: Double get() = 1.0
 
@@ -123,6 +126,8 @@ public val Double.Companion.algebra: DoubleField get() = DoubleField
  */
 @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
 public object FloatField : ExtendedField<Float>, Norm<Float, Float> {
+    override val bufferFactory: MutableBufferFactory<Float> = MutableBufferFactory(::FloatBuffer)
+
     override inline val zero: Float get() = 0.0f
     override inline val one: Float get() = 1.0f
 
@@ -177,11 +182,10 @@ public val Float.Companion.algebra: FloatField get() = FloatField
  */
 @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
 public object IntRing : Ring<Int>, Norm<Int, Int>, NumericAlgebra<Int> {
-    override inline val zero: Int
+    override val bufferFactory: MutableBufferFactory<Int> = MutableBufferFactory(::IntBuffer)
-        get() = 0
 
-    override inline val one: Int
+    override inline val zero: Int get() = 0
-        get() = 1
+    override inline val one: Int get() = 1
 
     override fun number(value: Number): Int = value.toInt()
     override inline fun add(left: Int, right: Int): Int = left + right
@@ -201,11 +205,10 @@ public val Int.Companion.algebra: IntRing get() = IntRing
  */
 @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
 public object ShortRing : Ring<Short>, Norm<Short, Short>, NumericAlgebra<Short> {
-    override inline val zero: Short
+    override val bufferFactory: MutableBufferFactory<Short> = MutableBufferFactory(::ShortBuffer)
-        get() = 0
 
-    override inline val one: Short
+    override inline val zero: Short get() = 0
-        get() = 1
+    override inline val one: Short get() = 1
 
     override fun number(value: Number): Short = value.toShort()
     override inline fun add(left: Short, right: Short): Short = (left + right).toShort()
@@ -225,11 +228,10 @@ public val Short.Companion.algebra: ShortRing get() = ShortRing
  */
 @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
 public object ByteRing : Ring<Byte>, Norm<Byte, Byte>, NumericAlgebra<Byte> {
-    override inline val zero: Byte
+    override val bufferFactory: MutableBufferFactory<Byte> = MutableBufferFactory(::ByteBuffer)
-        get() = 0
 
-    override inline val one: Byte
+    override inline val zero: Byte get() = 0
-        get() = 1
+    override inline val one: Byte get() = 1
 
     override fun number(value: Number): Byte = value.toByte()
     override inline fun add(left: Byte, right: Byte): Byte = (left + right).toByte()
@@ -249,11 +251,10 @@ public val Byte.Companion.algebra: ByteRing get() = ByteRing
  */
 @Suppress("EXTENSION_SHADOWED_BY_MEMBER", "OVERRIDE_BY_INLINE", "NOTHING_TO_INLINE")
 public object LongRing : Ring<Long>, Norm<Long, Long>, NumericAlgebra<Long> {
-    override inline val zero: Long
+    override val bufferFactory: MutableBufferFactory<Long> = MutableBufferFactory(::LongBuffer)
-        get() = 0L
 
-    override inline val one: Long
+    override inline val zero: Long get() = 0L
-        get() = 1L
+    override inline val one: Long get() = 1L
 
     override fun number(value: Number): Long = value.toLong()
     override inline fun add(left: Long, right: Long): Long = left + right

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

@@ -14,14 +14,34 @@ import kotlin.reflect.KClass
  *
  * @param T the type of buffer.
  */
-public typealias BufferFactory<T> = (Int, (Int) -> T) -> Buffer<T>
+public fun interface BufferFactory<T> {
+    public operator fun invoke(size: Int, builder: (Int) -> T): Buffer<T>
+
+    public companion object{
+        public inline fun <reified T : Any> auto(): BufferFactory<T> =
+            BufferFactory(Buffer.Companion::auto)
+
+        public fun <T> boxing(): BufferFactory<T> =
+            BufferFactory(Buffer.Companion::boxing)
+    }
+}
 
 /**
  * Function that produces [MutableBuffer] from its size and function that supplies values.
  *
  * @param T the type of buffer.
  */
-public typealias MutableBufferFactory<T> = (Int, (Int) -> T) -> MutableBuffer<T>
+public fun interface MutableBufferFactory<T> : BufferFactory<T> {
+    override fun invoke(size: Int, builder: (Int) -> T): MutableBuffer<T>
+
+    public companion object {
+        public inline fun <reified T : Any> auto(): MutableBufferFactory<T> =
+            MutableBufferFactory(MutableBuffer.Companion::auto)
+
+        public fun <T> boxing(): MutableBufferFactory<T> =
+            MutableBufferFactory(MutableBuffer.Companion::boxing)
+    }
+}
 
 /**
  * A generic read-only random-access structure for both primitives and objects.

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

@@ -0,0 +1,57 @@
+/*
+ * Copyright 2018-2021 KMath contributors.
+ * 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.structures
+
+import kotlin.jvm.JvmInline
+
+/**
+ * Specialized [MutableBuffer] implementation over [ByteArray].
+ *
+ * @property array the underlying array.
+ */
+@JvmInline
+public value class ByteBuffer(public val array: ByteArray) : MutableBuffer<Byte> {
+    override val size: Int get() = array.size
+
+    override operator fun get(index: Int): Byte = array[index]
+
+    override operator fun set(index: Int, value: Byte) {
+        array[index] = value
+    }
+
+    override operator fun iterator(): ByteIterator = array.iterator()
+    override fun copy(): MutableBuffer<Byte> = ByteBuffer(array.copyOf())
+}
+
+/**
+ * Creates a new [ByteBuffer] with the specified [size], where each element is calculated by calling the specified
+ * [init] function.
+ *
+ * The function [init] is called for each array element sequentially starting from the first one.
+ * It should return the value for a buffer element given its index.
+ */
+public inline fun ByteBuffer(size: Int, init: (Int) -> Byte): ByteBuffer = ByteBuffer(ByteArray(size) { init(it) })
+
+/**
+ * Returns a new [ByteBuffer] of given elements.
+ */
+public fun ByteBuffer(vararg bytes: Byte): ByteBuffer = ByteBuffer(bytes)
+
+/**
+ * Returns a new [ByteArray] containing all the elements of this [Buffer].
+ */
+public fun Buffer<Byte>.toByteArray(): ByteArray = when (this) {
+    is ByteBuffer -> array.copyOf()
+    else -> ByteArray(size, ::get)
+}
+
+/**
+ * Returns [ByteBuffer] over this array.
+ *
+ * @receiver the array.
+ * @return the new buffer.
+ */
+public fun ByteArray.asBuffer(): ByteBuffer = ByteBuffer(this)

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

@@ -0,0 +1,59 @@
+/*
+ * Copyright 2018-2021 KMath contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the LICENSE file.
+ */
+
+@file:OptIn(UnstableKMathAPI::class)
+
+package space.kscience.kmath.expressions
+
+import space.kscience.kmath.misc.UnstableKMathAPI
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.structures.DoubleBuffer
+import kotlin.contracts.InvocationKind
+import kotlin.contracts.contract
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertFails
+
+internal inline fun diff(
+    order: Int,
+    vararg parameters: Pair<Symbol, Double>,
+    block: DSField<Double, DoubleField>.() -> Unit,
+) {
+    contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) }
+    DSField(DoubleField, order, mapOf(*parameters), ::DoubleBuffer).block()
+}
+
+internal class DSTest {
+    private val x by symbol
+    private val y by symbol
+
+    @Test
+    fun dsAlgebraTest() {
+        diff(2, x to 1.0, y to 1.0) {
+            val x = bindSymbol(x)//by binding()
+            val y = bindSymbol("y")
+            val z = x * (-sin(x * y) + y) + 2.0
+            println(z.derivative(x))
+            println(z.derivative(y, x))
+            assertEquals(z.derivative(x, y), z.derivative(y, x))
+            // check improper order cause failure
+            assertFails { z.derivative(x, x, y) }
+        }
+    }
+
+    @Test
+    fun dsExpressionTest() {
+        val f = DSFieldExpression(DoubleField, ::DoubleBuffer) {
+            val x by binding
+            val y by binding
+            x.pow(2) + 2 * x * y + y.pow(2) + 1
+        }
+
+        assertEquals(10.0, f(x to 1.0, y to 2.0))
+        assertEquals(6.0, f.derivative(x)(x to 1.0, y to 2.0))
+        assertEquals(2.0, f.derivative(x, x)(x to 1.234, y to -2.0))
+        assertEquals(2.0, f.derivative(x, y)(x to 1.0, y to 2.0))
+    }
+}

kmath-histograms/src/commonMain/kotlin/space/kscience/kmath/histogram/UniformHistogramGroupND.kt

@@ -21,14 +21,14 @@ public typealias HyperSquareBin<V> = DomainBin<Double, HyperSquareDomain, V>
 
 /**
  * Multivariate histogram space for hyper-square real-field bins.
- * @param bufferFactory is an optional parameter used to optimize buffer production.
+ * @param valueBufferFactory is an optional parameter used to optimize buffer production.
  */
 public class UniformHistogramGroupND<V : Any, A : Field<V>>(
     override val valueAlgebraND: FieldOpsND<V, A>,
     private val lower: Buffer<Double>,
     private val upper: Buffer<Double>,
     private val binNums: IntArray = IntArray(lower.size) { 20 },
-    private val bufferFactory: BufferFactory<V> = Buffer.Companion::boxing,
+    private val valueBufferFactory: BufferFactory<V> = valueAlgebraND.elementAlgebra.bufferFactory,
 ) : HistogramGroupND<Double, HyperSquareDomain, V> {
 
     init {
@@ -94,7 +94,7 @@ public class UniformHistogramGroupND<V : Any, A : Field<V>>(
             }
         }
         hBuilder.apply(builder)
-        val values: BufferND<V> = ndCounter.mapToBuffer(bufferFactory) { it.value }
+        val values: BufferND<V> = ndCounter.mapToBuffer(valueBufferFactory) { it.value }
         return HistogramND(this, values)
     }
 
@@ -114,12 +114,12 @@ public class UniformHistogramGroupND<V : Any, A : Field<V>>(
 public fun <V : Any, A : Field<V>> Histogram.Companion.uniformNDFromRanges(
     valueAlgebraND: FieldOpsND<V, A>,
     vararg ranges: ClosedFloatingPointRange<Double>,
-    bufferFactory: BufferFactory<V> = Buffer.Companion::boxing,
+    bufferFactory: BufferFactory<V> = valueAlgebraND.elementAlgebra.bufferFactory,
 ): UniformHistogramGroupND<V, A> = UniformHistogramGroupND(
     valueAlgebraND,
     ranges.map(ClosedFloatingPointRange<Double>::start).asBuffer(),
     ranges.map(ClosedFloatingPointRange<Double>::endInclusive).asBuffer(),
-    bufferFactory = bufferFactory
+    valueBufferFactory = bufferFactory
 )
 
 public fun Histogram.Companion.uniformDoubleNDFromRanges(
@@ -140,7 +140,7 @@ public fun Histogram.Companion.uniformDoubleNDFromRanges(
 public fun <V : Any, A : Field<V>> Histogram.Companion.uniformNDFromRanges(
     valueAlgebraND: FieldOpsND<V, A>,
     vararg ranges: Pair<ClosedFloatingPointRange<Double>, Int>,
-    bufferFactory: BufferFactory<V> = Buffer.Companion::boxing,
+    bufferFactory: BufferFactory<V> = valueAlgebraND.elementAlgebra.bufferFactory,
 ): UniformHistogramGroupND<V, A> = UniformHistogramGroupND(
     valueAlgebraND,
     ListBuffer(
@@ -154,7 +154,7 @@ public fun <V : Any, A : Field<V>> Histogram.Companion.uniformNDFromRanges(
             .map(ClosedFloatingPointRange<Double>::endInclusive)
     ),
     ranges.map(Pair<ClosedFloatingPointRange<Double>, Int>::second).toIntArray(),
-    bufferFactory = bufferFactory
+    valueBufferFactory = bufferFactory
 )
 
 public fun Histogram.Companion.uniformDoubleNDFromRanges(

kmath-multik/src/main/kotlin/space/kscience/kmath/multik/MultikDoubleAlgebra.kt

@@ -6,6 +6,7 @@
 package space.kscience.kmath.multik
 
 import org.jetbrains.kotlinx.multik.ndarray.data.DataType
+import space.kscience.kmath.misc.PerformancePitfall
 import space.kscience.kmath.nd.StructureND
 import space.kscience.kmath.operations.DoubleField
 import space.kscience.kmath.operations.ExponentialOperations
@@ -22,10 +23,13 @@ public object MultikDoubleAlgebra : MultikDivisionTensorAlgebra<Double, DoubleFi
 
     override fun tan(arg: StructureND<Double>): MultikTensor<Double> = sin(arg) / cos(arg)
 
+    @PerformancePitfall
     override fun asin(arg: StructureND<Double>): MultikTensor<Double>  = arg.map { asin(it) }
 
+    @PerformancePitfall
     override fun acos(arg: StructureND<Double>): MultikTensor<Double> = arg.map { acos(it) }
 
+    @PerformancePitfall
     override fun atan(arg: StructureND<Double>): MultikTensor<Double> = arg.map { atan(it) }
 
     override fun exp(arg: StructureND<Double>): MultikTensor<Double> = multikMath.mathEx.exp(arg.asMultik().array).wrap()
@@ -42,10 +46,13 @@ public object MultikDoubleAlgebra : MultikDivisionTensorAlgebra<Double, DoubleFi
         return (expPlus - expMinus) / (expPlus + expMinus)
     }
 
+    @PerformancePitfall
     override fun asinh(arg: StructureND<Double>): MultikTensor<Double> = arg.map { asinh(it) }
 
+    @PerformancePitfall
     override fun acosh(arg: StructureND<Double>): MultikTensor<Double> = arg.map { acosh(it) }
 
+    @PerformancePitfall
     override fun atanh(arg: StructureND<Double>): MultikTensor<Double> = arg.map { atanh(it) }
 }
 

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

@@ -48,7 +48,8 @@ internal object InternalUtils {
                 cache.copyInto(
                     logFactorials,
                     BEGIN_LOG_FACTORIALS,
-                    BEGIN_LOG_FACTORIALS, endCopy
+                    BEGIN_LOG_FACTORIALS,
+                    endCopy,
                 )
             } else
             // All values to be computed

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

@@ -35,7 +35,7 @@ public fun interface Sampler<out T : Any> {
 public fun <T : Any> Sampler<T>.sampleBuffer(
     generator: RandomGenerator,
     size: Int,
-    bufferFactory: BufferFactory<T> = Buffer.Companion::boxing,
+    bufferFactory: BufferFactory<T> = BufferFactory.boxing(),
 ): Chain<Buffer<T>> {
     require(size > 1)
     //creating temporary storage once