diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSAlgebra.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSAlgebra.kt new file mode 100644 index 000000000..d9fc46b47 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSAlgebra.kt @@ -0,0 +1,437 @@ +/* + * 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> { + public val derivativeAlgebra: DSAlgebra + public val data: Buffer +} + +/** + * 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 +public fun > DS.getPartialDerivative(vararg orders: Int): T = + data[derivativeAlgebra.compiler.getPartialDerivativeIndex(*orders)] + +/** + * The value part of the derivative structure. + * + * @see getPartialDerivative + */ +@UnstableKMathAPI +public val > DS.value: T get() = data[0] + +@UnstableKMathAPI +public abstract class DSAlgebra>( + public val algebra: A, + public val bufferFactory: MutableBufferFactory, + public val order: Int, + bindings: Map, +) : ExpressionAlgebra> { + + @OptIn(UnstableKMathAPI::class) + private fun bufferForVariable(index: Int, value: T): Buffer { + val buffer = bufferFactory(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(numberOfVariables).apply { + set(index, 1) + }) + + buffer[indexOfDerivative] = algebra.one + } + return buffer + } + + @UnstableKMathAPI + protected inner class DSImpl internal constructor( + override val data: Buffer, + ) : DS { + override val derivativeAlgebra: DSAlgebra get() = this@DSAlgebra + } + + protected fun DS(data: Buffer): DS = 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 { + 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 { + 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 { + override val derivativeAlgebra: DSAlgebra get() = this@DSAlgebra + override val data: Buffer = bufferForVariable(index, value) + } + + + public val numberOfVariables: Int = bindings.size + + /** + * Get the compiler for number of free parameters and order. + * + * @return cached rules set. + */ + @PublishedApi + internal val compiler: DSCompiler by lazy { + // get the cached compilers + val cache: Array?>>? = 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(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? = if (p == 0) null else newCache[p - 1][o]!! + val derivativeCompiler: DSCompiler? = if (o == 0) null else newCache[p][o - 1]!! + + newCache[p][o] = DSCompiler( + algebra, + bufferFactory, + p, + o, + valueCompiler, + derivativeCompiler, + ) + } + } + } + + return@lazy newCache[numberOfVariables][order]!! + } + + private val variables: Map = bindings.entries.mapIndexed { index, (key, value) -> + key to DSSymbol( + index, + key, + value, + ) + }.toMap() + + + public override fun const(value: T): DS { + val buffer = bufferFactory(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.derivative(symbols: List): 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.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( + algebra: A, + bufferFactory: MutableBufferFactory, + order: Int, + bindings: Map, +) : DSAlgebra(algebra, bufferFactory, order, bindings), + Ring>, ScaleOperations>, + NumericAlgebra>, + NumbersAddOps> where A : Ring, A : NumericAlgebra, A : ScaleOperations { + + override fun bindSymbolOrNull(value: String): DSSymbol? = + super.bindSymbolOrNull(value) + + override fun DS.unaryMinus(): DS = mapData { -it } + + /** + * Create a copy of given [Buffer] and modify it according to [block] + */ + protected inline fun DS.transformDataBuffer(block: A.(MutableBuffer) -> Unit): DS { + require(derivativeAlgebra == this@DSRing) { "All derivative operations should be done in the same algebra" } + val newData = bufferFactory(compiler.size) { data[it] } + algebra.block(newData) + return DS(newData) + } + + protected fun DS.mapData(block: A.(T) -> T): DS { + require(derivativeAlgebra == this@DSRing) { "All derivative operations should be done in the same algebra" } + val newData: Buffer = data.map(bufferFactory) { + algebra.block(it) + } + return DS(newData) + } + + protected fun DS.mapDataIndexed(block: (Int, T) -> T): DS { + require(derivativeAlgebra == this@DSRing) { "All derivative operations should be done in the same algebra" } + val newData: Buffer = data.mapIndexed(bufferFactory, block) + return DS(newData) + } + + override val zero: DS by lazy { + const(algebra.zero) + } + + override val one: DS by lazy { + const(algebra.one) + } + + override fun number(value: Number): DS = const(algebra.number(value)) + + override fun add(left: DS, right: DS): DS = 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, value: Double): DS = a.mapData { + it.times(value) + } + + override fun multiply( + left: DS, + right: DS, + ): DS = left.transformDataBuffer { result -> + compiler.multiply(left.data, 0, right.data, 0, result, 0) + } +// +// override fun DS.minus(arg: DS): DS = transformDataBuffer { result -> +// subtract(data, 0, arg.data, 0, result, 0) +// } + + override operator fun DS.plus(other: Number): DS = transformDataBuffer { + it[0] += number(other) + } + +// +// override operator fun DS.minus(other: Number): DS = +// this + (-other.toDouble()) + + override operator fun Number.plus(other: DS): DS = other + this + override operator fun Number.minus(other: DS): DS = other - this +} + +@UnstableKMathAPI +public class DerivativeStructureRingExpression( + public val algebra: A, + public val bufferFactory: MutableBufferFactory, + public val function: DSRing.() -> DS, +) : DifferentiableExpression where A : Ring, A : ScaleOperations, A : NumericAlgebra { + override operator fun invoke(arguments: Map): T = + DSRing(algebra, bufferFactory, 0, arguments).function().value + + override fun derivativeOrNull(symbols: List): Expression = Expression { arguments -> + with( + DSRing( + algebra, + bufferFactory, + symbols.size, + arguments + ) + ) { 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>( + algebra: A, + bufferFactory: MutableBufferFactory, + order: Int, + bindings: Map, +) : DSRing(algebra, bufferFactory, order, bindings), ExtendedField> { + override fun number(value: Number): DS = const(algebra.number(value)) + + override fun divide(left: DS, right: DS): DS = left.transformDataBuffer { result -> + compiler.divide(left.data, 0, right.data, 0, result, 0) + } + + override fun sin(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.sin(arg.data, 0, result, 0) + } + + override fun cos(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.cos(arg.data, 0, result, 0) + } + + override fun tan(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.tan(arg.data, 0, result, 0) + } + + override fun asin(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.asin(arg.data, 0, result, 0) + } + + override fun acos(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.acos(arg.data, 0, result, 0) + } + + override fun atan(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.atan(arg.data, 0, result, 0) + } + + override fun sinh(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.sinh(arg.data, 0, result, 0) + } + + override fun cosh(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.cosh(arg.data, 0, result, 0) + } + + override fun tanh(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.tanh(arg.data, 0, result, 0) + } + + override fun asinh(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.asinh(arg.data, 0, result, 0) + } + + override fun acosh(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.acosh(arg.data, 0, result, 0) + } + + override fun atanh(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.atanh(arg.data, 0, result, 0) + } + + override fun power(arg: DS, pow: Number): DS = 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): DS = arg.transformDataBuffer { result -> + compiler.sqrt(arg.data, 0, result, 0) + } + + public fun power(arg: DS, pow: DS): DS = arg.transformDataBuffer { result -> + compiler.pow(arg.data, 0, pow.data, 0, result, 0) + } + + override fun exp(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.exp(arg.data, 0, result, 0) + } + + override fun ln(arg: DS): DS = arg.transformDataBuffer { result -> + compiler.ln(arg.data, 0, result, 0) + } +} + +@UnstableKMathAPI +public class DerivativeStructureFieldExpression>( + public val algebra: A, + public val bufferFactory: MutableBufferFactory, + public val function: DSField.() -> DS, +) : DifferentiableExpression { + override operator fun invoke(arguments: Map): T = + DSField(algebra, bufferFactory, 0, arguments).function().value + + override fun derivativeOrNull(symbols: List): Expression = Expression { arguments -> + DSField( + algebra, + bufferFactory, + symbols.size, + arguments, + ).run { function().derivative(symbols) } + } +} diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DerivativeStructure.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DerivativeStructure.kt deleted file mode 100644 index 01c045cdb..000000000 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DerivativeStructure.kt +++ /dev/null @@ -1,157 +0,0 @@ -/* - * 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.Ring -import space.kscience.kmath.structures.Buffer -import space.kscience.kmath.structures.asBuffer - -/** - * 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 open class DerivativeStructure> @PublishedApi internal constructor( - private val derivativeAlgebra: DerivativeStructureAlgebra, - @PublishedApi internal val data: Buffer, -) { - - public val compiler: DSCompiler get() = derivativeAlgebra.compiler - - /** - * The number of free parameters. - */ - public val freeParameters: Int get() = compiler.freeParameters - - /** - * The derivation order. - */ - public val order: Int get() = compiler.order - - /** - * The value part of the derivative structure. - * - * @see getPartialDerivative - */ - public val value: T get() = data[0] - - /** - * 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 - */ - public fun getPartialDerivative(vararg orders: Int): T = data[compiler.getPartialDerivativeIndex(*orders)] - - - /** - * Test for the equality of two derivative structures. - * - * Derivative structures are considered equal if they have the same number - * of free parameters, the same derivation order, and the same derivatives. - * - * @return `true` if two derivative structures are equal. - */ - public override fun equals(other: Any?): Boolean { - if (this === other) return true - - if (other is DerivativeStructure<*, *>) { - return ((freeParameters == other.freeParameters) && - (order == other.order) && - data == other.data) - } - - return false - } - - public override fun hashCode(): Int = - 227 + 229 * freeParameters + 233 * order + 239 * data.hashCode() - - public companion object { - - /** - * 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( - derivativeAlgebra: DerivativeStructureAlgebra, - index: Int, - value: T, - ): DerivativeStructure { - val compiler = derivativeAlgebra.compiler - require(index < compiler.freeParameters) { "number is too large: $index >= ${compiler.freeParameters}" } - return DerivativeStructure(derivativeAlgebra, derivativeAlgebra.bufferForVariable(index, value)) - } - - /** - * Build an instance from all its derivatives. - * - * @param derivatives derivatives sorted according to [DSCompiler.getPartialDerivativeIndex]. - */ - public fun > ofDerivatives( - derivativeAlgebra: DerivativeStructureAlgebra, - vararg derivatives: T, - ): DerivativeStructure { - val compiler = derivativeAlgebra.compiler - require(derivatives.size == compiler.size) { "dimension mismatch: ${derivatives.size} and ${compiler.size}" } - val data = derivatives.asBuffer() - - return DerivativeStructure( - derivativeAlgebra, - data - ) - } - } -} - -@OptIn(UnstableKMathAPI::class) -private fun > DerivativeStructureAlgebra.bufferForVariable(index: Int, value: T): Buffer { - val buffer = bufferFactory(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(numberOfVariables).apply { - set(index, 1) - }) - - buffer[indexOfDerivative] = algebra.one - } - return buffer -} - -/** - * A class implementing both [DerivativeStructure] and [Symbol]. - */ -@UnstableKMathAPI -public class DerivativeStructureSymbol> internal constructor( - derivativeAlgebra: DerivativeStructureAlgebra, - index: Int, - symbol: Symbol, - value: T, -) : Symbol by symbol, DerivativeStructure( - derivativeAlgebra, derivativeAlgebra.bufferForVariable(index, value) -) { - override fun toString(): String = symbol.toString() - override fun equals(other: Any?): Boolean = (other as? Symbol) == symbol - override fun hashCode(): Int = symbol.hashCode() -} diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpression.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpression.kt deleted file mode 100644 index 638057921..000000000 --- a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpression.kt +++ /dev/null @@ -1,349 +0,0 @@ -/* - * 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 kotlin.math.max -import kotlin.math.min - -@UnstableKMathAPI -public abstract class DerivativeStructureAlgebra>( - public val algebra: A, - public val bufferFactory: MutableBufferFactory, - public val order: Int, - bindings: Map, -) : ExpressionAlgebra> { - - public val numberOfVariables: Int = bindings.size - - - /** - * Get the compiler for number of free parameters and order. - * - * @return cached rules set. - */ - @PublishedApi - internal val compiler: DSCompiler by lazy { - // get the cached compilers - val cache: Array?>>? = 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(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? = if (p == 0) null else newCache[p - 1][o]!! - val derivativeCompiler: DSCompiler? = if (o == 0) null else newCache[p][o - 1]!! - - newCache[p][o] = DSCompiler( - algebra, - bufferFactory, - p, - o, - valueCompiler, - derivativeCompiler, - ) - } - } - } - - return@lazy newCache[numberOfVariables][order]!! - } - - private val variables: Map> = - bindings.entries.mapIndexed { index, (key, value) -> - key to DerivativeStructureSymbol( - this, - index, - key, - value, - ) - }.toMap() - - - - public override fun const(value: T): DerivativeStructure { - val buffer = bufferFactory(compiler.size) { algebra.zero } - buffer[0] = value - - return DerivativeStructure( - this, - buffer - ) - } - - override fun bindSymbolOrNull(value: String): DerivativeStructureSymbol? = variables[StringSymbol(value)] - - override fun bindSymbol(value: String): DerivativeStructureSymbol = - bindSymbolOrNull(value) ?: error("Symbol '$value' is not supported in $this") - - public fun bindSymbolOrNull(symbol: Symbol): DerivativeStructureSymbol? = variables[symbol.identity] - - public fun bindSymbol(symbol: Symbol): DerivativeStructureSymbol = - bindSymbolOrNull(symbol.identity) ?: error("Symbol '${symbol}' is not supported in $this") - - public fun DerivativeStructure.derivative(symbols: List): 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 DerivativeStructure.derivative(vararg symbols: Symbol): T = derivative(symbols.toList()) - -} - - -/** - * A ring over [DerivativeStructure]. - * - * @property order The derivation order. - * @param bindings The map of bindings values. All bindings are considered free parameters. - */ -@UnstableKMathAPI -public open class DerivativeStructureRing( - algebra: A, - bufferFactory: MutableBufferFactory, - order: Int, - bindings: Map, -) : DerivativeStructureAlgebra(algebra, bufferFactory, order, bindings), - Ring>, ScaleOperations>, - NumericAlgebra>, - NumbersAddOps> where A : Ring, A : NumericAlgebra, A : ScaleOperations { - - override fun bindSymbolOrNull(value: String): DerivativeStructureSymbol? = - super.bindSymbolOrNull(value) - - override fun DerivativeStructure.unaryMinus(): DerivativeStructure { - val newData = algebra { data.map(bufferFactory) { -it } } - return DerivativeStructure(this@DerivativeStructureRing, newData) - } - - /** - * Create a copy of given [Buffer] and modify it according to [block] - */ - protected inline fun DerivativeStructure.transformDataBuffer(block: DSCompiler.(MutableBuffer) -> Unit): DerivativeStructure { - val newData = bufferFactory(compiler.size) { data[it] } - compiler.block(newData) - return DerivativeStructure(this@DerivativeStructureRing, newData) - } - - protected fun DerivativeStructure.mapData(block: (T) -> T): DerivativeStructure { - val newData: Buffer = data.map(bufferFactory, block) - return DerivativeStructure(this@DerivativeStructureRing, newData) - } - - protected fun DerivativeStructure.mapDataIndexed(block: (Int, T) -> T): DerivativeStructure { - val newData: Buffer = data.mapIndexed(bufferFactory, block) - return DerivativeStructure(this@DerivativeStructureRing, newData) - } - - override val zero: DerivativeStructure by lazy { - const(algebra.zero) - } - - override val one: DerivativeStructure by lazy { - const(algebra.one) - } - - override fun number(value: Number): DerivativeStructure = const(algebra.number(value)) - - override fun add(left: DerivativeStructure, right: DerivativeStructure): DerivativeStructure { - left.compiler.checkCompatibility(right.compiler) - return left.transformDataBuffer { result -> - add(left.data, 0, right.data, 0, result, 0) - } - } - - override fun scale(a: DerivativeStructure, value: Double): DerivativeStructure = algebra { - a.mapData { it.times(value) } - } - - override fun multiply( - left: DerivativeStructure, - right: DerivativeStructure, - ): DerivativeStructure { - left.compiler.checkCompatibility(right.compiler) - return left.transformDataBuffer { result -> - multiply(left.data, 0, right.data, 0, result, 0) - } - } - - override fun DerivativeStructure.minus(arg: DerivativeStructure): DerivativeStructure { - compiler.checkCompatibility(arg.compiler) - return transformDataBuffer { result -> - subtract(data, 0, arg.data, 0, result, 0) - } - } - - override operator fun DerivativeStructure.plus(other: Number): DerivativeStructure = algebra { - transformDataBuffer { - it[0] += number(other) - } - } - - override operator fun DerivativeStructure.minus(other: Number): DerivativeStructure = - this + (-other.toDouble()) - - override operator fun Number.plus(other: DerivativeStructure): DerivativeStructure = other + this - override operator fun Number.minus(other: DerivativeStructure): DerivativeStructure = other - this -} - -@UnstableKMathAPI -public class DerivativeStructureRingExpression( - public val algebra: A, - public val bufferFactory: MutableBufferFactory, - public val function: DerivativeStructureRing.() -> DerivativeStructure, -) : DifferentiableExpression where A : Ring, A : ScaleOperations, A : NumericAlgebra { - override operator fun invoke(arguments: Map): T = - DerivativeStructureRing(algebra, bufferFactory, 0, arguments).function().value - - override fun derivativeOrNull(symbols: List): Expression = Expression { arguments -> - with( - DerivativeStructureRing( - algebra, - bufferFactory, - symbols.size, - arguments - ) - ) { 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 DerivativeStructureField>( - algebra: A, - bufferFactory: MutableBufferFactory, - order: Int, - bindings: Map, -) : DerivativeStructureRing(algebra, bufferFactory, order, bindings), ExtendedField> { - override fun number(value: Number): DerivativeStructure = const(algebra.number(value)) - - override fun divide(left: DerivativeStructure, right: DerivativeStructure): DerivativeStructure { - left.compiler.checkCompatibility(right.compiler) - return left.transformDataBuffer { result -> - left.compiler.divide(left.data, 0, right.data, 0, result, 0) - } - } - - override fun sin(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - sin(arg.data, 0, result, 0) - } - - override fun cos(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - cos(arg.data, 0, result, 0) - } - - override fun tan(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - tan(arg.data, 0, result, 0) - } - - override fun asin(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - asin(arg.data, 0, result, 0) - } - - override fun acos(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - acos(arg.data, 0, result, 0) - } - - override fun atan(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - atan(arg.data, 0, result, 0) - } - - override fun sinh(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - sinh(arg.data, 0, result, 0) - } - - override fun cosh(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - cosh(arg.data, 0, result, 0) - } - - override fun tanh(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - tanh(arg.data, 0, result, 0) - } - - override fun asinh(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - asinh(arg.data, 0, result, 0) - } - - override fun acosh(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - acosh(arg.data, 0, result, 0) - } - - override fun atanh(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - atanh(arg.data, 0, result, 0) - } - - override fun power(arg: DerivativeStructure, pow: Number): DerivativeStructure = when (pow) { - is Int -> arg.transformDataBuffer { result -> - pow(arg.data, 0, pow, result, 0) - } - else -> arg.transformDataBuffer { result -> - pow(arg.data, 0, pow.toDouble(), result, 0) - } - } - - override fun sqrt(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - sqrt(arg.data, 0, result, 0) - } - - public fun power(arg: DerivativeStructure, pow: DerivativeStructure): DerivativeStructure { - arg.compiler.checkCompatibility(pow.compiler) - return arg.transformDataBuffer { result -> - pow(arg.data, 0, pow.data, 0, result, 0) - } - } - - override fun exp(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - exp(arg.data, 0, result, 0) - } - - override fun ln(arg: DerivativeStructure): DerivativeStructure = arg.transformDataBuffer { result -> - ln(arg.data, 0, result, 0) - } -} - -@UnstableKMathAPI -public class DerivativeStructureFieldExpression>( - public val algebra: A, - public val bufferFactory: MutableBufferFactory, - public val function: DerivativeStructureField.() -> DerivativeStructure, -) : DifferentiableExpression { - override operator fun invoke(arguments: Map): T = - DerivativeStructureField(algebra, bufferFactory, 0, arguments).function().value - - override fun derivativeOrNull(symbols: List): Expression = Expression { arguments -> - with( - DerivativeStructureField( - algebra, - bufferFactory, - symbols.size, - arguments, - ) - ) { function().derivative(symbols) } - } -} diff --git a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpressionTest.kt b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpressionTest.kt index 429fe310b..fdeda4512 100644 --- a/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpressionTest.kt +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpressionTest.kt @@ -19,10 +19,10 @@ import kotlin.test.assertFails internal inline fun diff( order: Int, vararg parameters: Pair, - block: DerivativeStructureField.() -> Unit, + block: DSField.() -> Unit, ) { contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) } - DerivativeStructureField(DoubleField, ::DoubleBuffer, order, mapOf(*parameters)).block() + DSField(DoubleField, ::DoubleBuffer, order, mapOf(*parameters)).block() } internal class AutoDiffTest {