From 7b1bdc21a4ba5eb5d36d32afeffde7055c26f88b Mon Sep 17 00:00:00 2001 From: Iaroslav Postovalov Date: Thu, 12 Aug 2021 16:37:53 +0300 Subject: [PATCH] Copy DerivativeStructure to multiplatform --- README.md | 2 +- docs/templates/README-TEMPLATE.md | 2 +- .../kscience/kmath/expressions/DSCompiler.kt | 1541 +++++++++++++++++ .../kmath/expressions/DerivativeStructure.kt | 186 ++ .../DerivativeStructureExpression.kt | 332 ++++ .../DerivativeStructureExpressionTest.kt | 59 + .../kscience/kmath/internal/InternalUtils.kt | 3 +- 7 files changed, 2122 insertions(+), 3 deletions(-) create mode 100644 kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSCompiler.kt create mode 100644 kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DerivativeStructure.kt create mode 100644 kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpression.kt create mode 100644 kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpressionTest.kt diff --git a/README.md b/README.md index 92260716e..aeedfefbb 100644 --- a/README.md +++ b/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. diff --git a/docs/templates/README-TEMPLATE.md b/docs/templates/README-TEMPLATE.md index b0c418697..2e64a3e09 100644 --- a/docs/templates/README-TEMPLATE.md +++ b/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. diff --git a/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSCompiler.kt b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSCompiler.kt new file mode 100644 index 000000000..bb88ce52c --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DSCompiler.kt @@ -0,0 +1,1541 @@ +/* + * 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.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 + +internal fun MutableBuffer.fill(element: T, fromIndex: Int = 0, toIndex: Int = size) { + for (i in fromIndex until toIndex) this[i] = element +} + +/** + * Class holding "compiled" computation rules for derivative structures. + * + * This class implements the computation rules 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. However, to avoid performances bottlenecks, the recursive rules are + * "compiled" once in an unfolded form. This class does this recursion unrolling and stores the computation rules as + * simple loops with pre-computed indirection arrays. + * + * This class maps all derivative computation into single dimension arrays that hold the value and partial derivatives. + * The class does not hold these arrays, which remains under the responsibility of the caller. For each combination of + * number of free parameters and derivation order, only one compiler is necessary, and this compiler will be used to + * perform computations on all arrays provided to it, which can represent hundreds or thousands of different parameters + * kept together with all their partial derivatives. + * + * The arrays on which compilers operate contain only the partial derivatives together with the 0th + * derivative, i.e., the value. The partial derivatives are stored in a compiler-specific order, which can be retrieved + * using methods [getPartialDerivativeIndex] and [getPartialDerivativeOrders]. The value is guaranteed to be stored as + * the first element (i.e., the [getPartialDerivativeIndex] method returns 0 when called with 0 for all derivation + * orders and [getPartialDerivativeOrders] returns an array filled with 0 when called with 0 as the index). + * + * Note that the ordering changes with number of parameters and derivation order. For example given 2 parameters x and + * y, df/dy is stored at index 2 when derivation order is set to 1 (in this case the array has three elements: f, + * df/dx and df/dy). If derivation order is set to 2, then df/dy will be stored at index 3 (in this case the array has + * six elements: f, df/dx, df/dxdx, df/dy, df/dxdy and df/dydy). + * + * Given this structure, users can perform some simple operations like adding, subtracting or multiplying constants and + * negating the elements by themselves, knowing if they want to mutate their array or create a new array. These simple + * operations are not provided by the compiler. The compiler provides only the more complex operations between several + * arrays. + * + * Derived from + * [Commons Math's `DSCompiler`](https://github.com/apache/commons-math/blob/924f6c357465b39beb50e3c916d5eb6662194175/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/differentiation/DSCompiler.java). + * + * @property freeParameters Number of free parameters. + * @property order Derivation order. + * @see DerivativeStructure + */ +internal class DSCompiler> internal constructor( + val algebra: A, + val bufferFactory: MutableBufferFactory, + val freeParameters: Int, + val order: Int, + valueCompiler: DSCompiler?, + derivativeCompiler: DSCompiler?, +) { + /** + * Number of partial derivatives (including the single 0 order derivative element). + */ + val sizes: Array by lazy { + compileSizes( + freeParameters, + order, + valueCompiler, + ) + } + + /** + * Indirection array for partial derivatives. + */ + val derivativesIndirection: Array by lazy { + compileDerivativesIndirection( + freeParameters, order, + valueCompiler, derivativeCompiler, + ) + } + + /** + * Indirection array of the lower derivative elements. + */ + val lowerIndirection: IntArray by lazy { + compileLowerIndirection( + freeParameters, order, + valueCompiler, derivativeCompiler, + ) + } + + /** + * Indirection arrays for multiplication. + */ + val multIndirection: Array> by lazy { + compileMultiplicationIndirection( + freeParameters, order, + valueCompiler, derivativeCompiler, lowerIndirection, + ) + } + + /** + * Indirection arrays for function composition. + */ + val compositionIndirection: Array> by lazy { + compileCompositionIndirection( + freeParameters, order, + valueCompiler, derivativeCompiler, + sizes, derivativesIndirection, + ) + } + + /** + * Get the array size required for holding partial derivatives' data. + * + * This number includes the single 0 order derivative element, which is + * guaranteed to be stored in the first element of the array. + */ + val size: Int + get() = sizes[freeParameters][order] + + /** + * Get the index of a partial derivative in the array. + * + * If all orders are set to 0, then the 0th order derivative is returned, which is the value of the + * function. + * + * The indices of derivatives are between 0 and [size] − 1. Their specific order is fixed for a given compiler, but + * otherwise not publicly specified. There are however some simple cases which have guaranteed indices: + * + * * the index of 0th order derivative is always 0 + * * if there is only 1 [freeParameters], then the + * derivatives are sorted in increasing derivation order (i.e., f at index 0, df/dp + * at index 1, d2f/dp2 at index 2 … + * dkf/dpk at index k), + * * if the [order] is 1, then the derivatives + * are sorted in increasing free parameter order (i.e., f at index 0, df/dx1 + * at index 1, df/dx2 at index 2 … df/dxk at index k), + * * all other cases are not publicly specified. + * + * This method is the inverse of method [getPartialDerivativeOrders]. + * + * @param orders derivation orders with respect to each parameter. + * @return index of the partial derivative. + * @see getPartialDerivativeOrders + */ + fun getPartialDerivativeIndex(vararg orders: Int): Int { + // safety check + require(orders.size == freeParameters) { "dimension mismatch: ${orders.size} and $freeParameters" } + return getPartialDerivativeIndex(freeParameters, order, sizes, *orders) + } + + /** + * Get the derivation orders for a specific index in the array. + * + * This method is the inverse of [getPartialDerivativeIndex]. + * + * @param index of the partial derivative + * @return orders derivation orders with respect to each parameter + * @see getPartialDerivativeIndex + */ + fun getPartialDerivativeOrders(index: Int): IntArray = derivativesIndirection[index] +} + +/** + * Compute natural logarithm of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for logarithm the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.ln( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Field, A : ExponentialOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + function[0] = ln(operand[operandOffset]) + + if (order > 0) { + val inv = one / operand[operandOffset] + var xk = inv + for (i in 1..order) { + function[i] = xk + xk *= (-i * inv) + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute integer power of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param n power to apply. + * @param result array where result must be stored (for power the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.pow( + operand: Buffer, + operandOffset: Int, + n: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Field, A : PowerOperations = algebra { + if (n == 0) { + // special case, x^0 = 1 for all x + result[resultOffset] = one + result.fill(zero, resultOffset + 1, resultOffset + size) + return + } + + // create the power function value and derivatives + // [x^n, nx^(n-1), n(n-1)x^(n-2), ... ] + val function = bufferFactory(1 + order) { zero } + + if (n > 0) { + // strictly positive power + val maxOrder: Int = min(order, n) + var xk = operand[operandOffset] pow n - maxOrder + for (i in maxOrder downTo 1) { + function[i] = xk + xk *= operand[operandOffset] + } + function[0] = xk + } else { + // strictly negative power + val inv = one / operand[operandOffset] + var xk = inv pow -n + + for (i in 0..order) { + function[i] = xk + xk *= inv + } + } + + var coefficient = number(n) + + for (i in 1..order) { + function[i] = function[i] * coefficient + coefficient *= (n - i).toDouble() + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute exponential of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for exponential the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.exp( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Ring, A : ScaleOperations, A : ExponentialOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + function.fill(exp(operand[operandOffset])) + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute square root of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for nth root the result array *cannot* be the input + * array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.sqrt( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Field, A : PowerOperations = algebra { + // create the function value and derivatives + // [x^(1/n), (1/n)x^((1/n)-1), (1-n)/n^2x^((1/n)-2), ... ] + val function = bufferFactory(1 + order) { zero } + function[0] = sqrt(operand[operandOffset]) + var xk: T = 0.5 * one / function[0] + val xReciprocal = one / operand[operandOffset] + + for (i in 1..order) { + function[i] = xk + xk *= xReciprocal * (0.5 - i) + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute cosine of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for cosine the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.cos( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int, +) where A : Ring, A : TrigonometricOperations, A : ScaleOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + function[0] = cos(operand[operandOffset]) + + if (order > 0) { + function[1] = -sin(operand[operandOffset]) + for (i in 2..order) { + function[i] = -function[i - 2] + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute power of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param p power to apply. + * @param result array where result must be stored (for power the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.pow( + operand: Buffer, + operandOffset: Int, + p: Double, + result: MutableBuffer, + resultOffset: Int +) where A : Ring, A : NumericAlgebra, A : PowerOperations, A : ScaleOperations = algebra { + // create the function value and derivatives + // [x^p, px^(p-1), p(p-1)x^(p-2), ... ] + val function = bufferFactory(1 + order) { zero } + var xk = operand[operandOffset] pow p - order + + for (i in order downTo 1) { + function[i] = xk + xk *= operand[operandOffset] + } + + function[0] = xk + var coefficient = p + + for (i in 1..order) { + function[i] = function[i] * coefficient + coefficient *= p - i + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute tangent of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for tangent the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.tan( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Ring, A : TrigonometricOperations, A : ScaleOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + val t = tan(operand[operandOffset]) + function[0] = t + + if (order > 0) { + + // the nth order derivative of tan has the form: + // dn(tan(x)/dxn = P_n(tan(x)) + // where P_n(t) is a degree n+1 polynomial with same parity as n+1 + // P_0(t) = t, P_1(t) = 1 + t^2, P_2(t) = 2 t (1 + t^2) ... + // the general recurrence relation for P_n is: + // P_n(x) = (1+t^2) P_(n-1)'(t) + // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array + val p = bufferFactory(order + 2) { zero } + p[1] = one + val t2 = t * t + for (n in 1..order) { + + // update and evaluate polynomial P_n(t) + var v = one + p[n + 1] = n * p[n] + var k = n + 1 + while (k >= 0) { + v = v * t2 + p[k] + if (k > 2) { + p[k - 2] = (k - 1) * p[k - 1] + (k - 3) * p[k - 3] + } else if (k == 2) { + p[0] = p[1] + } + k -= 2 + } + if (n and 0x1 == 0) { + v *= t + } + function[n] = v + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute power of a derivative structure. + * + * @param x array holding the base. + * @param xOffset offset of the base in its array. + * @param y array holding the exponent. + * @param yOffset offset of the exponent in its array. + * @param result array where result must be stored (for power the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.pow( + x: Buffer, + xOffset: Int, + y: Buffer, + yOffset: Int, + result: MutableBuffer, + resultOffset: Int, +) where A : Field, A : ExponentialOperations = algebra { + val logX = bufferFactory(size) { zero } + ln(x, xOffset, logX, 0) + val yLogX = bufferFactory(size) { zero } + multiply(logX, 0, y, yOffset, yLogX, 0) + exp(yLogX, 0, result, resultOffset) +} + +/** + * Compute sine of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for sine the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.sin( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Ring, A : ScaleOperations, A : TrigonometricOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + function[0] = sin(operand[operandOffset]) + if (order > 0) { + function[1] = cos(operand[operandOffset]) + for (i in 2..order) { + function[i] = -function[i - 2] + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute arc cosine of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for arc cosine the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.acos( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Field, A : TrigonometricOperations, A : PowerOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + val x = operand[operandOffset] + function[0] = acos(x) + if (order > 0) { + // the nth order derivative of acos has the form: + // dn(acos(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2) + // where P_n(x) is a degree n-1 polynomial with same parity as n-1 + // P_1(x) = -1, P_2(x) = -x, P_3(x) = -2x^2 - 1 ... + // the general recurrence relation for P_n is: + // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x) + // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array + val p = bufferFactory(order) { zero } + p[0] = -one + val x2 = x * x + val f = one / (one - x2) + var coeff = sqrt(f) + function[1] = coeff * p[0] + + for (n in 2..order) { + // update and evaluate polynomial P_n(x) + var v = zero + p[n - 1] = (n - 1) * p[n - 2] + var k = n - 1 + + while (k >= 0) { + v = v * x2 + p[k] + if (k > 2) { + p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3] + } else if (k == 2) { + p[0] = p[1] + } + k -= 2 + } + + if (n and 0x1 == 0) { + v *= x + } + + coeff *= f + function[n] = coeff * v + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute arc sine of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for arc sine the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.asin( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Field, A : TrigonometricOperations, A : PowerOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + val x = operand[operandOffset] + function[0] = asin(x) + if (order > 0) { + // the nth order derivative of asin has the form: + // dn(asin(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2) + // where P_n(x) is a degree n-1 polynomial with same parity as n-1 + // P_1(x) = 1, P_2(x) = x, P_3(x) = 2x^2 + 1 ... + // the general recurrence relation for P_n is: + // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x) + // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array + val p = bufferFactory(order) { zero } + p[0] = one + val x2 = x * x + val f = one / (one - x2) + var coeff = sqrt(f) + function[1] = coeff * p[0] + for (n in 2..order) { + + // update and evaluate polynomial P_n(x) + var v = zero + p[n - 1] = (n - 1) * p[n - 2] + var k = n - 1 + while (k >= 0) { + v = v * x2 + p[k] + if (k > 2) { + p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3] + } else if (k == 2) { + p[0] = p[1] + } + k -= 2 + } + if (n and 0x1 == 0) { + v *= x + } + coeff *= f + function[n] = coeff * v + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute arc tangent of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for arc tangent the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.atan( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Field, A : TrigonometricOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + val x = operand[operandOffset] + function[0] = atan(x) + + if (order > 0) { + // the nth order derivative of atan has the form: + // dn(atan(x)/dxn = Q_n(x) / (1 + x^2)^n + // where Q_n(x) is a degree n-1 polynomial with same parity as n-1 + // Q_1(x) = 1, Q_2(x) = -2x, Q_3(x) = 6x^2 - 2 ... + // the general recurrence relation for Q_n is: + // Q_n(x) = (1+x^2) Q_(n-1)'(x) - 2(n-1) x Q_(n-1)(x) + // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array + val q = bufferFactory(order) { zero } + q[0] = one + val x2 = x * x + val f = one / (one + x2) + var coeff = f + function[1] = coeff * q[0] + for (n in 2..order) { + + // update and evaluate polynomial Q_n(x) + var v = zero + q[n - 1] = -n * q[n - 2] + var k = n - 1 + while (k >= 0) { + v = v * x2 + q[k] + if (k > 2) { + q[k - 2] = (k - 1) * q[k - 1] + (k - 1 - 2 * n) * q[k - 3] + } else if (k == 2) { + q[0] = q[1] + } + k -= 2 + } + if (n and 0x1 == 0) { + v *= x + } + coeff *= f + function[n] = coeff * v + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute hyperbolic cosine of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for hyperbolic cosine the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.cosh( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Ring, A : ScaleOperations, A : ExponentialOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + function[0] = cosh(operand[operandOffset]) + + if (order > 0) { + function[1] = sinh(operand[operandOffset]) + for (i in 2..order) { + function[i] = function[i - 2] + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute hyperbolic tangent of a derivative structure. + * + * @param operand array holding the operand + * @param operandOffset offset of the operand in its array + * @param result array where result must be stored (for hyperbolic tangent the result array *cannot* be the input + * array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.tanh( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Field, A : ExponentialOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + val t = tanh(operand[operandOffset]) + function[0] = t + if (order > 0) { + + // the nth order derivative of tanh has the form: + // dn(tanh(x)/dxn = P_n(tanh(x)) + // where P_n(t) is a degree n+1 polynomial with same parity as n+1 + // P_0(t) = t, P_1(t) = 1 - t^2, P_2(t) = -2 t (1 - t^2) ... + // the general recurrence relation for P_n is: + // P_n(x) = (1-t^2) P_(n-1)'(t) + // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array + val p = bufferFactory(order + 2) { zero } + p[1] = one + val t2 = t * t + for (n in 1..order) { + + // update and evaluate polynomial P_n(t) + var v = zero + p[n + 1] = -n * p[n] + var k = n + 1 + while (k >= 0) { + v = v * t2 + p[k] + if (k > 2) { + p[k - 2] = (k - 1) * p[k - 1] - (k - 3) * p[k - 3] + } else if (k == 2) { + p[0] = p[1] + } + k -= 2 + } + if (n and 0x1 == 0) { + v *= t + } + function[n] = v + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute inverse hyperbolic cosine of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for inverse hyperbolic cosine the result array *cannot* be the input + * array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.acosh( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Field, A : ExponentialOperations, A : PowerOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + val x = operand[operandOffset] + function[0] = acosh(x) + + if (order > 0) { + // the nth order derivative of acosh has the form: + // dn(acosh(x)/dxn = P_n(x) / [x^2 - 1]^((2n-1)/2) + // where P_n(x) is a degree n-1 polynomial with same parity as n-1 + // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 + 1 ... + // the general recurrence relation for P_n is: + // P_n(x) = (x^2-1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x) + // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array + val p = bufferFactory(order) { zero } + p[0] = one + val x2 = x * x + val f = one / (x2 - one) + var coeff = sqrt(f) + function[1] = coeff * p[0] + for (n in 2..order) { + + // update and evaluate polynomial P_n(x) + var v = zero + p[n - 1] = (1 - n) * p[n - 2] + var k = n - 1 + while (k >= 0) { + v = v * x2 + p[k] + if (k > 2) { + p[k - 2] = (1 - k) * p[k - 1] + (k - 2 * n) * p[k - 3] + } else if (k == 2) { + p[0] = -p[1] + } + k -= 2 + } + if (n and 0x1 == 0) { + v *= x + } + + coeff *= f + function[n] = coeff * v + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Compute composition of a derivative structure by a function. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param f array of value and derivatives of the function at the current point (i.e. at `operand[operandOffset]`). + * @param result array where result must be stored (for composition the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.compose( + operand: Buffer, + operandOffset: Int, + f: Buffer, + result: MutableBuffer, + resultOffset: Int, +) where A : Ring, A : ScaleOperations = algebra { + for (i in compositionIndirection.indices) { + val mappingI = compositionIndirection[i] + var r = zero + for (j in mappingI.indices) { + val mappingIJ = mappingI[j] + var product = mappingIJ[0] * f[mappingIJ[1]] + for (k in 2 until mappingIJ.size) { + product *= operand[operandOffset + mappingIJ[k]] + } + r += product + } + result[resultOffset + i] = r + } +} + +/** + * Compute hyperbolic sine of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for hyperbolic sine the result array *cannot* be the input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.sinh( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Field, A : ExponentialOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + function[0] = sinh(operand[operandOffset]) + + if (order > 0) { + function[1] = cosh(operand[operandOffset]) + for (i in 2..order) { + function[i] = function[i - 2] + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Perform division of two derivative structures. + * + * @param lhs array holding left-hand side of division. + * @param lhsOffset offset of the left-hand side in its array. + * @param rhs array right-hand side of division. + * @param rhsOffset offset of the right-hand side in its array. + * @param result array where result must be stored (for division the result array *cannot* be one of the input arrays). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.divide( + lhs: Buffer, + lhsOffset: Int, + rhs: Buffer, + rhsOffset: Int, + result: MutableBuffer, + resultOffset: Int, +) where A : Field, A : PowerOperations, A : ScaleOperations = algebra { + val reciprocal = bufferFactory(size) { zero } + pow(rhs, lhsOffset, -1, reciprocal, 0) + multiply(lhs, lhsOffset, reciprocal, rhsOffset, result, resultOffset) +} + +/** + * Perform multiplication of two derivative structures. + * + * @param lhs array holding left-hand side of multiplication. + * @param lhsOffset offset of the left-hand side in its array. + * @param rhs array right-hand side of multiplication. + * @param rhsOffset offset of the right-hand side in its array. + * @param result array where result must be stored (for multiplication the result array *cannot* be one of the input + * arrays). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.multiply( + lhs: Buffer, + lhsOffset: Int, + rhs: Buffer, + rhsOffset: Int, + result: MutableBuffer, + resultOffset: Int, +) where A : Ring, A : ScaleOperations = algebra { + for (i in multIndirection.indices) { + val mappingI = multIndirection[i] + var r = zero + + for (j in mappingI.indices) { + r += mappingI[j][0] * lhs[lhsOffset + mappingI[j][1]] * rhs[rhsOffset + mappingI[j][2]] + } + + result[resultOffset + i] = r + } +} + +/** + * Perform subtraction of two derivative structures. + * + * @param lhs array holding left-hand side of subtraction. + * @param lhsOffset offset of the left-hand side in its array. + * @param rhs array right-hand side of subtraction. + * @param rhsOffset offset of the right-hand side in its array. + * @param result array where result must be stored (it may be one of the input arrays). + * @param resultOffset offset of the result in its array. + */ +internal fun > DSCompiler.subtract( + lhs: Buffer, + lhsOffset: Int, + rhs: Buffer, + rhsOffset: Int, + result: MutableBuffer, + resultOffset: Int, +) = algebra { + for (i in 0 until size) { + result[resultOffset + i] = lhs[lhsOffset + i] - rhs[rhsOffset + i] + } +} + +/** + * Compute inverse hyperbolic sine of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for inverse hyperbolic sine the result array *cannot* be the input + * array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.asinh( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int +) where A : Field, A : ExponentialOperations, A : PowerOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + val x = operand[operandOffset] + function[0] = asinh(x) + if (order > 0) { + // the nth order derivative of asinh has the form: + // dn(asinh(x)/dxn = P_n(x) / [x^2 + 1]^((2n-1)/2) + // where P_n(x) is a degree n-1 polynomial with same parity as n-1 + // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 - 1 ... + // the general recurrence relation for P_n is: + // P_n(x) = (x^2+1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x) + // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array + val p = bufferFactory(order) { zero } + p[0] = one + val x2 = x * x + val f = one / (one + x2) + var coeff = sqrt(f) + function[1] = coeff * p[0] + for (n in 2..order) { + + // update and evaluate polynomial P_n(x) + var v = zero + p[n - 1] = (1 - n) * p[n - 2] + var k = n - 1 + while (k >= 0) { + v = v * x2 + p[k] + if (k > 2) { + p[k - 2] = (k - 1) * p[k - 1] + (k - 2 * n) * p[k - 3] + } else if (k == 2) { + p[0] = p[1] + } + k -= 2 + } + if (n and 0x1 == 0) { + v *= x + } + coeff *= f + function[n] = coeff * v + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Perform addition of two derivative structures. + * + * @param lhs array holding left-hand side of addition. + * @param lhsOffset offset of the left-hand side in its array. + * @param rhs array right-hand side of addition. + * @param rhsOffset offset of the right-hand side in its array. + * @param result array where result must be stored (it may be one of the input arrays). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.add( + lhs: Buffer, + lhsOffset: Int, + rhs: Buffer, + rhsOffset: Int, + result: MutableBuffer, + resultOffset: Int, +) where A : Group = algebra { + for (i in 0 until size) { + result[resultOffset + i] = lhs[lhsOffset + i] + rhs[rhsOffset + i] + } +} + +/** + * Check rules set compatibility. + * + * @param compiler other compiler to check against instance. + */ +internal fun > DSCompiler.checkCompatibility(compiler: DSCompiler) { + require(freeParameters == compiler.freeParameters) { + "dimension mismatch: $freeParameters and ${compiler.freeParameters}" + } + require(order == compiler.order) { + "dimension mismatch: $order and ${compiler.order}" + } +} + +/** + * Compute inverse hyperbolic tangent of a derivative structure. + * + * @param operand array holding the operand. + * @param operandOffset offset of the operand in its array. + * @param result array where result must be stored (for inverse hyperbolic tangent the result array *cannot* be the + * input array). + * @param resultOffset offset of the result in its array. + */ +internal fun DSCompiler.atanh( + operand: Buffer, + operandOffset: Int, + result: MutableBuffer, + resultOffset: Int, +) where A : Field, A : ExponentialOperations = algebra { + // create the function value and derivatives + val function = bufferFactory(1 + order) { zero } + val x = operand[operandOffset] + function[0] = atanh(x) + + if (order > 0) { + // the nth order derivative of atanh has the form: + // dn(atanh(x)/dxn = Q_n(x) / (1 - x^2)^n + // where Q_n(x) is a degree n-1 polynomial with same parity as n-1 + // Q_1(x) = 1, Q_2(x) = 2x, Q_3(x) = 6x^2 + 2 ... + // the general recurrence relation for Q_n is: + // Q_n(x) = (1-x^2) Q_(n-1)'(x) + 2(n-1) x Q_(n-1)(x) + // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array + val q = bufferFactory(order) { zero } + q[0] = one + val x2 = x * x + val f = one / (one - x2) + var coeff = f + function[1] = coeff * q[0] + for (n in 2..order) { + + // update and evaluate polynomial Q_n(x) + var v = zero + q[n - 1] = n * q[n - 2] + var k = n - 1 + while (k >= 0) { + v = v * x2 + q[k] + if (k > 2) { + q[k - 2] = (k - 1) * q[k - 1] + (2 * n - k + 1) * q[k - 3] + } else if (k == 2) { + q[0] = q[1] + } + k -= 2 + } + if (n and 0x1 == 0) { + v *= x + } + coeff *= f + function[n] = coeff * v + } + } + + // apply function composition + compose(operand, operandOffset, function, result, resultOffset) +} + +/** + * Get the compiler for number of free parameters and order. + * + * @param parameters number of free parameters. + * @param order derivation order. + * @return cached rules set. + */ +internal fun > getCompiler( + algebra: A, + bufferFactory: MutableBufferFactory, + parameters: Int, + order: Int +): DSCompiler { + // get the cached compilers + val cache: Array?>>? = null + + // we need to create more compilers + val maxParameters: Int = max(parameters, 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 + + // create the array in increasing diagonal order + for (diag in 0..parameters + order) { + for (o in max(0, diag - parameters)..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 newCache[parameters][order]!! +} + +/** + * Compile the sizes array. + * + * @param parameters number of free parameters. + * @param order derivation order. + * @param valueCompiler compiler for the value part. + * @return sizes array. + */ +private fun > compileSizes( + parameters: Int, order: Int, + valueCompiler: DSCompiler?, +): Array { + val sizes = Array(parameters + 1) { + IntArray(order + 1) + } + + if (parameters == 0) { + sizes[0].fill(1) + } else { + checkNotNull(valueCompiler) + valueCompiler.sizes.copyInto(sizes, endIndex = parameters) + sizes[parameters][0] = 1 + for (i in 0 until order) { + sizes[parameters][i + 1] = sizes[parameters][i] + sizes[parameters - 1][i + 1] + } + } + return sizes +} + +/** + * Compile the derivatives' indirection array. + * + * @param parameters number of free parameters. + * @param order derivation order. + * @param valueCompiler compiler for the value part. + * @param derivativeCompiler compiler for the derivative part. + * @return derivatives indirection array. + */ +private fun > compileDerivativesIndirection( + parameters: Int, + order: Int, + valueCompiler: DSCompiler?, + derivativeCompiler: DSCompiler?, +): Array { + if (parameters == 0 || order == 0) { + return Array(1) { IntArray(parameters) } + } + + val vSize: Int = valueCompiler!!.derivativesIndirection.size + val dSize: Int = derivativeCompiler!!.derivativesIndirection.size + val derivativesIndirection = Array(vSize + dSize) { IntArray(parameters) } + + // set up the indices for the value part + for (i in 0 until vSize) { + // copy the first indices, the last one remaining set to 0 + valueCompiler.derivativesIndirection[i].copyInto(derivativesIndirection[i], endIndex = parameters - 1) + } + + // set up the indices for the derivative part + for (i in 0 until dSize) { + // copy the indices + derivativeCompiler.derivativesIndirection[i].copyInto(derivativesIndirection[vSize], 0, 0, parameters) + + // increment the derivation order for the last parameter + derivativesIndirection[vSize + i][parameters - 1]++ + } + + return derivativesIndirection +} + +/** + * Compile the lower derivatives' indirection array. + * + * This indirection array contains the indices of all elements except derivatives for last derivation order. + * + * @param parameters number of free parameters. + * @param order derivation order. + * @param valueCompiler compiler for the value part. + * @param derivativeCompiler compiler for the derivative part. + * @return lower derivatives' indirection array. + */ +private fun > compileLowerIndirection( + parameters: Int, + order: Int, + valueCompiler: DSCompiler?, + derivativeCompiler: DSCompiler?, +): IntArray { + if (parameters == 0 || order <= 1) return intArrayOf(0) + checkNotNull(valueCompiler) + checkNotNull(derivativeCompiler) + + // this is an implementation of definition 6 in Dan Kalman's paper. + val vSize: Int = valueCompiler.lowerIndirection.size + val dSize: Int = derivativeCompiler.lowerIndirection.size + val lowerIndirection = IntArray(vSize + dSize) + valueCompiler.lowerIndirection.copyInto(lowerIndirection, endIndex = vSize) + for (i in 0 until dSize) { + lowerIndirection[vSize + i] = valueCompiler.size + derivativeCompiler.lowerIndirection[i] + } + return lowerIndirection +} + +/** + * Compile the multiplication indirection array. + * + * This indirection array contains the indices of all pairs of elements involved when computing a multiplication. This + * allows a straightforward loop-based multiplication (see [multiply]). + * + * @param parameters number of free parameters. + * @param order derivation order. + * @param valueCompiler compiler for the value part. + * @param derivativeCompiler compiler for the derivative part. + * @param lowerIndirection lower derivatives' indirection array. + * @return multiplication indirection array. + */ +@Suppress("UNCHECKED_CAST") +private fun > compileMultiplicationIndirection( + parameters: Int, + order: Int, + valueCompiler: DSCompiler?, + derivativeCompiler: DSCompiler?, + lowerIndirection: IntArray, +): Array> { + if (parameters == 0 || order == 0) return arrayOf(arrayOf(intArrayOf(1, 0, 0))) + + // this is an implementation of definition 3 in Dan Kalman's paper. + val vSize = valueCompiler!!.multIndirection.size + val dSize = derivativeCompiler!!.multIndirection.size + val multIndirection: Array?> = arrayOfNulls(vSize + dSize) + valueCompiler.multIndirection.copyInto(multIndirection, endIndex = vSize) + + for (i in 0 until dSize) { + val dRow = derivativeCompiler.multIndirection[i] + val row: List = buildList(dRow.size * 2) { + for (j in dRow.indices) { + add(intArrayOf(dRow[j][0], lowerIndirection[dRow[j][1]], vSize + dRow[j][2])) + add(intArrayOf(dRow[j][0], vSize + dRow[j][1], lowerIndirection[dRow[j][2]])) + } + } + + // combine terms with similar derivation orders + val combined: List = buildList(row.size) { + for (j in row.indices) { + val termJ = row[j] + if (termJ[0] > 0) { + for (k in j + 1 until row.size) { + val termK = row[k] + + if (termJ[1] == termK[1] && termJ[2] == termK[2]) { + // combine termJ and termK + termJ[0] += termK[0] + // make sure we will skip termK later on in the outer loop + termK[0] = 0 + } + } + + add(termJ) + } + } + } + + multIndirection[vSize + i] = combined.toTypedArray() + } + + return multIndirection as Array> +} + +/** + * Compile the indirection array of function composition. + * + * This indirection array contains the indices of all sets of elements involved when computing a composition. This + * allows a straightforward loop-based composition (see [compose]). + * + * @param parameters number of free parameters. + * @param order derivation order. + * @param valueCompiler compiler for the value part. + * @param derivativeCompiler compiler for the derivative part. + * @param sizes sizes array. + * @param derivativesIndirection derivatives indirection array. + * @return multiplication indirection array. + */ +@Suppress("UNCHECKED_CAST") +private fun > compileCompositionIndirection( + parameters: Int, + order: Int, + valueCompiler: DSCompiler?, + derivativeCompiler: DSCompiler?, + sizes: Array, + derivativesIndirection: Array, +): Array> { + if (parameters == 0 || order == 0) { + return arrayOf(arrayOf(intArrayOf(1, 0))) + } + + val vSize = valueCompiler!!.compositionIndirection.size + val dSize = derivativeCompiler!!.compositionIndirection.size + val compIndirection: Array?> = arrayOfNulls(vSize + dSize) + + // the composition rules from the value part can be reused as is + valueCompiler.compositionIndirection.copyInto(compIndirection, endIndex = vSize) + + // the composition rules for the derivative part are deduced by differentiation the rules from the + // underlying compiler once with respect to the parameter this compiler handles and the underlying one + // did not handle + + // the composition rules for the derivative part are deduced by differentiation the rules from the + // underlying compiler once with respect to the parameter this compiler handles and the underlying one did + // not handle + for (i in 0 until dSize) { + val row: List = buildList { + for (term in derivativeCompiler.compositionIndirection[i]) { + + // handle term p * f_k(g(x)) * g_l1(x) * g_l2(x) * ... * g_lp(x) + + // derive the first factor in the term: f_k with respect to new parameter + val derivedTermF = IntArray(term.size + 1) + derivedTermF[0] = term[0] // p + derivedTermF[1] = term[1] + 1 // f_(k+1) + val orders = IntArray(parameters) + orders[parameters - 1] = 1 + derivedTermF[term.size] = getPartialDerivativeIndex( + parameters, + order, + sizes, + *orders + ) // g_1 + + for (j in 2 until term.size) { + // convert the indices as the mapping for the current order is different from the mapping with one + // less order + derivedTermF[j] = convertIndex( + term[j], parameters, + derivativeCompiler.derivativesIndirection, + parameters, order, sizes + ) + } + + derivedTermF.sort(2, derivedTermF.size) + add(derivedTermF) + + // derive the various g_l + for (l in 2 until term.size) { + val derivedTermG = IntArray(term.size) + derivedTermG[0] = term[0] + derivedTermG[1] = term[1] + + for (j in 2 until term.size) { + // convert the indices as the mapping for the current order + // is different from the mapping with one less order + derivedTermG[j] = convertIndex( + term[j], + parameters, + derivativeCompiler.derivativesIndirection, + parameters, + order, + sizes, + ) + + if (j == l) { + // derive this term + derivativesIndirection[derivedTermG[j]].copyInto(orders, endIndex = parameters) + orders[parameters - 1]++ + + derivedTermG[j] = getPartialDerivativeIndex( + parameters, + order, + sizes, + *orders, + ) + } + } + + derivedTermG.sort(2, derivedTermG.size) + add(derivedTermG) + } + } + } + + // combine terms with similar derivation orders + val combined: List = buildList(row.size) { + for (j in row.indices) { + val termJ = row[j] + + if (termJ[0] > 0) { + (j + 1 until row.size).map { k -> row[k] }.forEach { termK -> + var equals = termJ.size == termK.size + var l = 1 + + while (equals && l < termJ.size) { + equals = equals and (termJ[l] == termK[l]) + ++l + } + + if (equals) { + // combine termJ and termK + termJ[0] += termK[0] + // make sure we will skip termK later on in the outer loop + termK[0] = 0 + } + } + + add(termJ) + } + } + } + + compIndirection[vSize + i] = combined.toTypedArray() + } + + return compIndirection as Array> +} + +/** + * Get the index of a partial derivative in an array. + * + * @param parameters number of free parameters. + * @param order derivation order. + * @param sizes sizes array. + * @param orders derivation orders with respect to each parameter (the length of this array must match the number of + * parameters). + * @return index of the partial derivative. + */ +private fun getPartialDerivativeIndex( + parameters: Int, + order: Int, + sizes: Array, + vararg orders: Int, +): Int { + + // the value is obtained by diving into the recursive Dan Kalman's structure + // this is theorem 2 of his paper, with recursion replaced by iteration + var index = 0 + var m = order + var ordersSum = 0 + + for (i in parameters - 1 downTo 0) { + // derivative order for current free parameter + var derivativeOrder = orders[i] + + // safety check + ordersSum += derivativeOrder + require(ordersSum <= order) { "number is too large: $ordersSum > $order" } + + while (derivativeOrder-- > 0) { + // as long as we differentiate according to current free parameter, + // we have to skip the value part and dive into the derivative part, + // so we add the size of the value part to the base index + index += sizes[i][m--] + } + } + + return index +} + +/** + * Convert an index from one (parameters, order) structure to another. + * + * @param index index of a partial derivative in source derivative structure. + * @param srcP number of free parameters in source derivative structure. + * @param srcDerivativesIndirection derivatives indirection array for the source derivative structure. + * @param destP number of free parameters in destination derivative structure. + * @param destO derivation order in destination derivative structure. + * @param destSizes sizes array for the destination derivative structure. + * @return index of the partial derivative with the *same* characteristics in destination derivative structure. + */ +private fun convertIndex( + index: Int, + srcP: Int, + srcDerivativesIndirection: Array, + destP: Int, + destO: Int, + destSizes: Array, +): Int { + val orders = IntArray(destP) + srcDerivativesIndirection[index].copyInto(orders, endIndex = min(srcP, destP)) + return getPartialDerivativeIndex(destP, destO, destSizes, *orders) +} 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 new file mode 100644 index 000000000..a1a6354f0 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DerivativeStructure.kt @@ -0,0 +1,186 @@ +/* + * 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.NumericAlgebra +import space.kscience.kmath.operations.Ring +import space.kscience.kmath.operations.ScaleOperations +import space.kscience.kmath.structures.MutableBuffer + +/** + * 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 internal constructor( + internal val derivativeAlgebra: DerivativeStructureRing, + internal val compiler: DSCompiler, +) where A : Ring, A : NumericAlgebra, A : ScaleOperations { + /** + * Combined array holding all values. + */ + internal var data: MutableBuffer = + derivativeAlgebra.bufferFactory(compiler.size) { derivativeAlgebra.algebra.zero } + + /** + * Build an instance with all values and derivatives set to 0. + * + * @param parameters number of free parameters. + * @param order derivation order. + */ + public constructor ( + derivativeAlgebra: DerivativeStructureRing, + parameters: Int, + order: Int, + ) : this( + derivativeAlgebra, + getCompiler(derivativeAlgebra.algebra, derivativeAlgebra.bufferFactory, parameters, order), + ) + + /** + * Build an instance representing a constant value. + * + * @param parameters number of free parameters. + * @param order derivation order. + * @param value value of the constant. + * @see DerivativeStructure + */ + public constructor ( + derivativeAlgebra: DerivativeStructureRing, + parameters: Int, + order: Int, + value: T, + ) : this( + derivativeAlgebra, + parameters, + order, + ) { + data[0] = value + } + + /** + * 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. + * + * @param parameters number of free parameters. + * @param order derivation order. + * @param index index of the variable (from 0 to `parameters - 1`). + * @param value value of the variable. + */ + public constructor ( + derivativeAlgebra: DerivativeStructureRing, + parameters: Int, + order: Int, + index: Int, + value: T, + ) : this(derivativeAlgebra, parameters, order, value) { + require(index < parameters) { "number is too large: $index >= $parameters" } + + if (order > 0) { + // the derivative of the variable with respect to itself is 1. + data[getCompiler(derivativeAlgebra.algebra, derivativeAlgebra.bufferFactory, index, order).size] = + derivativeAlgebra.algebra.one + } + } + + /** + * Build an instance from all its derivatives. + * + * @param parameters number of free parameters. + * @param order derivation order. + * @param derivatives derivatives sorted according to [DSCompiler.getPartialDerivativeIndex]. + */ + public constructor ( + derivativeAlgebra: DerivativeStructureRing, + parameters: Int, + order: Int, + vararg derivatives: T, + ) : this( + derivativeAlgebra, + parameters, + order, + ) { + require(derivatives.size == data.size) { "dimension mismatch: ${derivatives.size} and ${data.size}" } + data = derivativeAlgebra.bufferFactory(data.size) { derivatives[it] } + } + + /** + * Copy constructor. + * + * @param ds instance to copy. + */ + internal constructor(ds: DerivativeStructure) : this(ds.derivativeAlgebra, ds.compiler) { + this.data = ds.data.copy() + } + + /** + * 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() +} 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 new file mode 100644 index 000000000..f91fb55e8 --- /dev/null +++ b/kmath-core/src/commonMain/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpression.kt @@ -0,0 +1,332 @@ +/* + * 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.MutableBufferFactory +import space.kscience.kmath.structures.indices + +/** + * A class implementing both [DerivativeStructure] and [Symbol]. + */ +@UnstableKMathAPI +public class DerivativeStructureSymbol( + derivativeAlgebra: DerivativeStructureRing, + size: Int, + order: Int, + index: Int, + symbol: Symbol, + value: T, +) : Symbol by symbol, DerivativeStructure( + derivativeAlgebra, + size, + order, + index, + value +) where A : Ring, A : NumericAlgebra, A : ScaleOperations { + override fun toString(): String = symbol.toString() + override fun equals(other: Any?): Boolean = (other as? Symbol) == symbol + override fun hashCode(): Int = symbol.hashCode() +} + +/** + * 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( + public val algebra: A, + public val bufferFactory: MutableBufferFactory, + public val order: Int, + bindings: Map, +) : Ring>, ScaleOperations>, + NumericAlgebra>, + ExpressionAlgebra>, + NumbersAddOps> where A : Ring, A : NumericAlgebra, A : ScaleOperations { + public val numberOfVariables: Int = bindings.size + + override val zero: DerivativeStructure by lazy { + DerivativeStructure( + this, + numberOfVariables, + order, + ) + } + + override val one: DerivativeStructure by lazy { + DerivativeStructure( + this, + numberOfVariables, + order, + algebra.one, + ) + } + + override fun number(value: Number): DerivativeStructure = const(algebra.number(value)) + + private val variables: Map> = + bindings.entries.mapIndexed { index, (key, value) -> + key to DerivativeStructureSymbol( + this, + numberOfVariables, + order, + index, + key, + value, + ) + }.toMap() + + public override fun const(value: T): DerivativeStructure = + DerivativeStructure(this, numberOfVariables, order, value) + + 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()) + + override fun DerivativeStructure.unaryMinus(): DerivativeStructure { + val ds = DerivativeStructure(this@DerivativeStructureRing, compiler) + for (i in ds.data.indices) { + ds.data[i] = algebra { -data[i] } + } + return ds + } + + override fun add(left: DerivativeStructure, right: DerivativeStructure): DerivativeStructure { + left.compiler.checkCompatibility(right.compiler) + val ds = DerivativeStructure(left) + left.compiler.add(left.data, 0, right.data, 0, ds.data, 0) + return ds + } + + override fun scale(a: DerivativeStructure, value: Double): DerivativeStructure { + val ds = DerivativeStructure(a) + for (i in ds.data.indices) { + ds.data[i] = algebra { ds.data[i].times(value) } + } + return ds + } + + override fun multiply( + left: DerivativeStructure, + right: DerivativeStructure + ): DerivativeStructure { + left.compiler.checkCompatibility(right.compiler) + val result = DerivativeStructure(this, left.compiler) + left.compiler.multiply(left.data, 0, right.data, 0, result.data, 0) + return result + } + + override fun DerivativeStructure.minus(arg: DerivativeStructure): DerivativeStructure { + compiler.checkCompatibility(arg.compiler) + val ds = DerivativeStructure(this) + compiler.subtract(data, 0, arg.data, 0, ds.data, 0) + return ds + } + + override operator fun DerivativeStructure.plus(other: Number): DerivativeStructure { + val ds = DerivativeStructure(this) + ds.data[0] = algebra { ds.data[0] + number(other) } + return ds + } + + 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) + val result = DerivativeStructure(this, left.compiler) + left.compiler.divide(left.data, 0, right.data, 0, result.data, 0) + return result + } + + override fun sin(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.sin(arg.data, 0, result.data, 0) + return result + } + + override fun cos(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.cos(arg.data, 0, result.data, 0) + return result + } + + override fun tan(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.tan(arg.data, 0, result.data, 0) + return result + } + + override fun asin(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.asin(arg.data, 0, result.data, 0) + return result + } + + override fun acos(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.acos(arg.data, 0, result.data, 0) + return result + } + + override fun atan(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.atan(arg.data, 0, result.data, 0) + return result + } + + override fun sinh(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.sinh(arg.data, 0, result.data, 0) + return result + } + + override fun cosh(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.cosh(arg.data, 0, result.data, 0) + return result + } + + override fun tanh(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.tanh(arg.data, 0, result.data, 0) + return result + } + + override fun asinh(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.asinh(arg.data, 0, result.data, 0) + return result + } + + override fun acosh(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.acosh(arg.data, 0, result.data, 0) + return result + } + + override fun atanh(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.atanh(arg.data, 0, result.data, 0) + return result + } + + override fun power(arg: DerivativeStructure, pow: Number): DerivativeStructure = when (pow) { + is Int -> { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.pow(arg.data, 0, pow, result.data, 0) + result + } + else -> { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.pow(arg.data, 0, pow.toDouble(), result.data, 0) + result + } + } + + override fun sqrt(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.sqrt(arg.data, 0, result.data, 0) + return result + } + + public fun power(arg: DerivativeStructure, pow: DerivativeStructure): DerivativeStructure { + arg.compiler.checkCompatibility(pow.compiler) + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.pow(arg.data, 0, pow.data, 0, result.data, 0) + return result + } + + override fun exp(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.exp(arg.data, 0, result.data, 0) + return result + } + + override fun ln(arg: DerivativeStructure): DerivativeStructure { + val result = DerivativeStructure(this, arg.compiler) + arg.compiler.ln(arg.data, 0, result.data, 0) + return result + } +} + +@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 new file mode 100644 index 000000000..429fe310b --- /dev/null +++ b/kmath-core/src/commonTest/kotlin/space/kscience/kmath/expressions/DerivativeStructureExpressionTest.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, + block: DerivativeStructureField.() -> Unit, +) { + contract { callsInPlace(block, InvocationKind.EXACTLY_ONCE) } + DerivativeStructureField(DoubleField, ::DoubleBuffer, order, mapOf(*parameters)).block() +} + +internal class AutoDiffTest { + private val x by symbol + private val y by symbol + + @Test + fun derivativeStructureFieldTest() { + 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 autoDifTest() { + val f = DerivativeStructureFieldExpression(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)) + } +} diff --git a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/internal/InternalUtils.kt b/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/internal/InternalUtils.kt index 3997a77b3..71fd15fe6 100644 --- a/kmath-stat/src/commonMain/kotlin/space/kscience/kmath/internal/InternalUtils.kt +++ b/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