commons-math3-3.6.1-类-DSCompiler-中英对照文档及源码赏析.md
摘要：中英对照文档、源码赏析、org.apache.commons.math3.analysis.differentiation.DSCompiler、DSCompiler、commons-math3、3.6.1、commons-math3-3.6.1、java
完整中文文档、中英对照文档下载请移步：commons-math3-中文文档、中英对照文档-CSDN下载
1. 开源组件说明
jar包名称：commons-math3-3.6.1.jar
Maven 依赖：
<dependency>
<groupId>org.apache.commons</groupId>
<artifactId>commons-math3</artifactId>
<version>3.6.1</version>
</dependency>
完整中文文档、中英对照文档下载请移步：commons-math3-中文文档、中英对照文档-CSDN下载
2. 类：org.apache.commons.math3.analysis.differentiation.DSCompiler
2.1. DSCompiler 中英对照文档
org.apache.commons.math3.analysis.differentiation
类 DSCompiler
- java.lang.Object
-
- org.apache.commons.math3.analysis.differentiation.DSCompiler
-
public class DSCompiler extends Object
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, Mathematics Magazine, vol. 75, no. 3, June 2002. However, in order to avoid performances bottlenecks, the recursive rules are "compiled" once in an unfold 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 theur 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
getPartialDerivativeIndexandgetPartialDerivativeOrders(int). The value is guaranteed to be stored as the first element (i.e. thegetPartialDerivativeIndexmethod returns 0 when called with 0 for all derivation orders andgetPartialDerivativeOrdersreturns 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.
This class is mainly used as the engine for scalar variable
DerivativeStructure. It can also be used directly to hold several variables in arrays for more complex data structures. User can for example store a vector of n variables depending on three x, y and z free parameters in one array as follows:类持有“编译”计算规则的衍生结构。此类实现丹卡尔曼纸双重递归多变量自动分化，数学杂志，Vol的计算规则。75，没有。3，2002年6月。然而，为了避免性能瓶颈，递归规则以展开形式“编译”一次。此类递归展开并将计算规则存储为具有预计算机间接阵列的简单环路。此类将所有导数计算映射到包含值和部分导数的单尺寸阵列中。该课程不持有这些阵列，这些阵列仍然受到来电者的责任。对于每种免费参数和派生顺序的每个组合，只有一个编译器是必要的，并且该编译器将用于对提供给它的所有阵列执行计算，这可以代表数百或数千个不同的参数与所有固定部分衍生物一起保持在一起。编译器操作的阵列仅包含部分衍生物以及第0衍生物，我。e。价值。部分导数以编译器特定顺序存储，可以使用FactPartialDerivativeIndex和GetPartialDivationers（int）来检索。该值保证将存储为第一元素（i。e。GetPartialDerivativeInd方法返回0当调用0对于所有导出订单时，GetPartialDerivorationders返回填充有0的数组，其中0作为索引调用0）。请注意，订购更改为参数和派生顺序。例如，当派生顺序设置为1（在这种情况下，阵列有三个元素：F，DF / DX和DF / Dy，DF / Dy存储在索引2中。如果推导顺序设置为2，则DF / DY将存储在索引3（在这种情况下，阵列具有六个元素：F，DF / DX，DF / DXDX，DF / DY，DF / DXDY和DF / DYDY）。鉴于此结构，用户可以执行一些简单的操作，例如添加，减去或乘以常量并通过自己否定元素，知道它们是否要突变它们的数组或创建新数组。编译器不提供这些简单的操作。编译器仅在几个阵列之间提供更复杂的操作。此类主要用作标量变量衍生体的引擎。它也可以直接用于在阵列中保持多个变量以获得更复杂的数据结构。用户可以例如存储n个变量的向量，具体取决于一个阵列中的三个x，y和z空闲参数，如下所示：
// parameter 0 is x, parameter 1 is y, parameter 2 is z int parameters = 3; DSCompiler compiler = DSCompiler.getCompiler(parameters, order); int size = compiler.getSize(); // pack all elements in a single array double[] array = new double[n * size]; for (int i = 0; i < n; ++i) { // we know value is guaranteed to be the first element array[i * size] = v[i]; // we don't know where first derivatives are stored, so we ask the compiler array[i * size + compiler.getPartialDerivativeIndex(1, 0, 0) = dvOnDx[i][0]; array[i * size + compiler.getPartialDerivativeIndex(0, 1, 0) = dvOnDy[i][0]; array[i * size + compiler.getPartialDerivativeIndex(0, 0, 1) = dvOnDz[i][0]; // we let all higher order derivatives set to 0 }Then in another function, user can perform some operations on all elements stored in the single array, such as a simple product of all variables:
然后在另一个函数中，用户可以对存储在单个阵列中的所有元素的一些操作，例如所有变量的简单产品：
// compute the product of all elements double[] product = new double[size]; prod[0] = 1.0; for (int i = 0; i < n; ++i) { double[] tmp = product.clone(); compiler.multiply(tmp, 0, array, i * size, product, 0); } // value double p = product[0]; // first derivatives double dPdX = product[compiler.getPartialDerivativeIndex(1, 0, 0)]; double dPdY = product[compiler.getPartialDerivativeIndex(0, 1, 0)]; double dPdZ = product[compiler.getPartialDerivativeIndex(0, 0, 1)]; // cross derivatives (assuming order was at least 2) double dPdXdX = product[compiler.getPartialDerivativeIndex(2, 0, 0)]; double dPdXdY = product[compiler.getPartialDerivativeIndex(1, 1, 0)]; double dPdXdZ = product[compiler.getPartialDerivativeIndex(1, 0, 1)]; double dPdYdY = product[compiler.getPartialDerivativeIndex(0, 2, 0)]; double dPdYdZ = product[compiler.getPartialDerivativeIndex(0, 1, 1)]; double dPdZdZ = product[compiler.getPartialDerivativeIndex(0, 0, 2)];
-
从以下版本开始:
-
3.1
3.1
另请参阅:
-
DerivativeStructure
衍生体系
-
3.1
-
-
方法概要
Methods 限定符和类型 方法和说明 voidacos(double[] operand, int operandOffset, double[] result, int resultOffset)Compute arc cosine of a derivative structure.
计算衍生结构的弧余弦。
voidacosh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute inverse hyperbolic cosine of a derivative structure.
计算衍生结构的反向双曲余弦。
voidadd(double[] lhs, int lhsOffset, double[] rhs, int rhsOffset, double[] result, int resultOffset)Perform addition of two derivative structures.
进行两个衍生结构的添加。
voidasin(double[] operand, int operandOffset, double[] result, int resultOffset)Compute arc sine of a derivative structure.
计算衍生结构的弧正弦。
voidasinh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute inverse hyperbolic sine of a derivative structure.
计算衍生结构的逆双曲正弦。
voidatan(double[] operand, int operandOffset, double[] result, int resultOffset)Compute arc tangent of a derivative structure.
计算衍生结构的弧形切线。
voidatan2(double[] y, int yOffset, double[] x, int xOffset, double[] result, int resultOffset)Compute two arguments arc tangent of a derivative structure.
计算衍生结构的两个参数arc切线。
voidatanh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute inverse hyperbolic tangent of a derivative structure.
计算衍生结构的反向双曲线切线。
voidcheckCompatibility(DSCompiler compiler)Check rules set compatibility.
检查规则集兼容性。
voidcompose(double[] operand, int operandOffset, double[] f, double[] result, int resultOffset)Compute composition of a derivative structure by a function.
通过功能计算衍生结构的组成。
voidcos(double[] operand, int operandOffset, double[] result, int resultOffset)Compute cosine of a derivative structure.
计算衍生物结构的余弦。
voidcosh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute hyperbolic cosine of a derivative structure.
计算衍生结构的双曲线余弦。
voiddivide(double[] lhs, int lhsOffset, double[] rhs, int rhsOffset, double[] result, int resultOffset)Perform division of two derivative structures.
执行两个衍生结构的划分。
voidexp(double[] operand, int operandOffset, double[] result, int resultOffset)Compute exponential of a derivative structure.
计算衍生结构的指数。
voidexpm1(double[] operand, int operandOffset, double[] result, int resultOffset)Compute exp(x) - 1 of a derivative structure.
计算exp（x） - 衍生结构的1。
static DSCompilergetCompiler(int parameters, int order)Get the compiler for number of free parameters and order.
获取编译器的免费参数和订单。
intgetFreeParameters()Get the number of free parameters.
获取免费参数的数量。
intgetOrder()Get the derivation order.
获取派生顺序。
intgetPartialDerivativeIndex(int... orders)Get the index of a partial derivative in the array.
获取数组中部分导数的索引。
int[]getPartialDerivativeOrders(int index)Get the derivation orders for a specific index in the array.
获取阵列中特定索引的派生顺序。
intgetSize()Get the array size required for holding partial derivatives data.
获取保持部分衍生产品数据所需的数组大小。
voidlinearCombination(double a1, double[] c1, int offset1, double a2, double[] c2, int offset2, double[] result, int resultOffset)Compute linear combination.
计算线性组合。
voidlinearCombination(double a1, double[] c1, int offset1, double a2, double[] c2, int offset2, double a3, double[] c3, int offset3, double[] result, int resultOffset)Compute linear combination.
计算线性组合。
voidlinearCombination(double a1, double[] c1, int offset1, double a2, double[] c2, int offset2, double a3, double[] c3, int offset3, double a4, double[] c4, int offset4, double[] result, int resultOffset)Compute linear combination.
计算线性组合。
voidlog(double[] operand, int operandOffset, double[] result, int resultOffset)Compute natural logarithm of a derivative structure.
计算衍生结构的自然对数。
voidlog10(double[] operand, int operandOffset, double[] result, int resultOffset)Computes base 10 logarithm of a derivative structure.
计算衍生结构的基础10对数。
voidlog1p(double[] operand, int operandOffset, double[] result, int resultOffset)Computes shifted logarithm of a derivative structure.
计算衍生结构的移位对数。
voidmultiply(double[] lhs, int lhsOffset, double[] rhs, int rhsOffset, double[] result, int resultOffset)Perform multiplication of two derivative structures.
执行两个衍生结构的乘法。
voidpow(double[] x, int xOffset, double[] y, int yOffset, double[] result, int resultOffset)Compute power of a derivative structure.
计算衍生结构的功率。
voidpow(double[] operand, int operandOffset, double p, double[] result, int resultOffset)Compute power of a derivative structure.
计算衍生结构的功率。
voidpow(double[] operand, int operandOffset, int n, double[] result, int resultOffset)Compute integer power of a derivative structure.
计算衍生结构的整数功率。
voidpow(double a, double[] operand, int operandOffset, double[] result, int resultOffset)Compute power of a double to a derivative structure.
计算双倍到衍生结构的功率。
voidremainder(double[] lhs, int lhsOffset, double[] rhs, int rhsOffset, double[] result, int resultOffset)Perform remainder of two derivative structures.
执行剩余的两个衍生结构。
voidrootN(double[] operand, int operandOffset, int n, double[] result, int resultOffset)Compute n th root of a derivative structure.
计算衍生结构的第n个根。
voidsin(double[] operand, int operandOffset, double[] result, int resultOffset)Compute sine of a derivative structure.
计算衍生结构的正弦。
voidsinh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute hyperbolic sine of a derivative structure.
计算衍生结构的双曲线正弦。
voidsubtract(double[] lhs, int lhsOffset, double[] rhs, int rhsOffset, double[] result, int resultOffset)Perform subtraction of two derivative structures.
执行两种衍生结构的减法。
voidtan(double[] operand, int operandOffset, double[] result, int resultOffset)Compute tangent of a derivative structure.
计算衍生结构的切线。
voidtanh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute hyperbolic tangent of a derivative structure.
计算衍生结构的双曲线切线。
doubletaylor(double[] ds, int dsOffset, double... delta)Evaluate Taylor expansion of a derivative structure.
评估达勒扩大衍生物结构。
-
-
-
方法详细说明
-
getCompiler
public static DSCompiler getCompiler(int parameters, int order) throws NumberIsTooLargeException
Get the compiler for number of free parameters and order.
获取编译器的免费参数和订单。
-
参数:
-
parameters- number of free parameters
参数 - 免费参数数
-
order- derivation order
订单 - 衍生顺序
返回:
-
cached rules set
缓存规则集
抛出:
-
NumberIsTooLargeException- if order is too large
NumberStoolarGeexception - 如果订单太大
-
-
getPartialDerivativeIndex
public int getPartialDerivativeIndex(int... orders) throws DimensionMismatchException, NumberIsTooLargeExceptionGet 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
getSize()- 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
free parameter, 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
derivation orderis 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(int)
获取数组中部分导数的索引。如果所有订单设置为0，则返回0th阶数，这是函数的值。衍生物的指数在0和getsize之间（） - 1.它们的特定顺序是针对给定编译器修复的，但否则不公开指定。然而，有一些简单的案例有保证指数：第0阶衍生物的索引始终为0如果只有1个免费参数，则衍生物以越来越多的推导顺序排序（即索引0，df / dp在索引1处的f。，D2F / DP2在索引2 ... DKF / DPK处Idex K），如果推导顺序为1，则衍生物以增加的自由参数顺序（即索引0，DF / DX1在索引1，DF处）进行排序。/ dx2在索引2 ... df / dxk处index k），所有其他情况都不公开指定此方法是方法getPartialDeriverationders（int）的倒数
-
参数:
-
orders- derivation orders with respect to each parameter
订单 - 派生关于每个参数的订单
返回:
-
index of the partial derivative
部分衍生物的指数
抛出:
-
DimensionMismatchException- if the numbers of parameters does not match the instance
DimensionMismatchException - 如果参数的数量不匹配实例
-
NumberIsTooLargeException- if sum of derivation orders is larger than the instance limits
NumberStoolarGeexception - 如果派生订单的总和大于实例限制
另请参阅:
-
getPartialDerivativeOrders(int)
getPartialDerivorationders（int）
-
getPartialDerivativeOrders
public int[] getPartialDerivativeOrders(int index)
Get the derivation orders for a specific index in the array.This method is the inverse of
getPartialDerivativeIndex(int...).
获取阵列中特定索引的派生顺序。此方法是GetPartialDerivativeIndIndex（int ...）的倒数。
-
参数:
-
index- of the partial derivative
索引 - 部分衍生物
返回:
-
orders derivation orders with respect to each parameter
订单导出关于每个参数的订单
另请参阅:
-
getPartialDerivativeIndex(int...)
GetPartialDerivativeIndex（int ...）
-
-
getFreeParameters
public int getFreeParameters()
Get the number of free parameters.
获取免费参数的数量。
-
返回:
-
number of free parameters
自由参数数量
-
number of free parameters
-
getOrder
public int getOrder()
Get the derivation order.
获取派生顺序。
-
返回:
-
derivation order
推导令
-
derivation order
-
getSize
public int getSize()
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.
获取保持部分衍生产品数据所需的数组大小。此数字包括单个0阶衍生元素，保证存储在阵列的第一元素中。
-
返回:
-
array size required for holding partial derivatives data
保持部分衍生物数据所需的数组大小
-
array size required for holding partial derivatives data
-
linearCombination
public void linearCombination(double a1, double[] c1, int offset1, double a2, double[] c2, int offset2, double[] result, int resultOffset)Compute linear combination. The derivative structure built will be a1 * ds1 + a2 * ds2
计算线性组合。构建的衍生结构将是A1 * DS1 + A2 * DS2
-
参数:
-
a1- first scale factor
A1 - 第一规模因子
-
c1- first base (unscaled) component
C1 - 第一基础（未坐标）组件
-
offset1- offset of first operand in its array
offset1 - 第一个操作数在其数组中的偏移量
-
a2- second scale factor
A2 - 二级规模因子
-
c2- second base (unscaled) component
C2 - 第二基础（未坐标）组件
-
offset2- offset of second operand in its array
offset2 - 第二个操作数在其数组中的偏移量
-
result- array where result must be stored (it may be one of the input arrays)
结果 - 必须存储结果的阵列（它可以是输入阵列之一）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
linearCombination
public void linearCombination(double a1, double[] c1, int offset1, double a2, double[] c2, int offset2, double a3, double[] c3, int offset3, double[] result, int resultOffset)Compute linear combination. The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
计算线性组合。构建的衍生结构将是A1 * DS1 + A2 * DS2 + A3 * DS3 + A4 * DS4
-
参数:
-
a1- first scale factor
A1 - 第一规模因子
-
c1- first base (unscaled) component
C1 - 第一基础（未坐标）组件
-
offset1- offset of first operand in its array
offset1 - 第一个操作数在其数组中的偏移量
-
a2- second scale factor
A2 - 二级规模因子
-
c2- second base (unscaled) component
C2 - 第二基础（未坐标）组件
-
offset2- offset of second operand in its array
offset2 - 第二个操作数在其数组中的偏移量
-
a3- third scale factor
A3 - 第三比例因子
-
c3- third base (unscaled) component
C3 - 第三基础（未坐标）组件
-
offset3- offset of third operand in its array
offset3 - 第三操作数在其数组中的偏移量
-
result- array where result must be stored (it may be one of the input arrays)
结果 - 必须存储结果的阵列（它可以是输入阵列之一）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
linearCombination
public void linearCombination(double a1, double[] c1, int offset1, double a2, double[] c2, int offset2, double a3, double[] c3, int offset3, double a4, double[] c4, int offset4, double[] result, int resultOffset)Compute linear combination. The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
计算线性组合。构建的衍生结构将是A1 * DS1 + A2 * DS2 + A3 * DS3 + A4 * DS4
-
参数:
-
a1- first scale factor
A1 - 第一规模因子
-
c1- first base (unscaled) component
C1 - 第一基础（未坐标）组件
-
offset1- offset of first operand in its array
offset1 - 第一个操作数在其数组中的偏移量
-
a2- second scale factor
A2 - 二级规模因子
-
c2- second base (unscaled) component
C2 - 第二基础（未坐标）组件
-
offset2- offset of second operand in its array
offset2 - 第二个操作数在其数组中的偏移量
-
a3- third scale factor
A3 - 第三比例因子
-
c3- third base (unscaled) component
C3 - 第三基础（未坐标）组件
-
offset3- offset of third operand in its array
offset3 - 第三操作数在其数组中的偏移量
-
a4- fourth scale factor
A4 - 第四比例因子
-
c4- fourth base (unscaled) component
C4 - 第四基（未坐标）组件
-
offset4- offset of fourth operand in its array
offset4 - 第四个操作数在其数组中的偏移量
-
result- array where result must be stored (it may be one of the input arrays)
结果 - 必须存储结果的阵列（它可以是输入阵列之一）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
add
public void add(double[] lhs, int lhsOffset, double[] rhs, int rhsOffset, double[] result, int resultOffset)Perform addition of two derivative structures.
进行两个衍生结构的添加。
-
参数:
-
lhs- array holding left hand side of addition
LHS - 阵列握住左侧
-
lhsOffset- offset of the left hand side in its array
lhsoffset - 左侧的阵列偏移
-
rhs- array right hand side of addition
RHS - 阵列右侧添加
-
rhsOffset- offset of the right hand side in its array
RHSOFFSET - 右侧偏移阵列
-
result- array where result must be stored (it may be one of the input arrays)
结果 - 必须存储结果的阵列（它可以是输入阵列之一）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
subtract
public void subtract(double[] lhs, int lhsOffset, double[] rhs, int rhsOffset, double[] result, int resultOffset)Perform subtraction of two derivative structures.
执行两种衍生结构的减法。
-
参数:
-
lhs- array holding left hand side of subtraction
LHS - 阵列握住减法的左侧
-
lhsOffset- offset of the left hand side in its array
lhsoffset - 左侧的阵列偏移
-
rhs- array right hand side of subtraction
RHS - 阵列减法的右侧
-
rhsOffset- offset of the right hand side in its array
RHSOFFSET - 右侧偏移阵列
-
result- array where result must be stored (it may be one of the input arrays)
结果 - 必须存储结果的阵列（它可以是输入阵列之一）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
multiply
public void multiply(double[] lhs, int lhsOffset, double[] rhs, int rhsOffset, double[] result, int resultOffset)Perform multiplication of two derivative structures.
执行两个衍生结构的乘法。
-
参数:
-
lhs- array holding left hand side of multiplication
LHS - 阵列握住左侧的乘法
-
lhsOffset- offset of the left hand side in its array
lhsoffset - 左侧的阵列偏移
-
rhs- array right hand side of multiplication
RHS - 阵列右手乘以乘法
-
rhsOffset- offset of the right hand side in its array
RHSOFFSET - 右侧偏移阵列
-
result- array where result must be stored (for multiplication the result array cannot be one of the input arrays)
结果 - 必须存储结果的阵列（对于乘法，结果数组不能是输入阵列之一）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
divide
public void divide(double[] lhs, int lhsOffset, double[] rhs, int rhsOffset, double[] result, int resultOffset)Perform division of two derivative structures.
执行两个衍生结构的划分。
-
参数:
-
lhs- array holding left hand side of division
LHS - 阵列握住左侧部门
-
lhsOffset- offset of the left hand side in its array
lhsoffset - 左侧的阵列偏移
-
rhs- array right hand side of division
RHS - 阵列右手边
-
rhsOffset- offset of the right hand side in its array
RHSOFFSET - 右侧偏移阵列
-
result- array where result must be stored (for division the result array cannot be one of the input arrays)
结果 - 必须存储结果的阵列（对于划分结果阵列不能是输入阵列之一）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
remainder
public void remainder(double[] lhs, int lhsOffset, double[] rhs, int rhsOffset, double[] result, int resultOffset)Perform remainder of two derivative structures.
执行剩余的两个衍生结构。
-
参数:
-
lhs- array holding left hand side of remainder
LHS - 阵列握住余下的左侧
-
lhsOffset- offset of the left hand side in its array
lhsoffset - 左侧的阵列偏移
-
rhs- array right hand side of remainder
RHS - 阵列右侧的剩余部分
-
rhsOffset- offset of the right hand side in its array
RHSOFFSET - 右侧偏移阵列
-
result- array where result must be stored (it may be one of the input arrays)
结果 - 必须存储结果的阵列（它可以是输入阵列之一）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
pow
public void pow(double a, double[] operand, int operandOffset, double[] result, int resultOffset)Compute power of a double to a derivative structure.
计算双倍到衍生结构的功率。
-
参数:
-
a- number to exponentiate
a - 代表的号码
-
operand- array holding the power
操作数 - 阵列持有电源
-
operandOffset- offset of the power in its array
OutmandOffset - 其数组中电源的偏移量
-
result- array where result must be stored (for power the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于Power结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
从以下版本开始:
-
3.3
3.3
-
-
pow
public void pow(double[] operand, int operandOffset, double p, double[] result, int resultOffset)Compute power of a derivative structure.
计算衍生结构的功率。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
p- power to apply
P - 申请权力
-
result- array where result must be stored (for power the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于Power结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
pow
public void pow(double[] operand, int operandOffset, int n, double[] result, int resultOffset)Compute integer power of a derivative structure.
计算衍生结构的整数功率。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
n- power to apply
n - 申请的力量
-
result- array where result must be stored (for power the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于Power结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
pow
public void pow(double[] x, int xOffset, double[] y, int yOffset, double[] result, int resultOffset)Compute power of a derivative structure.
计算衍生结构的功率。
-
参数:
-
x- array holding the base
X - 阵列保持基座
-
xOffset- offset of the base in its array
xoffset - 基数在其数组中的偏移量
-
y- array holding the exponent
Y - 阵列持有指数
-
yOffset- offset of the exponent in its array
yoffset - exconent在其数组中的偏移量
-
result- array where result must be stored (for power the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于Power结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
rootN
public void rootN(double[] operand, int operandOffset, int n, double[] result, int resultOffset)Compute n th root of a derivative structure.
计算衍生结构的第n个根。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
n- order of the root
n - root的顺序
-
result- array where result must be stored (for n th root the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于第n个根结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
exp
public void exp(double[] operand, int operandOffset, double[] result, int resultOffset)Compute exponential of a derivative structure.
计算衍生结构的指数。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for exponential the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于指数结果数组不能是输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
expm1
public void expm1(double[] operand, int operandOffset, double[] result, int resultOffset)Compute exp(x) - 1 of a derivative structure.
计算exp（x） - 衍生结构的1。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for exponential the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于指数结果数组不能是输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
log
public void log(double[] operand, int operandOffset, double[] result, int resultOffset)Compute natural logarithm of a derivative structure.
计算衍生结构的自然对数。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for logarithm the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于对数而结果数组不能是输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
log1p
public void log1p(double[] operand, int operandOffset, double[] result, int resultOffset)Computes shifted logarithm of a derivative structure.
计算衍生结构的移位对数。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for shifted logarithm the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于移位的对数，结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
log10
public void log10(double[] operand, int operandOffset, double[] result, int resultOffset)Computes base 10 logarithm of a derivative structure.
计算衍生结构的基础10对数。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for base 10 logarithm the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于基本10对数，结果阵列不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
cos
public void cos(double[] operand, int operandOffset, double[] result, int resultOffset)Compute cosine of a derivative structure.
计算衍生物结构的余弦。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for cosine the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于余弦，结果数组不能是输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
sin
public void sin(double[] operand, int operandOffset, double[] result, int resultOffset)Compute sine of a derivative structure.
计算衍生结构的正弦。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for sine the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于SINE结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
tan
public void tan(double[] operand, int operandOffset, double[] result, int resultOffset)Compute tangent of a derivative structure.
计算衍生结构的切线。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for tangent the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于切线，结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
acos
public void acos(double[] operand, int operandOffset, double[] result, int resultOffset)Compute arc cosine of a derivative structure.
计算衍生结构的弧余弦。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for arc cosine the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于弧余弦，结果阵列不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
asin
public void asin(double[] operand, int operandOffset, double[] result, int resultOffset)Compute arc sine of a derivative structure.
计算衍生结构的弧正弦。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for arc sine the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于Arc Singe结果阵列不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
atan
public void atan(double[] operand, int operandOffset, double[] result, int resultOffset)Compute arc tangent of a derivative structure.
计算衍生结构的弧形切线。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for arc tangent the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于ARC切线而结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
atan2
public void atan2(double[] y, int yOffset, double[] x, int xOffset, double[] result, int resultOffset)Compute two arguments arc tangent of a derivative structure.
计算衍生结构的两个参数arc切线。
-
参数:
-
y- array holding the first operand
Y - 阵列持有第一个操作数
-
yOffset- offset of the first operand in its array
Yoffset - 第一个操作数在其数组中的偏移量
-
x- array holding the second operand
X - 阵列持有第二个操作数
-
xOffset- offset of the second operand in its array
xoffset - 第二个操作数在其数组中的偏移量
-
result- array where result must be stored (for two arguments arc tangent the result array cannot be the input array)
结果 - 阵列必须存储的位置（对于两个参数arc切线而结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
cosh
public void cosh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute hyperbolic cosine of a derivative structure.
计算衍生结构的双曲线余弦。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for hyperbolic cosine the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于双曲线余弦，结果阵列不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
sinh
public void sinh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute hyperbolic sine of a derivative structure.
计算衍生结构的双曲线正弦。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for hyperbolic sine the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于双曲正弦，结果阵列不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
tanh
public void tanh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute hyperbolic tangent of a derivative structure.
计算衍生结构的双曲线切线。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for hyperbolic tangent the result array cannot be the input array)
结果 - 阵列必须存储的结果（对于双曲线切线，结果阵列不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
acosh
public void acosh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute inverse hyperbolic cosine of a derivative structure.
计算衍生结构的反向双曲余弦。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for inverse hyperbolic cosine the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于逆双曲余弦，结果阵列不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
asinh
public void asinh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute inverse hyperbolic sine of a derivative structure.
计算衍生结构的逆双曲正弦。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for inverse hyperbolic sine the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于逆双曲正弦，结果阵列不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
atanh
public void atanh(double[] operand, int operandOffset, double[] result, int resultOffset)Compute inverse hyperbolic tangent of a derivative structure.
计算衍生结构的反向双曲线切线。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
result- array where result must be stored (for inverse hyperbolic tangent the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于反态双曲线切线，结果数组不能为输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
compose
public void compose(double[] operand, int operandOffset, double[] f, double[] result, int resultOffset)Compute composition of a derivative structure by a function.
通过功能计算衍生结构的组成。
-
参数:
-
operand- array holding the operand
操作数 - 持有操作数
-
operandOffset- offset of the operand in its array
OutmandOffset - 操作数在其数组中的偏移量
-
f- array of value and derivatives of the function at the current point (i.e. atoperand[operandOffset]).
f - 当前点的函数的值阵列和衍生物（即，在操作数[OperandOffset]）。
-
result- array where result must be stored (for composition the result array cannot be the input array)
结果 - 必须存储结果的阵列（对于组合而结果数组不能是输入数组）
-
resultOffset- offset of the result in its array
结果OFFSET - 其数组中结果的偏移量
-
-
taylor
public double taylor(double[] ds, int dsOffset, double... delta) throws MathArithmeticExceptionEvaluate Taylor expansion of a derivative structure.
评估达勒扩大衍生物结构。
-
参数:
-
ds- array holding the derivative structure
DS - 阵列保持衍生结构
-
dsOffset- offset of the derivative structure in its array
Dsoffset - 衍生结构在其数组中的偏移量
-
delta- parameters offsets (Δx, Δy, ...)
Delta - 参数偏移（ΔX，ΔY，......）
返回:
-
value of the Taylor expansion at x + Δx, y + Δy, ...
X +ΔX，Y +ΔY的泰勒膨胀的值，......
抛出:
-
MathArithmeticException- if factorials becomes too large
matharithmeticexception - 如果阶乘变得太大
-
-
checkCompatibility
public void checkCompatibility(DSCompiler compiler) throws DimensionMismatchException
Check rules set compatibility.
检查规则集兼容性。
-
参数:
-
compiler- other compiler to check against instance
编译器 - 其他编译器要检查实例
抛出:
-
DimensionMismatchException- if number of free parameters or orders are inconsistent
DimensionMismatchException - 如果免费参数或订单的数量不一致
-
-
-
2.2. DSCompiler 源码赏析
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.differentiation;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.concurrent.atomic.AtomicReference;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MathArithmeticException;
import org.apache.commons.math3.exception.MathInternalError;
import org.apache.commons.math3.exception.NotPositiveException;
import org.apache.commons.math3.exception.NumberIsTooLargeException;
import org.apache.commons.math3.util.CombinatoricsUtils;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;
/** Class holding "compiled" computation rules for derivative structures.
* <p>This class implements the computation rules described in Dan Kalman's paper <a
* href="http://www1.american.edu/cas/mathstat/People/kalman/pdffiles/mmgautodiff.pdf">Doubly
* Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75,
* no. 3, June 2002. However, in order to avoid performances bottlenecks, the recursive
* rules are "compiled" once in an unfold form. This class does this recursion unrolling
* and stores the computation rules as simple loops with pre-computed indirection arrays.</p>
* <p>
* 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 theur partial derivatives.
* </p>
* <p>
* The arrays on which compilers operate contain only the partial derivatives together
* with the 0<sup>th</sup> derivative, i.e. the value. The partial derivatives are stored in
* a compiler-specific order, which can be retrieved using methods {@link
* #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} and {@link
* #getPartialDerivativeOrders(int)}. The value is guaranteed to be stored as the first element
* (i.e. the {@link #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} method returns
* 0 when called with 0 for all derivation orders and {@link #getPartialDerivativeOrders(int)
* getPartialDerivativeOrders} returns an array filled with 0 when called with 0 as the index).
* </p>
* <p>
* 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).
* </p>
* <p>
* 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.
* </p>
* <p>This class is mainly used as the engine for scalar variable {@link DerivativeStructure}.
* It can also be used directly to hold several variables in arrays for more complex data
* structures. User can for example store a vector of n variables depending on three x, y
* and z free parameters in one array as follows:</p> <pre>
* // parameter 0 is x, parameter 1 is y, parameter 2 is z
* int parameters = 3;
* DSCompiler compiler = DSCompiler.getCompiler(parameters, order);
* int size = compiler.getSize();
*
* // pack all elements in a single array
* double[] array = new double[n * size];
* for (int i = 0; i < n; ++i) {
*
* // we know value is guaranteed to be the first element
* array[i * size] = v[i];
*
* // we don't know where first derivatives are stored, so we ask the compiler
* array[i * size + compiler.getPartialDerivativeIndex(1, 0, 0) = dvOnDx[i][0];
* array[i * size + compiler.getPartialDerivativeIndex(0, 1, 0) = dvOnDy[i][0];
* array[i * size + compiler.getPartialDerivativeIndex(0, 0, 1) = dvOnDz[i][0];
*
* // we let all higher order derivatives set to 0
*
* }
* </pre>
* <p>Then in another function, user can perform some operations on all elements stored
* in the single array, such as a simple product of all variables:</p> <pre>
* // compute the product of all elements
* double[] product = new double[size];
* prod[0] = 1.0;
* for (int i = 0; i < n; ++i) {
* double[] tmp = product.clone();
* compiler.multiply(tmp, 0, array, i * size, product, 0);
* }
*
* // value
* double p = product[0];
*
* // first derivatives
* double dPdX = product[compiler.getPartialDerivativeIndex(1, 0, 0)];
* double dPdY = product[compiler.getPartialDerivativeIndex(0, 1, 0)];
* double dPdZ = product[compiler.getPartialDerivativeIndex(0, 0, 1)];
*
* // cross derivatives (assuming order was at least 2)
* double dPdXdX = product[compiler.getPartialDerivativeIndex(2, 0, 0)];
* double dPdXdY = product[compiler.getPartialDerivativeIndex(1, 1, 0)];
* double dPdXdZ = product[compiler.getPartialDerivativeIndex(1, 0, 1)];
* double dPdYdY = product[compiler.getPartialDerivativeIndex(0, 2, 0)];
* double dPdYdZ = product[compiler.getPartialDerivativeIndex(0, 1, 1)];
* double dPdZdZ = product[compiler.getPartialDerivativeIndex(0, 0, 2)];
* </pre>
* @see DerivativeStructure
* @since 3.1
*/
public class DSCompiler {
/** Array of all compilers created so far. */
private static AtomicReference<DSCompiler[][]> compilers =
new AtomicReference<DSCompiler[][]>(null);
/** Number of free parameters. */
private final int parameters;
/** Derivation order. */
private final int order;
/** Number of partial derivatives (including the single 0 order derivative element). */
private final int[][] sizes;
/** Indirection array for partial derivatives. */
private final int[][] derivativesIndirection;
/** Indirection array of the lower derivative elements. */
private final int[] lowerIndirection;
/** Indirection arrays for multiplication. */
private final int[][][] multIndirection;
/** Indirection arrays for function composition. */
private final int[][][] compIndirection;
/** Private constructor, reserved for the factory method {@link #getCompiler(int, int)}.
* @param parameters number of free parameters
* @param order derivation order
* @param valueCompiler compiler for the value part
* @param derivativeCompiler compiler for the derivative part
* @throws NumberIsTooLargeException if order is too large
*/
private DSCompiler(final int parameters, final int order,
final DSCompiler valueCompiler, final DSCompiler derivativeCompiler)
throws NumberIsTooLargeException {
this.parameters = parameters;
this.order = order;
this.sizes = compileSizes(parameters, order, valueCompiler);
this.derivativesIndirection =
compileDerivativesIndirection(parameters, order,
valueCompiler, derivativeCompiler);
this.lowerIndirection =
compileLowerIndirection(parameters, order,
valueCompiler, derivativeCompiler);
this.multIndirection =
compileMultiplicationIndirection(parameters, order,
valueCompiler, derivativeCompiler, lowerIndirection);
this.compIndirection =
compileCompositionIndirection(parameters, order,
valueCompiler, derivativeCompiler,
sizes, derivativesIndirection);
}
/** Get the compiler for number of free parameters and order.
* @param parameters number of free parameters
* @param order derivation order
* @return cached rules set
* @throws NumberIsTooLargeException if order is too large
*/
public static DSCompiler getCompiler(int parameters, int order)
throws NumberIsTooLargeException {
// get the cached compilers
final DSCompiler[][] cache = compilers.get();
if (cache != null && cache.length > parameters &&
cache[parameters].length > order && cache[parameters][order] != null) {
// the compiler has already been created
return cache[parameters][order];
}
// we need to create more compilers
final int maxParameters = FastMath.max(parameters, cache == null ? 0 : cache.length);
final int maxOrder = FastMath.max(order, cache == null ? 0 : cache[0].length);
final DSCompiler[][] newCache = new DSCompiler[maxParameters + 1][maxOrder + 1];
if (cache != null) {
// preserve the already created compilers
for (int i = 0; i < cache.length; ++i) {
System.arraycopy(cache[i], 0, newCache[i], 0, cache[i].length);
}
}
// create the array in increasing diagonal order
for (int diag = 0; diag <= parameters + order; ++diag) {
for (int o = FastMath.max(0, diag - parameters); o <= FastMath.min(order, diag); ++o) {
final int p = diag - o;
if (newCache[p][o] == null) {
final DSCompiler valueCompiler = (p == 0) ? null : newCache[p - 1][o];
final DSCompiler derivativeCompiler = (o == 0) ? null : newCache[p][o - 1];
newCache[p][o] = new DSCompiler(p, o, valueCompiler, derivativeCompiler);
}
}
}
// atomically reset the cached compilers array
compilers.compareAndSet(cache, newCache);
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 static int[][] compileSizes(final int parameters, final int order,
final DSCompiler valueCompiler) {
final int[][] sizes = new int[parameters + 1][order + 1];
if (parameters == 0) {
Arrays.fill(sizes[0], 1);
} else {
System.arraycopy(valueCompiler.sizes, 0, sizes, 0, parameters);
sizes[parameters][0] = 1;
for (int i = 0; i < order; ++i) {
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 static int[][] compileDerivativesIndirection(final int parameters, final int order,
final DSCompiler valueCompiler,
final DSCompiler derivativeCompiler) {
if (parameters == 0 || order == 0) {
return new int[1][parameters];
}
final int vSize = valueCompiler.derivativesIndirection.length;
final int dSize = derivativeCompiler.derivativesIndirection.length;
final int[][] derivativesIndirection = new int[vSize + dSize][parameters];
// set up the indices for the value part
for (int i = 0; i < vSize; ++i) {
// copy the first indices, the last one remaining set to 0
System.arraycopy(valueCompiler.derivativesIndirection[i], 0,
derivativesIndirection[i], 0,
parameters - 1);
}
// set up the indices for the derivative part
for (int i = 0; i < dSize; ++i) {
// copy the indices
System.arraycopy(derivativeCompiler.derivativesIndirection[i], 0,
derivativesIndirection[vSize + i], 0,
parameters);
// increment the derivation order for the last parameter
derivativesIndirection[vSize + i][parameters - 1]++;
}
return derivativesIndirection;
}
/** Compile the lower derivatives indirection array.
* <p>
* This indirection array contains the indices of all elements
* except derivatives for last derivation order.
* </p>
* @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 static int[] compileLowerIndirection(final int parameters, final int order,
final DSCompiler valueCompiler,
final DSCompiler derivativeCompiler) {
if (parameters == 0 || order <= 1) {
return new int[] { 0 };
}
// this is an implementation of definition 6 in Dan Kalman's paper.
final int vSize = valueCompiler.lowerIndirection.length;
final int dSize = derivativeCompiler.lowerIndirection.length;
final int[] lowerIndirection = new int[vSize + dSize];
System.arraycopy(valueCompiler.lowerIndirection, 0, lowerIndirection, 0, vSize);
for (int i = 0; i < dSize; ++i) {
lowerIndirection[vSize + i] = valueCompiler.getSize() + derivativeCompiler.lowerIndirection[i];
}
return lowerIndirection;
}
/** Compile the multiplication indirection array.
* <p>
* This indirection array contains the indices of all pairs of elements
* involved when computing a multiplication. This allows a straightforward
* loop-based multiplication (see {@link #multiply(double[], int, double[], int, double[], int)}).
* </p>
* @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
*/
private static int[][][] compileMultiplicationIndirection(final int parameters, final int order,
final DSCompiler valueCompiler,
final DSCompiler derivativeCompiler,
final int[] lowerIndirection) {
if ((parameters == 0) || (order == 0)) {
return new int[][][] { { { 1, 0, 0 } } };
}
// this is an implementation of definition 3 in Dan Kalman's paper.
final int vSize = valueCompiler.multIndirection.length;
final int dSize = derivativeCompiler.multIndirection.length;
final int[][][] multIndirection = new int[vSize + dSize][][];
System.arraycopy(valueCompiler.multIndirection, 0, multIndirection, 0, vSize);
for (int i = 0; i < dSize; ++i) {
final int[][] dRow = derivativeCompiler.multIndirection[i];
List<int[]> row = new ArrayList<int[]>(dRow.length * 2);
for (int j = 0; j < dRow.length; ++j) {
row.add(new int[] { dRow[j][0], lowerIndirection[dRow[j][1]], vSize + dRow[j][2] });
row.add(new int[] { dRow[j][0], vSize + dRow[j][1], lowerIndirection[dRow[j][2]] });
}
// combine terms with similar derivation orders
final List<int[]> combined = new ArrayList<int[]>(row.size());
for (int j = 0; j < row.size(); ++j) {
final int[] termJ = row.get(j);
if (termJ[0] > 0) {
for (int k = j + 1; k < row.size(); ++k) {
final int[] termK = row.get(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;
}
}
combined.add(termJ);
}
}
multIndirection[vSize + i] = combined.toArray(new int[combined.size()][]);
}
return multIndirection;
}
/** Compile the function composition indirection array.
* <p>
* This indirection array contains the indices of all sets of elements
* involved when computing a composition. This allows a straightforward
* loop-based composition (see {@link #compose(double[], int, double[], double[], int)}).
* </p>
* @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
* @throws NumberIsTooLargeException if order is too large
*/
private static int[][][] compileCompositionIndirection(final int parameters, final int order,
final DSCompiler valueCompiler,
final DSCompiler derivativeCompiler,
final int[][] sizes,
final int[][] derivativesIndirection)
throws NumberIsTooLargeException {
if ((parameters == 0) || (order == 0)) {
return new int[][][] { { { 1, 0 } } };
}
final int vSize = valueCompiler.compIndirection.length;
final int dSize = derivativeCompiler.compIndirection.length;
final int[][][] compIndirection = new int[vSize + dSize][][];
// the composition rules from the value part can be reused as is
System.arraycopy(valueCompiler.compIndirection, 0, compIndirection, 0, 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
for (int i = 0; i < dSize; ++i) {
List<int[]> row = new ArrayList<int[]>();
for (int[] term : derivativeCompiler.compIndirection[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
int[] derivedTermF = new int[term.length + 1];
derivedTermF[0] = term[0]; // p
derivedTermF[1] = term[1] + 1; // f_(k+1)
int[] orders = new int[parameters];
orders[parameters - 1] = 1;
derivedTermF[term.length] = getPartialDerivativeIndex(parameters, order, sizes, orders); // g_1
for (int j = 2; j < term.length; ++j) {
// 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);
}
Arrays.sort(derivedTermF, 2, derivedTermF.length);
row.add(derivedTermF);
// derive the various g_l
for (int l = 2; l < term.length; ++l) {
int[] derivedTermG = new int[term.length];
derivedTermG[0] = term[0];
derivedTermG[1] = term[1];
for (int j = 2; j < term.length; ++j) {
// 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
System.arraycopy(derivativesIndirection[derivedTermG[j]], 0, orders, 0, parameters);
orders[parameters - 1]++;
derivedTermG[j] = getPartialDerivativeIndex(parameters, order, sizes, orders);
}
}
Arrays.sort(derivedTermG, 2, derivedTermG.length);
row.add(derivedTermG);
}
}
// combine terms with similar derivation orders
final List<int[]> combined = new ArrayList<int[]>(row.size());
for (int j = 0; j < row.size(); ++j) {
final int[] termJ = row.get(j);
if (termJ[0] > 0) {
for (int k = j + 1; k < row.size(); ++k) {
final int[] termK = row.get(k);
boolean equals = termJ.length == termK.length;
for (int l = 1; equals && l < termJ.length; ++l) {
equals &= termJ[l] == termK[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;
}
}
combined.add(termJ);
}
}
compIndirection[vSize + i] = combined.toArray(new int[combined.size()][]);
}
return compIndirection;
}
/** Get the index of a partial derivative in the array.
* <p>
* If all orders are set to 0, then the 0<sup>th</sup> order derivative
* is returned, which is the value of the function.
* </p>
* <p>The indices of derivatives are between 0 and {@link #getSize() getSize()} - 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:
* </p>
* <ul>
* <li>the index of 0<sup>th</sup> order derivative is always 0</li>
* <li>if there is only 1 {@link #getFreeParameters() free parameter}, then the
* derivatives are sorted in increasing derivation order (i.e. f at index 0, df/dp
* at index 1, d<sup>2</sup>f/dp<sup>2</sup> at index 2 ...
* d<sup>k</sup>f/dp<sup>k</sup> at index k),</li>
* <li>if the {@link #getOrder() derivation order} is 1, then the derivatives
* are sorted in increasing free parameter order (i.e. f at index 0, df/dx<sub>1</sub>
* at index 1, df/dx<sub>2</sub> at index 2 ... df/dx<sub>k</sub> at index k),</li>
* <li>all other cases are not publicly specified</li>
* </ul>
* <p>
* This method is the inverse of method {@link #getPartialDerivativeOrders(int)}
* </p>
* @param orders derivation orders with respect to each parameter
* @return index of the partial derivative
* @exception DimensionMismatchException if the numbers of parameters does not
* match the instance
* @exception NumberIsTooLargeException if sum of derivation orders is larger
* than the instance limits
* @see #getPartialDerivativeOrders(int)
*/
public int getPartialDerivativeIndex(final int ... orders)
throws DimensionMismatchException, NumberIsTooLargeException {
// safety check
if (orders.length != getFreeParameters()) {
throw new DimensionMismatchException(orders.length, getFreeParameters());
}
return getPartialDerivativeIndex(parameters, order, sizes, orders);
}
/** 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 lenght of this array must match the number of parameters)
* @return index of the partial derivative
* @exception NumberIsTooLargeException if sum of derivation orders is larger
* than the instance limits
*/
private static int getPartialDerivativeIndex(final int parameters, final int order,
final int[][] sizes, final int ... orders)
throws NumberIsTooLargeException {
// 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
int index = 0;
int m = order;
int ordersSum = 0;
for (int i = parameters - 1; i >= 0; --i) {
// derivative order for current free parameter
int derivativeOrder = orders[i];
// safety check
ordersSum += derivativeOrder;
if (ordersSum > order) {
throw new NumberIsTooLargeException(ordersSum, order, true);
}
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 <em>same</em> characteristics
* in destination derivative structure
* @throws NumberIsTooLargeException if order is too large
*/
private static int convertIndex(final int index,
final int srcP, final int[][] srcDerivativesIndirection,
final int destP, final int destO, final int[][] destSizes)
throws NumberIsTooLargeException {
int[] orders = new int[destP];
System.arraycopy(srcDerivativesIndirection[index], 0, orders, 0, FastMath.min(srcP, destP));
return getPartialDerivativeIndex(destP, destO, destSizes, orders);
}
/** Get the derivation orders for a specific index in the array.
* <p>
* This method is the inverse of {@link #getPartialDerivativeIndex(int...)}.
* </p>
* @param index of the partial derivative
* @return orders derivation orders with respect to each parameter
* @see #getPartialDerivativeIndex(int...)
*/
public int[] getPartialDerivativeOrders(final int index) {
return derivativesIndirection[index];
}
/** Get the number of free parameters.
* @return number of free parameters
*/
public int getFreeParameters() {
return parameters;
}
/** Get the derivation order.
* @return derivation order
*/
public int getOrder() {
return order;
}
/** Get the array size required for holding partial derivatives data.
* <p>
* This number includes the single 0 order derivative element, which is
* guaranteed to be stored in the first element of the array.
* </p>
* @return array size required for holding partial derivatives data
*/
public int getSize() {
return sizes[parameters][order];
}
/** Compute linear combination.
* The derivative structure built will be a1 * ds1 + a2 * ds2
* @param a1 first scale factor
* @param c1 first base (unscaled) component
* @param offset1 offset of first operand in its array
* @param a2 second scale factor
* @param c2 second base (unscaled) component
* @param offset2 offset of second operand 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
*/
public void linearCombination(final double a1, final double[] c1, final int offset1,
final double a2, final double[] c2, final int offset2,
final double[] result, final int resultOffset) {
for (int i = 0; i < getSize(); ++i) {
result[resultOffset + i] =
MathArrays.linearCombination(a1, c1[offset1 + i], a2, c2[offset2 + i]);
}
}
/** Compute linear combination.
* The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
* @param a1 first scale factor
* @param c1 first base (unscaled) component
* @param offset1 offset of first operand in its array
* @param a2 second scale factor
* @param c2 second base (unscaled) component
* @param offset2 offset of second operand in its array
* @param a3 third scale factor
* @param c3 third base (unscaled) component
* @param offset3 offset of third operand 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
*/
public void linearCombination(final double a1, final double[] c1, final int offset1,
final double a2, final double[] c2, final int offset2,
final double a3, final double[] c3, final int offset3,
final double[] result, final int resultOffset) {
for (int i = 0; i < getSize(); ++i) {
result[resultOffset + i] =
MathArrays.linearCombination(a1, c1[offset1 + i],
a2, c2[offset2 + i],
a3, c3[offset3 + i]);
}
}
/** Compute linear combination.
* The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
* @param a1 first scale factor
* @param c1 first base (unscaled) component
* @param offset1 offset of first operand in its array
* @param a2 second scale factor
* @param c2 second base (unscaled) component
* @param offset2 offset of second operand in its array
* @param a3 third scale factor
* @param c3 third base (unscaled) component
* @param offset3 offset of third operand in its array
* @param a4 fourth scale factor
* @param c4 fourth base (unscaled) component
* @param offset4 offset of fourth operand 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
*/
public void linearCombination(final double a1, final double[] c1, final int offset1,
final double a2, final double[] c2, final int offset2,
final double a3, final double[] c3, final int offset3,
final double a4, final double[] c4, final int offset4,
final double[] result, final int resultOffset) {
for (int i = 0; i < getSize(); ++i) {
result[resultOffset + i] =
MathArrays.linearCombination(a1, c1[offset1 + i],
a2, c2[offset2 + i],
a3, c3[offset3 + i],
a4, c4[offset4 + i]);
}
}
/** 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
*/
public void add(final double[] lhs, final int lhsOffset,
final double[] rhs, final int rhsOffset,
final double[] result, final int resultOffset) {
for (int i = 0; i < getSize(); ++i) {
result[resultOffset + i] = lhs[lhsOffset + i] + rhs[rhsOffset + i];
}
}
/** 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
*/
public void subtract(final double[] lhs, final int lhsOffset,
final double[] rhs, final int rhsOffset,
final double[] result, final int resultOffset) {
for (int i = 0; i < getSize(); ++i) {
result[resultOffset + i] = lhs[lhsOffset + i] - rhs[rhsOffset + i];
}
}
/** 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 <em>cannot</em> be one of
* the input arrays)
* @param resultOffset offset of the result in its array
*/
public void multiply(final double[] lhs, final int lhsOffset,
final double[] rhs, final int rhsOffset,
final double[] result, final int resultOffset) {
for (int i = 0; i < multIndirection.length; ++i) {
final int[][] mappingI = multIndirection[i];
double r = 0;
for (int j = 0; j < mappingI.length; ++j) {
r += mappingI[j][0] *
lhs[lhsOffset + mappingI[j][1]] *
rhs[rhsOffset + mappingI[j][2]];
}
result[resultOffset + i] = r;
}
}
/** 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 <em>cannot</em> be one of
* the input arrays)
* @param resultOffset offset of the result in its array
*/
public void divide(final double[] lhs, final int lhsOffset,
final double[] rhs, final int rhsOffset,
final double[] result, final int resultOffset) {
final double[] reciprocal = new double[getSize()];
pow(rhs, lhsOffset, -1, reciprocal, 0);
multiply(lhs, lhsOffset, reciprocal, 0, result, resultOffset);
}
/** Perform remainder of two derivative structures.
* @param lhs array holding left hand side of remainder
* @param lhsOffset offset of the left hand side in its array
* @param rhs array right hand side of remainder
* @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
*/
public void remainder(final double[] lhs, final int lhsOffset,
final double[] rhs, final int rhsOffset,
final double[] result, final int resultOffset) {
// compute k such that lhs % rhs = lhs - k rhs
final double rem = FastMath.IEEEremainder(lhs[lhsOffset], rhs[rhsOffset]);
final double k = FastMath.rint((lhs[lhsOffset] - rem) / rhs[rhsOffset]);
// set up value
result[resultOffset] = rem;
// set up partial derivatives
for (int i = 1; i < getSize(); ++i) {
result[resultOffset + i] = lhs[lhsOffset + i] - k * rhs[rhsOffset + i];
}
}
/** Compute power of a double to a derivative structure.
* @param a number to exponentiate
* @param operand array holding the power
* @param operandOffset offset of the power in its array
* @param result array where result must be stored (for
* power the result array <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
* @since 3.3
*/
public void pow(final double a,
final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
// [a^x, ln(a) a^x, ln(a)^2 a^x,, ln(a)^3 a^x, ... ]
final double[] function = new double[1 + order];
if (a == 0) {
if (operand[operandOffset] == 0) {
function[0] = 1;
double infinity = Double.POSITIVE_INFINITY;
for (int i = 1; i < function.length; ++i) {
infinity = -infinity;
function[i] = infinity;
}
} else if (operand[operandOffset] < 0) {
Arrays.fill(function, Double.NaN);
}
} else {
function[0] = FastMath.pow(a, operand[operandOffset]);
final double lnA = FastMath.log(a);
for (int i = 1; i < function.length; ++i) {
function[i] = lnA * function[i - 1];
}
}
// 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void pow(final double[] operand, final int operandOffset, final double p,
final double[] result, final int resultOffset) {
// create the function value and derivatives
// [x^p, px^(p-1), p(p-1)x^(p-2), ... ]
double[] function = new double[1 + order];
double xk = FastMath.pow(operand[operandOffset], p - order);
for (int i = order; i > 0; --i) {
function[i] = xk;
xk *= operand[operandOffset];
}
function[0] = xk;
double coefficient = p;
for (int i = 1; i <= order; ++i) {
function[i] *= coefficient;
coefficient *= p - i;
}
// 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void pow(final double[] operand, final int operandOffset, final int n,
final double[] result, final int resultOffset) {
if (n == 0) {
// special case, x^0 = 1 for all x
result[resultOffset] = 1.0;
Arrays.fill(result, resultOffset + 1, resultOffset + getSize(), 0);
return;
}
// create the power function value and derivatives
// [x^n, nx^(n-1), n(n-1)x^(n-2), ... ]
double[] function = new double[1 + order];
if (n > 0) {
// strictly positive power
final int maxOrder = FastMath.min(order, n);
double xk = FastMath.pow(operand[operandOffset], n - maxOrder);
for (int i = maxOrder; i > 0; --i) {
function[i] = xk;
xk *= operand[operandOffset];
}
function[0] = xk;
} else {
// strictly negative power
final double inv = 1.0 / operand[operandOffset];
double xk = FastMath.pow(inv, -n);
for (int i = 0; i <= order; ++i) {
function[i] = xk;
xk *= inv;
}
}
double coefficient = n;
for (int i = 1; i <= order; ++i) {
function[i] *= coefficient;
coefficient *= n - i;
}
// 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void pow(final double[] x, final int xOffset,
final double[] y, final int yOffset,
final double[] result, final int resultOffset) {
final double[] logX = new double[getSize()];
log(x, xOffset, logX, 0);
final double[] yLogX = new double[getSize()];
multiply(logX, 0, y, yOffset, yLogX, 0);
exp(yLogX, 0, result, resultOffset);
}
/** Compute n<sup>th</sup> root of a derivative structure.
* @param operand array holding the operand
* @param operandOffset offset of the operand in its array
* @param n order of the root
* @param result array where result must be stored (for
* n<sup>th</sup> root the result array <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void rootN(final double[] operand, final int operandOffset, final int n,
final double[] result, final int resultOffset) {
// create the function value and derivatives
// [x^(1/n), (1/n)x^((1/n)-1), (1-n)/n^2x^((1/n)-2), ... ]
double[] function = new double[1 + order];
double xk;
if (n == 2) {
function[0] = FastMath.sqrt(operand[operandOffset]);
xk = 0.5 / function[0];
} else if (n == 3) {
function[0] = FastMath.cbrt(operand[operandOffset]);
xk = 1.0 / (3.0 * function[0] * function[0]);
} else {
function[0] = FastMath.pow(operand[operandOffset], 1.0 / n);
xk = 1.0 / (n * FastMath.pow(function[0], n - 1));
}
final double nReciprocal = 1.0 / n;
final double xReciprocal = 1.0 / operand[operandOffset];
for (int i = 1; i <= order; ++i) {
function[i] = xk;
xk *= xReciprocal * (nReciprocal - i);
}
// 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void exp(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
Arrays.fill(function, FastMath.exp(operand[operandOffset]));
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** Compute exp(x) - 1 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void expm1(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
function[0] = FastMath.expm1(operand[operandOffset]);
Arrays.fill(function, 1, 1 + order, FastMath.exp(operand[operandOffset]));
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void log(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
function[0] = FastMath.log(operand[operandOffset]);
if (order > 0) {
double inv = 1.0 / operand[operandOffset];
double xk = inv;
for (int i = 1; i <= order; ++i) {
function[i] = xk;
xk *= -i * inv;
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** Computes shifted 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
* shifted logarithm the result array <em>cannot</em> be the input array)
* @param resultOffset offset of the result in its array
*/
public void log1p(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
function[0] = FastMath.log1p(operand[operandOffset]);
if (order > 0) {
double inv = 1.0 / (1.0 + operand[operandOffset]);
double xk = inv;
for (int i = 1; i <= order; ++i) {
function[i] = xk;
xk *= -i * inv;
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** Computes base 10 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
* base 10 logarithm the result array <em>cannot</em> be the input array)
* @param resultOffset offset of the result in its array
*/
public void log10(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
function[0] = FastMath.log10(operand[operandOffset]);
if (order > 0) {
double inv = 1.0 / operand[operandOffset];
double xk = inv / FastMath.log(10.0);
for (int i = 1; i <= order; ++i) {
function[i] = xk;
xk *= -i * inv;
}
}
// 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void cos(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
function[0] = FastMath.cos(operand[operandOffset]);
if (order > 0) {
function[1] = -FastMath.sin(operand[operandOffset]);
for (int i = 2; i <= order; ++i) {
function[i] = -function[i - 2];
}
}
// apply function composition
compose(operand, operandOffset, function, 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void sin(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
function[0] = FastMath.sin(operand[operandOffset]);
if (order > 0) {
function[1] = FastMath.cos(operand[operandOffset]);
for (int i = 2; i <= order; ++i) {
function[i] = -function[i - 2];
}
}
// 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void tan(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
final double[] function = new double[1 + order];
final double t = FastMath.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
final double[] p = new double[order + 2];
p[1] = 1;
final double t2 = t * t;
for (int n = 1; n <= order; ++n) {
// update and evaluate polynomial P_n(t)
double v = 0;
p[n + 1] = n * p[n];
for (int k = n + 1; k >= 0; k -= 2) {
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];
}
}
if ((n & 0x1) == 0) {
v *= t;
}
function[n] = v;
}
}
// 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void acos(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
final double x = operand[operandOffset];
function[0] = FastMath.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
final double[] p = new double[order];
p[0] = -1;
final double x2 = x * x;
final double f = 1.0 / (1 - x2);
double coeff = FastMath.sqrt(f);
function[1] = coeff * p[0];
for (int n = 2; n <= order; ++n) {
// update and evaluate polynomial P_n(x)
double v = 0;
p[n - 1] = (n - 1) * p[n - 2];
for (int k = n - 1; k >= 0; k -= 2) {
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];
}
}
if ((n & 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void asin(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
final double x = operand[operandOffset];
function[0] = FastMath.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
final double[] p = new double[order];
p[0] = 1;
final double x2 = x * x;
final double f = 1.0 / (1 - x2);
double coeff = FastMath.sqrt(f);
function[1] = coeff * p[0];
for (int n = 2; n <= order; ++n) {
// update and evaluate polynomial P_n(x)
double v = 0;
p[n - 1] = (n - 1) * p[n - 2];
for (int k = n - 1; k >= 0; k -= 2) {
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];
}
}
if ((n & 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void atan(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
final double x = operand[operandOffset];
function[0] = FastMath.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
final double[] q = new double[order];
q[0] = 1;
final double x2 = x * x;
final double f = 1.0 / (1 + x2);
double coeff = f;
function[1] = coeff * q[0];
for (int n = 2; n <= order; ++n) {
// update and evaluate polynomial Q_n(x)
double v = 0;
q[n - 1] = -n * q[n - 2];
for (int k = n - 1; k >= 0; k -= 2) {
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];
}
}
if ((n & 0x1) == 0) {
v *= x;
}
coeff *= f;
function[n] = coeff * v;
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** Compute two arguments arc tangent of a derivative structure.
* @param y array holding the first operand
* @param yOffset offset of the first operand in its array
* @param x array holding the second operand
* @param xOffset offset of the second operand in its array
* @param result array where result must be stored (for
* two arguments arc tangent the result array <em>cannot</em>
* be the input array)
* @param resultOffset offset of the result in its array
*/
public void atan2(final double[] y, final int yOffset,
final double[] x, final int xOffset,
final double[] result, final int resultOffset) {
// compute r = sqrt(x^2+y^2)
double[] tmp1 = new double[getSize()];
multiply(x, xOffset, x, xOffset, tmp1, 0); // x^2
double[] tmp2 = new double[getSize()];
multiply(y, yOffset, y, yOffset, tmp2, 0); // y^2
add(tmp1, 0, tmp2, 0, tmp2, 0); // x^2 + y^2
rootN(tmp2, 0, 2, tmp1, 0); // r = sqrt(x^2 + y^2)
if (x[xOffset] >= 0) {
// compute atan2(y, x) = 2 atan(y / (r + x))
add(tmp1, 0, x, xOffset, tmp2, 0); // r + x
divide(y, yOffset, tmp2, 0, tmp1, 0); // y /(r + x)
atan(tmp1, 0, tmp2, 0); // atan(y / (r + x))
for (int i = 0; i < tmp2.length; ++i) {
result[resultOffset + i] = 2 * tmp2[i]; // 2 * atan(y / (r + x))
}
} else {
// compute atan2(y, x) = +/- pi - 2 atan(y / (r - x))
subtract(tmp1, 0, x, xOffset, tmp2, 0); // r - x
divide(y, yOffset, tmp2, 0, tmp1, 0); // y /(r - x)
atan(tmp1, 0, tmp2, 0); // atan(y / (r - x))
result[resultOffset] =
((tmp2[0] <= 0) ? -FastMath.PI : FastMath.PI) - 2 * tmp2[0]; // +/-pi - 2 * atan(y / (r - x))
for (int i = 1; i < tmp2.length; ++i) {
result[resultOffset + i] = -2 * tmp2[i]; // +/-pi - 2 * atan(y / (r - x))
}
}
// fix value to take special cases (+0/+0, +0/-0, -0/+0, -0/-0, +/-infinity) correctly
result[resultOffset] = FastMath.atan2(y[yOffset], x[xOffset]);
}
/** 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void cosh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
function[0] = FastMath.cosh(operand[operandOffset]);
if (order > 0) {
function[1] = FastMath.sinh(operand[operandOffset]);
for (int i = 2; i <= order; ++i) {
function[i] = function[i - 2];
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void sinh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
function[0] = FastMath.sinh(operand[operandOffset]);
if (order > 0) {
function[1] = FastMath.cosh(operand[operandOffset]);
for (int i = 2; i <= order; ++i) {
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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void tanh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
final double[] function = new double[1 + order];
final double t = FastMath.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
final double[] p = new double[order + 2];
p[1] = 1;
final double t2 = t * t;
for (int n = 1; n <= order; ++n) {
// update and evaluate polynomial P_n(t)
double v = 0;
p[n + 1] = -n * p[n];
for (int k = n + 1; k >= 0; k -= 2) {
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];
}
}
if ((n & 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void acosh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
final double x = operand[operandOffset];
function[0] = FastMath.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
final double[] p = new double[order];
p[0] = 1;
final double x2 = x * x;
final double f = 1.0 / (x2 - 1);
double coeff = FastMath.sqrt(f);
function[1] = coeff * p[0];
for (int n = 2; n <= order; ++n) {
// update and evaluate polynomial P_n(x)
double v = 0;
p[n - 1] = (1 - n) * p[n - 2];
for (int k = n - 1; k >= 0; k -= 2) {
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];
}
}
if ((n & 0x1) == 0) {
v *= x;
}
coeff *= f;
function[n] = coeff * v;
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void asinh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
final double x = operand[operandOffset];
function[0] = FastMath.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
final double[] p = new double[order];
p[0] = 1;
final double x2 = x * x;
final double f = 1.0 / (1 + x2);
double coeff = FastMath.sqrt(f);
function[1] = coeff * p[0];
for (int n = 2; n <= order; ++n) {
// update and evaluate polynomial P_n(x)
double v = 0;
p[n - 1] = (1 - n) * p[n - 2];
for (int k = n - 1; k >= 0; k -= 2) {
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];
}
}
if ((n & 0x1) == 0) {
v *= x;
}
coeff *= f;
function[n] = coeff * v;
}
}
// apply function composition
compose(operand, operandOffset, function, result, resultOffset);
}
/** 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 <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void atanh(final double[] operand, final int operandOffset,
final double[] result, final int resultOffset) {
// create the function value and derivatives
double[] function = new double[1 + order];
final double x = operand[operandOffset];
function[0] = FastMath.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
final double[] q = new double[order];
q[0] = 1;
final double x2 = x * x;
final double f = 1.0 / (1 - x2);
double coeff = f;
function[1] = coeff * q[0];
for (int n = 2; n <= order; ++n) {
// update and evaluate polynomial Q_n(x)
double v = 0;
q[n - 1] = n * q[n - 2];
for (int k = n - 1; k >= 0; k -= 2) {
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];
}
}
if ((n & 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 {@code operand[operandOffset]}).
* @param result array where result must be stored (for
* composition the result array <em>cannot</em> be the input
* array)
* @param resultOffset offset of the result in its array
*/
public void compose(final double[] operand, final int operandOffset, final double[] f,
final double[] result, final int resultOffset) {
for (int i = 0; i < compIndirection.length; ++i) {
final int[][] mappingI = compIndirection[i];
double r = 0;
for (int j = 0; j < mappingI.length; ++j) {
final int[] mappingIJ = mappingI[j];
double product = mappingIJ[0] * f[mappingIJ[1]];
for (int k = 2; k < mappingIJ.length; ++k) {
product *= operand[operandOffset + mappingIJ[k]];
}
r += product;
}
result[resultOffset + i] = r;
}
}
/** Evaluate Taylor expansion of a derivative structure.
* @param ds array holding the derivative structure
* @param dsOffset offset of the derivative structure in its array
* @param delta parameters offsets (Δx, Δy, ...)
* @return value of the Taylor expansion at x + Δx, y + Δy, ...
* @throws MathArithmeticException if factorials becomes too large
*/
public double taylor(final double[] ds, final int dsOffset, final double ... delta)
throws MathArithmeticException {
double value = 0;
for (int i = getSize() - 1; i >= 0; --i) {
final int[] orders = getPartialDerivativeOrders(i);
double term = ds[dsOffset + i];
for (int k = 0; k < orders.length; ++k) {
if (orders[k] > 0) {
try {
term *= FastMath.pow(delta[k], orders[k]) /
CombinatoricsUtils.factorial(orders[k]);
} catch (NotPositiveException e) {
// this cannot happen
throw new MathInternalError(e);
}
}
}
value += term;
}
return value;
}
/** Check rules set compatibility.
* @param compiler other compiler to check against instance
* @exception DimensionMismatchException if number of free parameters or orders are inconsistent
*/
public void checkCompatibility(final DSCompiler compiler)
throws DimensionMismatchException {
if (parameters != compiler.parameters) {
throw new DimensionMismatchException(parameters, compiler.parameters);
}
if (order != compiler.order) {
throw new DimensionMismatchException(order, compiler.order);
}
}
}
2364
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
