001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 package org.apache.commons.math3.analysis.differentiation; 018 019 import java.util.ArrayList; 020 import java.util.Arrays; 021 import java.util.List; 022 import java.util.concurrent.atomic.AtomicReference; 023 024 import org.apache.commons.math3.exception.DimensionMismatchException; 025 import org.apache.commons.math3.exception.NumberIsTooLargeException; 026 import org.apache.commons.math3.util.ArithmeticUtils; 027 import org.apache.commons.math3.util.FastMath; 028 import org.apache.commons.math3.util.MathArrays; 029 030 /** Class holding "compiled" computation rules for derivative structures. 031 * <p>This class implements the computation rules described in Dan Kalman's paper <a 032 * href="http://www.math.american.edu/People/kalman/pdffiles/mmgautodiff.pdf">Doubly 033 * Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75, 034 * no. 3, June 2002. However, in order to avoid performances bottlenecks, the recursive 035 * rules are "compiled" once in an unfold form. This class does this recursion unrolling 036 * and stores the computation rules as simple loops with pre-computed indirection arrays.</p> 037 * <p> 038 * This class maps all derivative computation into single dimension arrays that hold the 039 * value and partial derivatives. The class does not hold these arrays, which remains under 040 * the responsibility of the caller. For each combination of number of free parameters and 041 * derivation order, only one compiler is necessary, and this compiler will be used to 042 * perform computations on all arrays provided to it, which can represent hundreds or 043 * thousands of different parameters kept together with all theur partial derivatives. 044 * </p> 045 * <p> 046 * The arrays on which compilers operate contain only the partial derivatives together 047 * with the 0<sup>th</sup> derivative, i.e. the value. The partial derivatives are stored in 048 * a compiler-specific order, which can be retrieved using methods {@link 049 * #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} and {@link 050 * #getPartialDerivativeOrders(int)}. The value is guaranteed to be stored as the first element 051 * (i.e. the {@link #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} method returns 052 * 0 when called with 0 for all derivation orders and {@link #getPartialDerivativeOrders(int) 053 * getPartialDerivativeOrders} returns an array filled with 0 when called with 0 as the index). 054 * </p> 055 * <p> 056 * Note that the ordering changes with number of parameters and derivation order. For example 057 * given 2 parameters x and y, df/dy is stored at index 2 when derivation order is set to 1 (in 058 * this case the array has three elements: f, df/dx and df/dy). If derivation order is set to 059 * 2, then df/dy will be stored at index 3 (in this case the array has six elements: f, df/dx, 060 * df/dxdx, df/dy, df/dxdy and df/dydy). 061 * </p> 062 * <p> 063 * Given this structure, users can perform some simple operations like adding, subtracting 064 * or multiplying constants and negating the elements by themselves, knowing if they want to 065 * mutate their array or create a new array. These simple operations are not provided by 066 * the compiler. The compiler provides only the more complex operations between several arrays. 067 * </p> 068 * <p>This class is mainly used as the engine for scalar variable {@link DerivativeStructure}. 069 * It can also be used directly to hold several variables in arrays for more complex data 070 * structures. User can for example store a vector of n variables depending on three x, y 071 * and z free parameters in one array as follows: 072 * <pre> 073 * // parameter 0 is x, parameter 1 is y, parameter 2 is z 074 * int parameters = 3; 075 * DSCompiler compiler = DSCompiler.getCompiler(parameters, order); 076 * int size = compiler.getSize(); 077 * 078 * // pack all elements in a single array 079 * double[] array = new double[n * size]; 080 * for (int i = 0; i < n; ++i) { 081 * 082 * // we know value is guaranteed to be the first element 083 * array[i * size] = v[i]; 084 * 085 * // we don't know where first derivatives are stored, so we ask the compiler 086 * array[i * size + compiler.getPartialDerivativeIndex(1, 0, 0) = dvOnDx[i][0]; 087 * array[i * size + compiler.getPartialDerivativeIndex(0, 1, 0) = dvOnDy[i][0]; 088 * array[i * size + compiler.getPartialDerivativeIndex(0, 0, 1) = dvOnDz[i][0]; 089 * 090 * // we let all higher order derivatives set to 0 091 * 092 * } 093 * </pre> 094 * Then in another function, user can perform some operations on all elements stored 095 * in the single array, such as a simple product of all variables: 096 * <pre> 097 * // compute the product of all elements 098 * double[] product = new double[size]; 099 * prod[0] = 1.0; 100 * for (int i = 0; i < n; ++i) { 101 * double[] tmp = product.clone(); 102 * compiler.multiply(tmp, 0, array, i * size, product, 0); 103 * } 104 * 105 * // value 106 * double p = product[0]; 107 * 108 * // first derivatives 109 * double dPdX = product[compiler.getPartialDerivativeIndex(1, 0, 0)]; 110 * double dPdY = product[compiler.getPartialDerivativeIndex(0, 1, 0)]; 111 * double dPdZ = product[compiler.getPartialDerivativeIndex(0, 0, 1)]; 112 * 113 * // cross derivatives (assuming order was at least 2) 114 * double dPdXdX = product[compiler.getPartialDerivativeIndex(2, 0, 0)]; 115 * double dPdXdY = product[compiler.getPartialDerivativeIndex(1, 1, 0)]; 116 * double dPdXdZ = product[compiler.getPartialDerivativeIndex(1, 0, 1)]; 117 * double dPdYdY = product[compiler.getPartialDerivativeIndex(0, 2, 0)]; 118 * double dPdYdZ = product[compiler.getPartialDerivativeIndex(0, 1, 1)]; 119 * double dPdZdZ = product[compiler.getPartialDerivativeIndex(0, 0, 2)]; 120 * </p> 121 * @see DerivativeStructure 122 * @version $Id: DSCompiler.java 1421949 2012-12-14 15:53:32Z luc $ 123 * @since 3.1 124 */ 125 public class DSCompiler { 126 127 /** Array of all compilers created so far. */ 128 private static AtomicReference<DSCompiler[][]> compilers = 129 new AtomicReference<DSCompiler[][]>(null); 130 131 /** Number of free parameters. */ 132 private final int parameters; 133 134 /** Derivation order. */ 135 private final int order; 136 137 /** Number of partial derivatives (including the single 0 order derivative element). */ 138 private final int[][] sizes; 139 140 /** Indirection array for partial derivatives. */ 141 private final int[][] derivativesIndirection; 142 143 /** Indirection array of the lower derivative elements. */ 144 private final int[] lowerIndirection; 145 146 /** Indirection arrays for multiplication. */ 147 private final int[][][] multIndirection; 148 149 /** Indirection arrays for function composition. */ 150 private final int[][][] compIndirection; 151 152 /** Private constructor, reserved for the factory method {@link #getCompiler(int, int)}. 153 * @param parameters number of free parameters 154 * @param order derivation order 155 * @param valueCompiler compiler for the value part 156 * @param derivativeCompiler compiler for the derivative part 157 */ 158 private DSCompiler(final int parameters, final int order, 159 final DSCompiler valueCompiler, final DSCompiler derivativeCompiler) { 160 161 this.parameters = parameters; 162 this.order = order; 163 this.sizes = compileSizes(parameters, order, valueCompiler); 164 this.derivativesIndirection = 165 compileDerivativesIndirection(parameters, order, 166 valueCompiler, derivativeCompiler); 167 this.lowerIndirection = 168 compileLowerIndirection(parameters, order, 169 valueCompiler, derivativeCompiler); 170 this.multIndirection = 171 compileMultiplicationIndirection(parameters, order, 172 valueCompiler, derivativeCompiler, lowerIndirection); 173 this.compIndirection = 174 compileCompositionIndirection(parameters, order, 175 valueCompiler, derivativeCompiler, 176 sizes, derivativesIndirection); 177 178 } 179 180 /** Get the compiler for number of free parameters and order. 181 * @param parameters number of free parameters 182 * @param order derivation order 183 * @return cached rules set 184 */ 185 public static DSCompiler getCompiler(int parameters, int order) { 186 187 // get the cached compilers 188 final DSCompiler[][] cache = compilers.get(); 189 if (cache != null && cache.length > parameters && cache[parameters].length > order) { 190 if (cache[parameters][order] != null) { 191 // the compiler has already been created 192 return cache[parameters][order]; 193 } 194 } 195 196 // we need to create more compilers 197 final int maxParameters = FastMath.max(parameters, cache == null ? 0 : cache.length); 198 final int maxOrder = FastMath.max(order, cache == null ? 0 : cache[0].length); 199 final DSCompiler[][] newCache = new DSCompiler[maxParameters + 1][maxOrder + 1]; 200 201 if (cache != null) { 202 // preserve the already created compilers 203 for (int i = 0; i < cache.length; ++i) { 204 System.arraycopy(cache[i], 0, newCache[i], 0, cache[i].length); 205 } 206 } 207 208 // create the array in increasing diagonal order 209 for (int diag = 0; diag <= parameters + order; ++diag) { 210 for (int o = FastMath.max(0, diag - parameters); o <= FastMath.min(order, diag); ++o) { 211 final int p = diag - o; 212 if (newCache[p][o] == null) { 213 final DSCompiler valueCompiler = (p == 0) ? null : newCache[p - 1][o]; 214 final DSCompiler derivativeCompiler = (o == 0) ? null : newCache[p][o - 1]; 215 newCache[p][o] = new DSCompiler(p, o, valueCompiler, derivativeCompiler); 216 } 217 } 218 } 219 220 // atomically reset the cached compilers array 221 compilers.compareAndSet(cache, newCache); 222 223 return newCache[parameters][order]; 224 225 } 226 227 /** Compile the sizes array. 228 * @param parameters number of free parameters 229 * @param order derivation order 230 * @param valueCompiler compiler for the value part 231 * @return sizes array 232 */ 233 private static int[][] compileSizes(final int parameters, final int order, 234 final DSCompiler valueCompiler) { 235 236 final int[][] sizes = new int[parameters + 1][order + 1]; 237 if (parameters == 0) { 238 Arrays.fill(sizes[0], 1); 239 } else { 240 System.arraycopy(valueCompiler.sizes, 0, sizes, 0, parameters); 241 sizes[parameters][0] = 1; 242 for (int i = 0; i < order; ++i) { 243 sizes[parameters][i + 1] = sizes[parameters][i] + sizes[parameters - 1][i + 1]; 244 } 245 } 246 247 return sizes; 248 249 } 250 251 /** Compile the derivatives indirection array. 252 * @param parameters number of free parameters 253 * @param order derivation order 254 * @param valueCompiler compiler for the value part 255 * @param derivativeCompiler compiler for the derivative part 256 * @return derivatives indirection array 257 */ 258 private static int[][] compileDerivativesIndirection(final int parameters, final int order, 259 final DSCompiler valueCompiler, 260 final DSCompiler derivativeCompiler) { 261 262 if (parameters == 0 || order == 0) { 263 return new int[1][parameters]; 264 } 265 266 final int vSize = valueCompiler.derivativesIndirection.length; 267 final int dSize = derivativeCompiler.derivativesIndirection.length; 268 final int[][] derivativesIndirection = new int[vSize + dSize][parameters]; 269 270 // set up the indices for the value part 271 for (int i = 0; i < vSize; ++i) { 272 // copy the first indices, the last one remaining set to 0 273 System.arraycopy(valueCompiler.derivativesIndirection[i], 0, 274 derivativesIndirection[i], 0, 275 parameters - 1); 276 } 277 278 // set up the indices for the derivative part 279 for (int i = 0; i < dSize; ++i) { 280 281 // copy the indices 282 System.arraycopy(derivativeCompiler.derivativesIndirection[i], 0, 283 derivativesIndirection[vSize + i], 0, 284 parameters); 285 286 // increment the derivation order for the last parameter 287 derivativesIndirection[vSize + i][parameters - 1]++; 288 289 } 290 291 return derivativesIndirection; 292 293 } 294 295 /** Compile the lower derivatives indirection array. 296 * <p> 297 * This indirection array contains the indices of all elements 298 * except derivatives for last derivation order. 299 * </p> 300 * @param parameters number of free parameters 301 * @param order derivation order 302 * @param valueCompiler compiler for the value part 303 * @param derivativeCompiler compiler for the derivative part 304 * @return lower derivatives indirection array 305 */ 306 private static int[] compileLowerIndirection(final int parameters, final int order, 307 final DSCompiler valueCompiler, 308 final DSCompiler derivativeCompiler) { 309 310 if (parameters == 0 || order <= 1) { 311 return new int[] { 0 }; 312 } 313 314 // this is an implementation of definition 6 in Dan Kalman's paper. 315 final int vSize = valueCompiler.lowerIndirection.length; 316 final int dSize = derivativeCompiler.lowerIndirection.length; 317 final int[] lowerIndirection = new int[vSize + dSize]; 318 System.arraycopy(valueCompiler.lowerIndirection, 0, lowerIndirection, 0, vSize); 319 for (int i = 0; i < dSize; ++i) { 320 lowerIndirection[vSize + i] = valueCompiler.getSize() + derivativeCompiler.lowerIndirection[i]; 321 } 322 323 return lowerIndirection; 324 325 } 326 327 /** Compile the multiplication indirection array. 328 * <p> 329 * This indirection array contains the indices of all pairs of elements 330 * involved when computing a multiplication. This allows a straightforward 331 * loop-based multiplication (see {@link #multiply(double[], int, double[], int, double[], int)}). 332 * </p> 333 * @param parameters number of free parameters 334 * @param order derivation order 335 * @param valueCompiler compiler for the value part 336 * @param derivativeCompiler compiler for the derivative part 337 * @param lowerIndirection lower derivatives indirection array 338 * @return multiplication indirection array 339 */ 340 private static int[][][] compileMultiplicationIndirection(final int parameters, final int order, 341 final DSCompiler valueCompiler, 342 final DSCompiler derivativeCompiler, 343 final int[] lowerIndirection) { 344 345 if ((parameters == 0) || (order == 0)) { 346 return new int[][][] { { { 1, 0, 0 } } }; 347 } 348 349 // this is an implementation of definition 3 in Dan Kalman's paper. 350 final int vSize = valueCompiler.multIndirection.length; 351 final int dSize = derivativeCompiler.multIndirection.length; 352 final int[][][] multIndirection = new int[vSize + dSize][][]; 353 354 System.arraycopy(valueCompiler.multIndirection, 0, multIndirection, 0, vSize); 355 356 for (int i = 0; i < dSize; ++i) { 357 final int[][] dRow = derivativeCompiler.multIndirection[i]; 358 List<int[]> row = new ArrayList<int[]>(); 359 for (int j = 0; j < dRow.length; ++j) { 360 row.add(new int[] { dRow[j][0], lowerIndirection[dRow[j][1]], vSize + dRow[j][2] }); 361 row.add(new int[] { dRow[j][0], vSize + dRow[j][1], lowerIndirection[dRow[j][2]] }); 362 } 363 364 // combine terms with similar derivation orders 365 final List<int[]> combined = new ArrayList<int[]>(row.size()); 366 for (int j = 0; j < row.size(); ++j) { 367 final int[] termJ = row.get(j); 368 if (termJ[0] > 0) { 369 for (int k = j + 1; k < row.size(); ++k) { 370 final int[] termK = row.get(k); 371 if (termJ[1] == termK[1] && termJ[2] == termK[2]) { 372 // combine termJ and termK 373 termJ[0] += termK[0]; 374 // make sure we will skip termK later on in the outer loop 375 termK[0] = 0; 376 } 377 } 378 combined.add(termJ); 379 } 380 } 381 382 multIndirection[vSize + i] = combined.toArray(new int[combined.size()][]); 383 384 } 385 386 return multIndirection; 387 388 } 389 390 /** Compile the function composition indirection array. 391 * <p> 392 * This indirection array contains the indices of all sets of elements 393 * involved when computing a composition. This allows a straightforward 394 * loop-based composition (see {@link #compose(double[], int, double[], double[], int)}). 395 * </p> 396 * @param parameters number of free parameters 397 * @param order derivation order 398 * @param valueCompiler compiler for the value part 399 * @param derivativeCompiler compiler for the derivative part 400 * @param sizes sizes array 401 * @param derivativesIndirection derivatives indirection array 402 * @return multiplication indirection array 403 */ 404 private static int[][][] compileCompositionIndirection(final int parameters, final int order, 405 final DSCompiler valueCompiler, 406 final DSCompiler derivativeCompiler, 407 final int[][] sizes, 408 final int[][] derivativesIndirection) { 409 410 if ((parameters == 0) || (order == 0)) { 411 return new int[][][] { { { 1, 0 } } }; 412 } 413 414 final int vSize = valueCompiler.compIndirection.length; 415 final int dSize = derivativeCompiler.compIndirection.length; 416 final int[][][] compIndirection = new int[vSize + dSize][][]; 417 418 // the composition rules from the value part can be reused as is 419 System.arraycopy(valueCompiler.compIndirection, 0, compIndirection, 0, vSize); 420 421 // the composition rules for the derivative part are deduced by 422 // differentiation the rules from the underlying compiler once 423 // with respect to the parameter this compiler handles and the 424 // underlying one did not handle 425 for (int i = 0; i < dSize; ++i) { 426 List<int[]> row = new ArrayList<int[]>(); 427 for (int[] term : derivativeCompiler.compIndirection[i]) { 428 429 // handle term p * f_k(g(x)) * g_l1(x) * g_l2(x) * ... * g_lp(x) 430 431 // derive the first factor in the term: f_k with respect to new parameter 432 int[] derivedTermF = new int[term.length + 1]; 433 derivedTermF[0] = term[0]; // p 434 derivedTermF[1] = term[1] + 1; // f_(k+1) 435 int[] orders = new int[parameters]; 436 orders[parameters - 1] = 1; 437 derivedTermF[term.length] = getPartialDerivativeIndex(parameters, order, sizes, orders); // g_1 438 for (int j = 2; j < term.length; ++j) { 439 // convert the indices as the mapping for the current order 440 // is different from the mapping with one less order 441 derivedTermF[j] = convertIndex(term[j], parameters, 442 derivativeCompiler.derivativesIndirection, 443 parameters, order, sizes); 444 } 445 Arrays.sort(derivedTermF, 2, derivedTermF.length); 446 row.add(derivedTermF); 447 448 // derive the various g_l 449 for (int l = 2; l < term.length; ++l) { 450 int[] derivedTermG = new int[term.length]; 451 derivedTermG[0] = term[0]; 452 derivedTermG[1] = term[1]; 453 for (int j = 2; j < term.length; ++j) { 454 // convert the indices as the mapping for the current order 455 // is different from the mapping with one less order 456 derivedTermG[j] = convertIndex(term[j], parameters, 457 derivativeCompiler.derivativesIndirection, 458 parameters, order, sizes); 459 if (j == l) { 460 // derive this term 461 System.arraycopy(derivativesIndirection[derivedTermG[j]], 0, orders, 0, parameters); 462 orders[parameters - 1]++; 463 derivedTermG[j] = getPartialDerivativeIndex(parameters, order, sizes, orders); 464 } 465 } 466 Arrays.sort(derivedTermG, 2, derivedTermG.length); 467 row.add(derivedTermG); 468 } 469 470 } 471 472 // combine terms with similar derivation orders 473 final List<int[]> combined = new ArrayList<int[]>(row.size()); 474 for (int j = 0; j < row.size(); ++j) { 475 final int[] termJ = row.get(j); 476 if (termJ[0] > 0) { 477 for (int k = j + 1; k < row.size(); ++k) { 478 final int[] termK = row.get(k); 479 boolean equals = termJ.length == termK.length; 480 for (int l = 1; equals && l < termJ.length; ++l) { 481 equals &= termJ[l] == termK[l]; 482 } 483 if (equals) { 484 // combine termJ and termK 485 termJ[0] += termK[0]; 486 // make sure we will skip termK later on in the outer loop 487 termK[0] = 0; 488 } 489 } 490 combined.add(termJ); 491 } 492 } 493 494 compIndirection[vSize + i] = combined.toArray(new int[combined.size()][]); 495 496 } 497 498 return compIndirection; 499 500 } 501 502 /** Get the index of a partial derivative in the array. 503 * <p> 504 * If all orders are set to 0, then the 0<sup>th</sup> order derivative 505 * is returned, which is the value of the function. 506 * </p> 507 * <p>The indices of derivatives are between 0 and {@link #getSize() getSize()} - 1. 508 * Their specific order is fixed for a given compiler, but otherwise not 509 * publicly specified. There are however some simple cases which have guaranteed 510 * indices: 511 * </p> 512 * <ul> 513 * <li>the index of 0<sup>th</sup> order derivative is always 0</li> 514 * <li>if there is only 1 {@link #getFreeParameters() free parameter}, then the 515 * derivatives are sorted in increasing derivation order (i.e. f at index 0, df/dp 516 * at index 1, d<sup>2</sup>f/dp<sup>2</sup> at index 2 ... 517 * d<sup>k</sup>f/dp<sup>k</sup> at index k),</li> 518 * <li>if the {@link #getOrder() derivation order} is 1, then the derivatives 519 * are sorted in incresing free parameter order (i.e. f at index 0, df/dx<sub>1</sub> 520 * at index 1, df/dx<sub>2</sub> at index 2 ... df/dx<sub>k</sub> at index k),</li> 521 * <li>all other cases are not publicly specified</li> 522 * </ul> 523 * <p> 524 * This method is the inverse of method {@link #getPartialDerivativeOrders(int)} 525 * </p> 526 * @param orders derivation orders with respect to each parameter 527 * @return index of the partial derivative 528 * @exception DimensionMismatchException if the numbers of parameters does not 529 * match the instance 530 * @exception NumberIsTooLargeException if sum of derivation orders is larger 531 * than the instance limits 532 * @see #getPartialDerivativeOrders(int) 533 */ 534 public int getPartialDerivativeIndex(final int ... orders) 535 throws DimensionMismatchException, NumberIsTooLargeException { 536 537 // safety check 538 if (orders.length != getFreeParameters()) { 539 throw new DimensionMismatchException(orders.length, getFreeParameters()); 540 } 541 542 return getPartialDerivativeIndex(parameters, order, sizes, orders); 543 544 } 545 546 /** Get the index of a partial derivative in an array. 547 * @param parameters number of free parameters 548 * @param order derivation order 549 * @param sizes sizes array 550 * @param orders derivation orders with respect to each parameter 551 * (the lenght of this array must match the number of parameters) 552 * @return index of the partial derivative 553 * @exception NumberIsTooLargeException if sum of derivation orders is larger 554 * than the instance limits 555 */ 556 private static int getPartialDerivativeIndex(final int parameters, final int order, 557 final int[][] sizes, final int ... orders) 558 throws NumberIsTooLargeException { 559 560 // the value is obtained by diving into the recursive Dan Kalman's structure 561 // this is theorem 2 of his paper, with recursion replaced by iteration 562 int index = 0; 563 int m = order; 564 int ordersSum = 0; 565 for (int i = parameters - 1; i >= 0; --i) { 566 567 // derivative order for current free parameter 568 int derivativeOrder = orders[i]; 569 570 // safety check 571 ordersSum += derivativeOrder; 572 if (ordersSum > order) { 573 throw new NumberIsTooLargeException(ordersSum, order, true); 574 } 575 576 while (derivativeOrder-- > 0) { 577 // as long as we differentiate according to current free parameter, 578 // we have to skip the value part and dive into the derivative part 579 // so we add the size of the value part to the base index 580 index += sizes[i][m--]; 581 } 582 583 } 584 585 return index; 586 587 } 588 589 /** Convert an index from one (parameters, order) structure to another. 590 * @param index index of a partial derivative in source derivative structure 591 * @param srcP number of free parameters in source derivative structure 592 * @param srcDerivativesIndirection derivatives indirection array for the source 593 * derivative structure 594 * @param destP number of free parameters in destination derivative structure 595 * @param destO derivation order in destination derivative structure 596 * @param destSizes sizes array for the destination derivative structure 597 * @return index of the partial derivative with the <em>same</em> characteristics 598 * in destination derivative structure 599 */ 600 private static int convertIndex(final int index, 601 final int srcP, final int[][] srcDerivativesIndirection, 602 final int destP, final int destO, final int[][] destSizes) { 603 int[] orders = new int[destP]; 604 System.arraycopy(srcDerivativesIndirection[index], 0, orders, 0, FastMath.min(srcP, destP)); 605 return getPartialDerivativeIndex(destP, destO, destSizes, orders); 606 } 607 608 /** Get the derivation orders for a specific index in the array. 609 * <p> 610 * This method is the inverse of {@link #getPartialDerivativeIndex(int...)}. 611 * </p> 612 * @param index of the partial derivative 613 * @return orders derivation orders with respect to each parameter 614 * @see #getPartialDerivativeIndex(int...) 615 */ 616 public int[] getPartialDerivativeOrders(final int index) { 617 return derivativesIndirection[index]; 618 } 619 620 /** Get the number of free parameters. 621 * @return number of free parameters 622 */ 623 public int getFreeParameters() { 624 return parameters; 625 } 626 627 /** Get the derivation order. 628 * @return derivation order 629 */ 630 public int getOrder() { 631 return order; 632 } 633 634 /** Get the array size required for holding partial derivatives data. 635 * <p> 636 * This number includes the single 0 order derivative element, which is 637 * guaranteed to be stored in the first element of the array. 638 * </p> 639 * @return array size required for holding partial derivatives data 640 */ 641 public int getSize() { 642 return sizes[parameters][order]; 643 } 644 645 /** Compute linear combination. 646 * The derivative structure built will be a1 * ds1 + a2 * ds2 647 * @param a1 first scale factor 648 * @param c1 first base (unscaled) component 649 * @param offset1 offset of first operand in its array 650 * @param a2 second scale factor 651 * @param c2 second base (unscaled) component 652 * @param offset2 offset of second operand in its array 653 * @param result array where result must be stored (it may be 654 * one of the input arrays) 655 * @param resultOffset offset of the result in its array 656 */ 657 public void linearCombination(final double a1, final double[] c1, final int offset1, 658 final double a2, final double[] c2, final int offset2, 659 final double[] result, final int resultOffset) { 660 for (int i = 0; i < getSize(); ++i) { 661 result[resultOffset + i] = 662 MathArrays.linearCombination(a1, c1[offset1 + i], a2, c2[offset2 + i]); 663 } 664 } 665 666 /** Compute linear combination. 667 * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4 668 * @param a1 first scale factor 669 * @param c1 first base (unscaled) component 670 * @param offset1 offset of first operand in its array 671 * @param a2 second scale factor 672 * @param c2 second base (unscaled) component 673 * @param offset2 offset of second operand in its array 674 * @param a3 third scale factor 675 * @param c3 third base (unscaled) component 676 * @param offset3 offset of third operand in its array 677 * @param result array where result must be stored (it may be 678 * one of the input arrays) 679 * @param resultOffset offset of the result in its array 680 */ 681 public void linearCombination(final double a1, final double[] c1, final int offset1, 682 final double a2, final double[] c2, final int offset2, 683 final double a3, final double[] c3, final int offset3, 684 final double[] result, final int resultOffset) { 685 for (int i = 0; i < getSize(); ++i) { 686 result[resultOffset + i] = 687 MathArrays.linearCombination(a1, c1[offset1 + i], 688 a2, c2[offset2 + i], 689 a3, c3[offset3 + i]); 690 } 691 } 692 693 /** Compute linear combination. 694 * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4 695 * @param a1 first scale factor 696 * @param c1 first base (unscaled) component 697 * @param offset1 offset of first operand in its array 698 * @param a2 second scale factor 699 * @param c2 second base (unscaled) component 700 * @param offset2 offset of second operand in its array 701 * @param a3 third scale factor 702 * @param c3 third base (unscaled) component 703 * @param offset3 offset of third operand in its array 704 * @param a4 fourth scale factor 705 * @param c4 fourth base (unscaled) component 706 * @param offset4 offset of fourth operand in its array 707 * @param result array where result must be stored (it may be 708 * one of the input arrays) 709 * @param resultOffset offset of the result in its array 710 */ 711 public void linearCombination(final double a1, final double[] c1, final int offset1, 712 final double a2, final double[] c2, final int offset2, 713 final double a3, final double[] c3, final int offset3, 714 final double a4, final double[] c4, final int offset4, 715 final double[] result, final int resultOffset) { 716 for (int i = 0; i < getSize(); ++i) { 717 result[resultOffset + i] = 718 MathArrays.linearCombination(a1, c1[offset1 + i], 719 a2, c2[offset2 + i], 720 a3, c3[offset3 + i], 721 a4, c4[offset4 + i]); 722 } 723 } 724 725 /** Perform addition of two derivative structures. 726 * @param lhs array holding left hand side of addition 727 * @param lhsOffset offset of the left hand side in its array 728 * @param rhs array right hand side of addition 729 * @param rhsOffset offset of the right hand side in its array 730 * @param result array where result must be stored (it may be 731 * one of the input arrays) 732 * @param resultOffset offset of the result in its array 733 */ 734 public void add(final double[] lhs, final int lhsOffset, 735 final double[] rhs, final int rhsOffset, 736 final double[] result, final int resultOffset) { 737 for (int i = 0; i < getSize(); ++i) { 738 result[resultOffset + i] = lhs[lhsOffset + i] + rhs[rhsOffset + i]; 739 } 740 } 741 /** Perform subtraction of two derivative structures. 742 * @param lhs array holding left hand side of subtraction 743 * @param lhsOffset offset of the left hand side in its array 744 * @param rhs array right hand side of subtraction 745 * @param rhsOffset offset of the right hand side in its array 746 * @param result array where result must be stored (it may be 747 * one of the input arrays) 748 * @param resultOffset offset of the result in its array 749 */ 750 public void subtract(final double[] lhs, final int lhsOffset, 751 final double[] rhs, final int rhsOffset, 752 final double[] result, final int resultOffset) { 753 for (int i = 0; i < getSize(); ++i) { 754 result[resultOffset + i] = lhs[lhsOffset + i] - rhs[rhsOffset + i]; 755 } 756 } 757 758 /** Perform multiplication of two derivative structures. 759 * @param lhs array holding left hand side of multiplication 760 * @param lhsOffset offset of the left hand side in its array 761 * @param rhs array right hand side of multiplication 762 * @param rhsOffset offset of the right hand side in its array 763 * @param result array where result must be stored (for 764 * multiplication the result array <em>cannot</em> be one of 765 * the input arrays) 766 * @param resultOffset offset of the result in its array 767 */ 768 public void multiply(final double[] lhs, final int lhsOffset, 769 final double[] rhs, final int rhsOffset, 770 final double[] result, final int resultOffset) { 771 for (int i = 0; i < multIndirection.length; ++i) { 772 final int[][] mappingI = multIndirection[i]; 773 double r = 0; 774 for (int j = 0; j < mappingI.length; ++j) { 775 r += mappingI[j][0] * 776 lhs[lhsOffset + mappingI[j][1]] * 777 rhs[rhsOffset + mappingI[j][2]]; 778 } 779 result[resultOffset + i] = r; 780 } 781 } 782 783 /** Perform division of two derivative structures. 784 * @param lhs array holding left hand side of division 785 * @param lhsOffset offset of the left hand side in its array 786 * @param rhs array right hand side of division 787 * @param rhsOffset offset of the right hand side in its array 788 * @param result array where result must be stored (for 789 * division the result array <em>cannot</em> be one of 790 * the input arrays) 791 * @param resultOffset offset of the result in its array 792 */ 793 public void divide(final double[] lhs, final int lhsOffset, 794 final double[] rhs, final int rhsOffset, 795 final double[] result, final int resultOffset) { 796 final double[] reciprocal = new double[getSize()]; 797 pow(rhs, lhsOffset, -1, reciprocal, 0); 798 multiply(lhs, lhsOffset, reciprocal, 0, result, resultOffset); 799 } 800 801 /** Perform remainder of two derivative structures. 802 * @param lhs array holding left hand side of remainder 803 * @param lhsOffset offset of the left hand side in its array 804 * @param rhs array right hand side of remainder 805 * @param rhsOffset offset of the right hand side in its array 806 * @param result array where result must be stored (it may be 807 * one of the input arrays) 808 * @param resultOffset offset of the result in its array 809 */ 810 public void remainder(final double[] lhs, final int lhsOffset, 811 final double[] rhs, final int rhsOffset, 812 final double[] result, final int resultOffset) { 813 814 // compute k such that lhs % rhs = lhs - k rhs 815 final double rem = lhs[lhsOffset] % rhs[rhsOffset]; 816 final double k = FastMath.rint((lhs[lhsOffset] - rem) / rhs[rhsOffset]); 817 818 // set up value 819 result[resultOffset] = rem; 820 821 // set up partial derivatives 822 for (int i = 1; i < getSize(); ++i) { 823 result[resultOffset + i] = lhs[lhsOffset + i] - k * rhs[rhsOffset + i]; 824 } 825 826 } 827 828 /** Compute power of a derivative structure. 829 * @param operand array holding the operand 830 * @param operandOffset offset of the operand in its array 831 * @param p power to apply 832 * @param result array where result must be stored (for 833 * power the result array <em>cannot</em> be the input 834 * array) 835 * @param resultOffset offset of the result in its array 836 */ 837 public void pow(final double[] operand, final int operandOffset, final double p, 838 final double[] result, final int resultOffset) { 839 840 // create the function value and derivatives 841 // [x^p, px^(p-1), p(p-1)x^(p-2), ... ] 842 double[] function = new double[1 + order]; 843 double xk = FastMath.pow(operand[operandOffset], p - order); 844 for (int i = order; i > 0; --i) { 845 function[i] = xk; 846 xk *= operand[operandOffset]; 847 } 848 function[0] = xk; 849 double coefficient = p; 850 for (int i = 1; i <= order; ++i) { 851 function[i] *= coefficient; 852 coefficient *= p - i; 853 } 854 855 // apply function composition 856 compose(operand, operandOffset, function, result, resultOffset); 857 858 } 859 860 /** Compute integer power of a derivative structure. 861 * @param operand array holding the operand 862 * @param operandOffset offset of the operand in its array 863 * @param n power to apply 864 * @param result array where result must be stored (for 865 * power the result array <em>cannot</em> be the input 866 * array) 867 * @param resultOffset offset of the result in its array 868 */ 869 public void pow(final double[] operand, final int operandOffset, final int n, 870 final double[] result, final int resultOffset) { 871 872 if (n == 0) { 873 // special case, x^0 = 1 for all x 874 result[resultOffset] = 1.0; 875 Arrays.fill(result, resultOffset + 1, resultOffset + getSize(), 0); 876 return; 877 } 878 879 // create the power function value and derivatives 880 // [x^n, nx^(n-1), n(n-1)x^(n-2), ... ] 881 double[] function = new double[1 + order]; 882 883 if (n > 0) { 884 // strictly positive power 885 final int maxOrder = FastMath.min(order, n); 886 double xk = FastMath.pow(operand[operandOffset], n - maxOrder); 887 for (int i = maxOrder; i > 0; --i) { 888 function[i] = xk; 889 xk *= operand[operandOffset]; 890 } 891 function[0] = xk; 892 } else { 893 // strictly negative power 894 final double inv = 1.0 / operand[operandOffset]; 895 double xk = FastMath.pow(inv, -n); 896 for (int i = 0; i <= order; ++i) { 897 function[i] = xk; 898 xk *= inv; 899 } 900 } 901 902 double coefficient = n; 903 for (int i = 1; i <= order; ++i) { 904 function[i] *= coefficient; 905 coefficient *= n - i; 906 } 907 908 // apply function composition 909 compose(operand, operandOffset, function, result, resultOffset); 910 911 } 912 913 /** Compute power of a derivative structure. 914 * @param x array holding the base 915 * @param xOffset offset of the base in its array 916 * @param y array holding the exponent 917 * @param yOffset offset of the exponent in its array 918 * @param result array where result must be stored (for 919 * power the result array <em>cannot</em> be the input 920 * array) 921 * @param resultOffset offset of the result in its array 922 */ 923 public void pow(final double[] x, final int xOffset, 924 final double[] y, final int yOffset, 925 final double[] result, final int resultOffset) { 926 final double[] logX = new double[getSize()]; 927 log(x, xOffset, logX, 0); 928 final double[] yLogX = new double[getSize()]; 929 multiply(logX, 0, y, yOffset, yLogX, 0); 930 exp(yLogX, 0, result, resultOffset); 931 } 932 933 /** Compute n<sup>th</sup> root of a derivative structure. 934 * @param operand array holding the operand 935 * @param operandOffset offset of the operand in its array 936 * @param n order of the root 937 * @param result array where result must be stored (for 938 * n<sup>th</sup> root the result array <em>cannot</em> be the input 939 * array) 940 * @param resultOffset offset of the result in its array 941 */ 942 public void rootN(final double[] operand, final int operandOffset, final int n, 943 final double[] result, final int resultOffset) { 944 945 // create the function value and derivatives 946 // [x^(1/n), (1/n)x^((1/n)-1), (1-n)/n^2x^((1/n)-2), ... ] 947 double[] function = new double[1 + order]; 948 double xk; 949 if (n == 2) { 950 function[0] = FastMath.sqrt(operand[operandOffset]); 951 xk = 0.5 / function[0]; 952 } else if (n == 3) { 953 function[0] = FastMath.cbrt(operand[operandOffset]); 954 xk = 1.0 / (3.0 * function[0] * function[0]); 955 } else { 956 function[0] = FastMath.pow(operand[operandOffset], 1.0 / n); 957 xk = 1.0 / (n * FastMath.pow(function[0], n - 1)); 958 } 959 final double nReciprocal = 1.0 / n; 960 final double xReciprocal = 1.0 / operand[operandOffset]; 961 for (int i = 1; i <= order; ++i) { 962 function[i] = xk; 963 xk *= xReciprocal * (nReciprocal - i); 964 } 965 966 // apply function composition 967 compose(operand, operandOffset, function, result, resultOffset); 968 969 } 970 971 /** Compute exponential of a derivative structure. 972 * @param operand array holding the operand 973 * @param operandOffset offset of the operand in its array 974 * @param result array where result must be stored (for 975 * exponential the result array <em>cannot</em> be the input 976 * array) 977 * @param resultOffset offset of the result in its array 978 */ 979 public void exp(final double[] operand, final int operandOffset, 980 final double[] result, final int resultOffset) { 981 982 // create the function value and derivatives 983 double[] function = new double[1 + order]; 984 Arrays.fill(function, FastMath.exp(operand[operandOffset])); 985 986 // apply function composition 987 compose(operand, operandOffset, function, result, resultOffset); 988 989 } 990 991 /** Compute exp(x) - 1 of a derivative structure. 992 * @param operand array holding the operand 993 * @param operandOffset offset of the operand in its array 994 * @param result array where result must be stored (for 995 * exponential the result array <em>cannot</em> be the input 996 * array) 997 * @param resultOffset offset of the result in its array 998 */ 999 public void expm1(final double[] operand, final int operandOffset, 1000 final double[] result, final int resultOffset) { 1001 1002 // create the function value and derivatives 1003 double[] function = new double[1 + order]; 1004 function[0] = FastMath.expm1(operand[operandOffset]); 1005 Arrays.fill(function, 1, 1 + order, FastMath.exp(operand[operandOffset])); 1006 1007 // apply function composition 1008 compose(operand, operandOffset, function, result, resultOffset); 1009 1010 } 1011 1012 /** Compute natural logarithm of a derivative structure. 1013 * @param operand array holding the operand 1014 * @param operandOffset offset of the operand in its array 1015 * @param result array where result must be stored (for 1016 * logarithm the result array <em>cannot</em> be the input 1017 * array) 1018 * @param resultOffset offset of the result in its array 1019 */ 1020 public void log(final double[] operand, final int operandOffset, 1021 final double[] result, final int resultOffset) { 1022 1023 // create the function value and derivatives 1024 double[] function = new double[1 + order]; 1025 function[0] = FastMath.log(operand[operandOffset]); 1026 if (order > 0) { 1027 double inv = 1.0 / operand[operandOffset]; 1028 double xk = inv; 1029 for (int i = 1; i <= order; ++i) { 1030 function[i] = xk; 1031 xk *= -i * inv; 1032 } 1033 } 1034 1035 // apply function composition 1036 compose(operand, operandOffset, function, result, resultOffset); 1037 1038 } 1039 1040 /** Computes shifted logarithm of a derivative structure. 1041 * @param operand array holding the operand 1042 * @param operandOffset offset of the operand in its array 1043 * @param result array where result must be stored (for 1044 * shifted logarithm the result array <em>cannot</em> be the input array) 1045 * @param resultOffset offset of the result in its array 1046 */ 1047 public void log1p(final double[] operand, final int operandOffset, 1048 final double[] result, final int resultOffset) { 1049 1050 // create the function value and derivatives 1051 double[] function = new double[1 + order]; 1052 function[0] = FastMath.log1p(operand[operandOffset]); 1053 if (order > 0) { 1054 double inv = 1.0 / (1.0 + operand[operandOffset]); 1055 double xk = inv; 1056 for (int i = 1; i <= order; ++i) { 1057 function[i] = xk; 1058 xk *= -i * inv; 1059 } 1060 } 1061 1062 // apply function composition 1063 compose(operand, operandOffset, function, result, resultOffset); 1064 1065 } 1066 1067 /** Computes base 10 logarithm of a derivative structure. 1068 * @param operand array holding the operand 1069 * @param operandOffset offset of the operand in its array 1070 * @param result array where result must be stored (for 1071 * base 10 logarithm the result array <em>cannot</em> be the input array) 1072 * @param resultOffset offset of the result in its array 1073 */ 1074 public void log10(final double[] operand, final int operandOffset, 1075 final double[] result, final int resultOffset) { 1076 1077 // create the function value and derivatives 1078 double[] function = new double[1 + order]; 1079 function[0] = FastMath.log10(operand[operandOffset]); 1080 if (order > 0) { 1081 double inv = 1.0 / operand[operandOffset]; 1082 double xk = inv / FastMath.log(10.0); 1083 for (int i = 1; i <= order; ++i) { 1084 function[i] = xk; 1085 xk *= -i * inv; 1086 } 1087 } 1088 1089 // apply function composition 1090 compose(operand, operandOffset, function, result, resultOffset); 1091 1092 } 1093 1094 /** Compute cosine of a derivative structure. 1095 * @param operand array holding the operand 1096 * @param operandOffset offset of the operand in its array 1097 * @param result array where result must be stored (for 1098 * cosine the result array <em>cannot</em> be the input 1099 * array) 1100 * @param resultOffset offset of the result in its array 1101 */ 1102 public void cos(final double[] operand, final int operandOffset, 1103 final double[] result, final int resultOffset) { 1104 1105 // create the function value and derivatives 1106 double[] function = new double[1 + order]; 1107 function[0] = FastMath.cos(operand[operandOffset]); 1108 if (order > 0) { 1109 function[1] = -FastMath.sin(operand[operandOffset]); 1110 for (int i = 2; i <= order; ++i) { 1111 function[i] = -function[i - 2]; 1112 } 1113 } 1114 1115 // apply function composition 1116 compose(operand, operandOffset, function, result, resultOffset); 1117 1118 } 1119 1120 /** Compute sine of a derivative structure. 1121 * @param operand array holding the operand 1122 * @param operandOffset offset of the operand in its array 1123 * @param result array where result must be stored (for 1124 * sine the result array <em>cannot</em> be the input 1125 * array) 1126 * @param resultOffset offset of the result in its array 1127 */ 1128 public void sin(final double[] operand, final int operandOffset, 1129 final double[] result, final int resultOffset) { 1130 1131 // create the function value and derivatives 1132 double[] function = new double[1 + order]; 1133 function[0] = FastMath.sin(operand[operandOffset]); 1134 if (order > 0) { 1135 function[1] = FastMath.cos(operand[operandOffset]); 1136 for (int i = 2; i <= order; ++i) { 1137 function[i] = -function[i - 2]; 1138 } 1139 } 1140 1141 // apply function composition 1142 compose(operand, operandOffset, function, result, resultOffset); 1143 1144 } 1145 1146 /** Compute tangent of a derivative structure. 1147 * @param operand array holding the operand 1148 * @param operandOffset offset of the operand in its array 1149 * @param result array where result must be stored (for 1150 * tangent the result array <em>cannot</em> be the input 1151 * array) 1152 * @param resultOffset offset of the result in its array 1153 */ 1154 public void tan(final double[] operand, final int operandOffset, 1155 final double[] result, final int resultOffset) { 1156 1157 // create the function value and derivatives 1158 final double[] function = new double[1 + order]; 1159 final double t = FastMath.tan(operand[operandOffset]); 1160 function[0] = t; 1161 1162 if (order > 0) { 1163 1164 // the nth order derivative of tan has the form: 1165 // dn(tan(x)/dxn = P_n(tan(x)) 1166 // where P_n(t) is a degree n+1 polynomial with same parity as n+1 1167 // P_0(t) = t, P_1(t) = 1 + t^2, P_2(t) = 2 t (1 + t^2) ... 1168 // the general recurrence relation for P_n is: 1169 // P_n(x) = (1+t^2) P_(n-1)'(t) 1170 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array 1171 final double[] p = new double[order + 2]; 1172 p[1] = 1; 1173 final double t2 = t * t; 1174 for (int n = 1; n <= order; ++n) { 1175 1176 // update and evaluate polynomial P_n(t) 1177 double v = 0; 1178 p[n + 1] = n * p[n]; 1179 for (int k = n + 1; k >= 0; k -= 2) { 1180 v = v * t2 + p[k]; 1181 if (k > 2) { 1182 p[k - 2] = (k - 1) * p[k - 1] + (k - 3) * p[k - 3]; 1183 } else if (k == 2) { 1184 p[0] = p[1]; 1185 } 1186 } 1187 if ((n & 0x1) == 0) { 1188 v *= t; 1189 } 1190 1191 function[n] = v; 1192 1193 } 1194 } 1195 1196 // apply function composition 1197 compose(operand, operandOffset, function, result, resultOffset); 1198 1199 } 1200 1201 /** Compute arc cosine of a derivative structure. 1202 * @param operand array holding the operand 1203 * @param operandOffset offset of the operand in its array 1204 * @param result array where result must be stored (for 1205 * arc cosine the result array <em>cannot</em> be the input 1206 * array) 1207 * @param resultOffset offset of the result in its array 1208 */ 1209 public void acos(final double[] operand, final int operandOffset, 1210 final double[] result, final int resultOffset) { 1211 1212 // create the function value and derivatives 1213 double[] function = new double[1 + order]; 1214 final double x = operand[operandOffset]; 1215 function[0] = FastMath.acos(x); 1216 if (order > 0) { 1217 // the nth order derivative of acos has the form: 1218 // dn(acos(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2) 1219 // where P_n(x) is a degree n-1 polynomial with same parity as n-1 1220 // P_1(x) = -1, P_2(x) = -x, P_3(x) = -2x^2 - 1 ... 1221 // the general recurrence relation for P_n is: 1222 // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x) 1223 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array 1224 final double[] p = new double[order]; 1225 p[0] = -1; 1226 final double x2 = x * x; 1227 final double f = 1.0 / (1 - x2); 1228 double coeff = FastMath.sqrt(f); 1229 function[1] = coeff * p[0]; 1230 for (int n = 2; n <= order; ++n) { 1231 1232 // update and evaluate polynomial P_n(x) 1233 double v = 0; 1234 p[n - 1] = (n - 1) * p[n - 2]; 1235 for (int k = n - 1; k >= 0; k -= 2) { 1236 v = v * x2 + p[k]; 1237 if (k > 2) { 1238 p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3]; 1239 } else if (k == 2) { 1240 p[0] = p[1]; 1241 } 1242 } 1243 if ((n & 0x1) == 0) { 1244 v *= x; 1245 } 1246 1247 coeff *= f; 1248 function[n] = coeff * v; 1249 1250 } 1251 } 1252 1253 // apply function composition 1254 compose(operand, operandOffset, function, result, resultOffset); 1255 1256 } 1257 1258 /** Compute arc sine of a derivative structure. 1259 * @param operand array holding the operand 1260 * @param operandOffset offset of the operand in its array 1261 * @param result array where result must be stored (for 1262 * arc sine the result array <em>cannot</em> be the input 1263 * array) 1264 * @param resultOffset offset of the result in its array 1265 */ 1266 public void asin(final double[] operand, final int operandOffset, 1267 final double[] result, final int resultOffset) { 1268 1269 // create the function value and derivatives 1270 double[] function = new double[1 + order]; 1271 final double x = operand[operandOffset]; 1272 function[0] = FastMath.asin(x); 1273 if (order > 0) { 1274 // the nth order derivative of asin has the form: 1275 // dn(asin(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2) 1276 // where P_n(x) is a degree n-1 polynomial with same parity as n-1 1277 // P_1(x) = 1, P_2(x) = x, P_3(x) = 2x^2 + 1 ... 1278 // the general recurrence relation for P_n is: 1279 // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x) 1280 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array 1281 final double[] p = new double[order]; 1282 p[0] = 1; 1283 final double x2 = x * x; 1284 final double f = 1.0 / (1 - x2); 1285 double coeff = FastMath.sqrt(f); 1286 function[1] = coeff * p[0]; 1287 for (int n = 2; n <= order; ++n) { 1288 1289 // update and evaluate polynomial P_n(x) 1290 double v = 0; 1291 p[n - 1] = (n - 1) * p[n - 2]; 1292 for (int k = n - 1; k >= 0; k -= 2) { 1293 v = v * x2 + p[k]; 1294 if (k > 2) { 1295 p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3]; 1296 } else if (k == 2) { 1297 p[0] = p[1]; 1298 } 1299 } 1300 if ((n & 0x1) == 0) { 1301 v *= x; 1302 } 1303 1304 coeff *= f; 1305 function[n] = coeff * v; 1306 1307 } 1308 } 1309 1310 // apply function composition 1311 compose(operand, operandOffset, function, result, resultOffset); 1312 1313 } 1314 1315 /** Compute arc tangent of a derivative structure. 1316 * @param operand array holding the operand 1317 * @param operandOffset offset of the operand in its array 1318 * @param result array where result must be stored (for 1319 * arc tangent the result array <em>cannot</em> be the input 1320 * array) 1321 * @param resultOffset offset of the result in its array 1322 */ 1323 public void atan(final double[] operand, final int operandOffset, 1324 final double[] result, final int resultOffset) { 1325 1326 // create the function value and derivatives 1327 double[] function = new double[1 + order]; 1328 final double x = operand[operandOffset]; 1329 function[0] = FastMath.atan(x); 1330 if (order > 0) { 1331 // the nth order derivative of atan has the form: 1332 // dn(atan(x)/dxn = Q_n(x) / (1 + x^2)^n 1333 // where Q_n(x) is a degree n-1 polynomial with same parity as n-1 1334 // Q_1(x) = 1, Q_2(x) = -2x, Q_3(x) = 6x^2 - 2 ... 1335 // the general recurrence relation for Q_n is: 1336 // Q_n(x) = (1+x^2) Q_(n-1)'(x) - 2(n-1) x Q_(n-1)(x) 1337 // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array 1338 final double[] q = new double[order]; 1339 q[0] = 1; 1340 final double x2 = x * x; 1341 final double f = 1.0 / (1 + x2); 1342 double coeff = f; 1343 function[1] = coeff * q[0]; 1344 for (int n = 2; n <= order; ++n) { 1345 1346 // update and evaluate polynomial Q_n(x) 1347 double v = 0; 1348 q[n - 1] = -n * q[n - 2]; 1349 for (int k = n - 1; k >= 0; k -= 2) { 1350 v = v * x2 + q[k]; 1351 if (k > 2) { 1352 q[k - 2] = (k - 1) * q[k - 1] + (k - 1 - 2 * n) * q[k - 3]; 1353 } else if (k == 2) { 1354 q[0] = q[1]; 1355 } 1356 } 1357 if ((n & 0x1) == 0) { 1358 v *= x; 1359 } 1360 1361 coeff *= f; 1362 function[n] = coeff * v; 1363 1364 } 1365 } 1366 1367 // apply function composition 1368 compose(operand, operandOffset, function, result, resultOffset); 1369 1370 } 1371 1372 /** Compute two arguments arc tangent of a derivative structure. 1373 * @param y array holding the first operand 1374 * @param yOffset offset of the first operand in its array 1375 * @param x array holding the second operand 1376 * @param xOffset offset of the second operand in its array 1377 * @param result array where result must be stored (for 1378 * two arguments arc tangent the result array <em>cannot</em> 1379 * be the input array) 1380 * @param resultOffset offset of the result in its array 1381 */ 1382 public void atan2(final double[] y, final int yOffset, 1383 final double[] x, final int xOffset, 1384 final double[] result, final int resultOffset) { 1385 1386 // compute r = sqrt(x^2+y^2) 1387 double[] tmp1 = new double[getSize()]; 1388 multiply(x, xOffset, x, xOffset, tmp1, 0); // x^2 1389 double[] tmp2 = new double[getSize()]; 1390 multiply(y, yOffset, y, yOffset, tmp2, 0); // y^2 1391 add(tmp1, 0, tmp2, 0, tmp2, 0); // x^2 + y^2 1392 rootN(tmp2, 0, 2, tmp1, 0); // r = sqrt(x^2 + y^2) 1393 1394 if (x[xOffset] >= 0) { 1395 1396 // compute atan2(y, x) = 2 atan(y / (r + x)) 1397 add(tmp1, 0, x, xOffset, tmp2, 0); // r + x 1398 divide(y, yOffset, tmp2, 0, tmp1, 0); // y /(r + x) 1399 atan(tmp1, 0, tmp2, 0); // atan(y / (r + x)) 1400 for (int i = 0; i < tmp2.length; ++i) { 1401 result[resultOffset + i] = 2 * tmp2[i]; // 2 * atan(y / (r + x)) 1402 } 1403 1404 } else { 1405 1406 // compute atan2(y, x) = +/- pi - 2 atan(y / (r - x)) 1407 subtract(tmp1, 0, x, xOffset, tmp2, 0); // r - x 1408 divide(y, yOffset, tmp2, 0, tmp1, 0); // y /(r - x) 1409 atan(tmp1, 0, tmp2, 0); // atan(y / (r - x)) 1410 result[resultOffset] = 1411 ((tmp2[0] <= 0) ? -FastMath.PI : FastMath.PI) - 2 * tmp2[0]; // +/-pi - 2 * atan(y / (r - x)) 1412 for (int i = 1; i < tmp2.length; ++i) { 1413 result[resultOffset + i] = -2 * tmp2[i]; // +/-pi - 2 * atan(y / (r - x)) 1414 } 1415 1416 } 1417 1418 } 1419 1420 /** Compute hyperbolic cosine of a derivative structure. 1421 * @param operand array holding the operand 1422 * @param operandOffset offset of the operand in its array 1423 * @param result array where result must be stored (for 1424 * hyperbolic cosine the result array <em>cannot</em> be the input 1425 * array) 1426 * @param resultOffset offset of the result in its array 1427 */ 1428 public void cosh(final double[] operand, final int operandOffset, 1429 final double[] result, final int resultOffset) { 1430 1431 // create the function value and derivatives 1432 double[] function = new double[1 + order]; 1433 function[0] = FastMath.cosh(operand[operandOffset]); 1434 if (order > 0) { 1435 function[1] = FastMath.sinh(operand[operandOffset]); 1436 for (int i = 2; i <= order; ++i) { 1437 function[i] = function[i - 2]; 1438 } 1439 } 1440 1441 // apply function composition 1442 compose(operand, operandOffset, function, result, resultOffset); 1443 1444 } 1445 1446 /** Compute hyperbolic sine of a derivative structure. 1447 * @param operand array holding the operand 1448 * @param operandOffset offset of the operand in its array 1449 * @param result array where result must be stored (for 1450 * hyperbolic sine the result array <em>cannot</em> be the input 1451 * array) 1452 * @param resultOffset offset of the result in its array 1453 */ 1454 public void sinh(final double[] operand, final int operandOffset, 1455 final double[] result, final int resultOffset) { 1456 1457 // create the function value and derivatives 1458 double[] function = new double[1 + order]; 1459 function[0] = FastMath.sinh(operand[operandOffset]); 1460 if (order > 0) { 1461 function[1] = FastMath.cosh(operand[operandOffset]); 1462 for (int i = 2; i <= order; ++i) { 1463 function[i] = function[i - 2]; 1464 } 1465 } 1466 1467 // apply function composition 1468 compose(operand, operandOffset, function, result, resultOffset); 1469 1470 } 1471 1472 /** Compute hyperbolic tangent of a derivative structure. 1473 * @param operand array holding the operand 1474 * @param operandOffset offset of the operand in its array 1475 * @param result array where result must be stored (for 1476 * hyperbolic tangent the result array <em>cannot</em> be the input 1477 * array) 1478 * @param resultOffset offset of the result in its array 1479 */ 1480 public void tanh(final double[] operand, final int operandOffset, 1481 final double[] result, final int resultOffset) { 1482 1483 // create the function value and derivatives 1484 final double[] function = new double[1 + order]; 1485 final double t = FastMath.tanh(operand[operandOffset]); 1486 function[0] = t; 1487 1488 if (order > 0) { 1489 1490 // the nth order derivative of tanh has the form: 1491 // dn(tanh(x)/dxn = P_n(tanh(x)) 1492 // where P_n(t) is a degree n+1 polynomial with same parity as n+1 1493 // P_0(t) = t, P_1(t) = 1 - t^2, P_2(t) = -2 t (1 - t^2) ... 1494 // the general recurrence relation for P_n is: 1495 // P_n(x) = (1-t^2) P_(n-1)'(t) 1496 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array 1497 final double[] p = new double[order + 2]; 1498 p[1] = 1; 1499 final double t2 = t * t; 1500 for (int n = 1; n <= order; ++n) { 1501 1502 // update and evaluate polynomial P_n(t) 1503 double v = 0; 1504 p[n + 1] = -n * p[n]; 1505 for (int k = n + 1; k >= 0; k -= 2) { 1506 v = v * t2 + p[k]; 1507 if (k > 2) { 1508 p[k - 2] = (k - 1) * p[k - 1] - (k - 3) * p[k - 3]; 1509 } else if (k == 2) { 1510 p[0] = p[1]; 1511 } 1512 } 1513 if ((n & 0x1) == 0) { 1514 v *= t; 1515 } 1516 1517 function[n] = v; 1518 1519 } 1520 } 1521 1522 // apply function composition 1523 compose(operand, operandOffset, function, result, resultOffset); 1524 1525 } 1526 1527 /** Compute inverse hyperbolic cosine of a derivative structure. 1528 * @param operand array holding the operand 1529 * @param operandOffset offset of the operand in its array 1530 * @param result array where result must be stored (for 1531 * inverse hyperbolic cosine the result array <em>cannot</em> be the input 1532 * array) 1533 * @param resultOffset offset of the result in its array 1534 */ 1535 public void acosh(final double[] operand, final int operandOffset, 1536 final double[] result, final int resultOffset) { 1537 1538 // create the function value and derivatives 1539 double[] function = new double[1 + order]; 1540 final double x = operand[operandOffset]; 1541 function[0] = FastMath.acosh(x); 1542 if (order > 0) { 1543 // the nth order derivative of acosh has the form: 1544 // dn(acosh(x)/dxn = P_n(x) / [x^2 - 1]^((2n-1)/2) 1545 // where P_n(x) is a degree n-1 polynomial with same parity as n-1 1546 // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 + 1 ... 1547 // the general recurrence relation for P_n is: 1548 // P_n(x) = (x^2-1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x) 1549 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array 1550 final double[] p = new double[order]; 1551 p[0] = 1; 1552 final double x2 = x * x; 1553 final double f = 1.0 / (x2 - 1); 1554 double coeff = FastMath.sqrt(f); 1555 function[1] = coeff * p[0]; 1556 for (int n = 2; n <= order; ++n) { 1557 1558 // update and evaluate polynomial P_n(x) 1559 double v = 0; 1560 p[n - 1] = (1 - n) * p[n - 2]; 1561 for (int k = n - 1; k >= 0; k -= 2) { 1562 v = v * x2 + p[k]; 1563 if (k > 2) { 1564 p[k - 2] = (1 - k) * p[k - 1] + (k - 2 * n) * p[k - 3]; 1565 } else if (k == 2) { 1566 p[0] = -p[1]; 1567 } 1568 } 1569 if ((n & 0x1) == 0) { 1570 v *= x; 1571 } 1572 1573 coeff *= f; 1574 function[n] = coeff * v; 1575 1576 } 1577 } 1578 1579 // apply function composition 1580 compose(operand, operandOffset, function, result, resultOffset); 1581 1582 } 1583 1584 /** Compute inverse hyperbolic sine of a derivative structure. 1585 * @param operand array holding the operand 1586 * @param operandOffset offset of the operand in its array 1587 * @param result array where result must be stored (for 1588 * inverse hyperbolic sine the result array <em>cannot</em> be the input 1589 * array) 1590 * @param resultOffset offset of the result in its array 1591 */ 1592 public void asinh(final double[] operand, final int operandOffset, 1593 final double[] result, final int resultOffset) { 1594 1595 // create the function value and derivatives 1596 double[] function = new double[1 + order]; 1597 final double x = operand[operandOffset]; 1598 function[0] = FastMath.asinh(x); 1599 if (order > 0) { 1600 // the nth order derivative of asinh has the form: 1601 // dn(asinh(x)/dxn = P_n(x) / [x^2 + 1]^((2n-1)/2) 1602 // where P_n(x) is a degree n-1 polynomial with same parity as n-1 1603 // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 - 1 ... 1604 // the general recurrence relation for P_n is: 1605 // P_n(x) = (x^2+1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x) 1606 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array 1607 final double[] p = new double[order]; 1608 p[0] = 1; 1609 final double x2 = x * x; 1610 final double f = 1.0 / (1 + x2); 1611 double coeff = FastMath.sqrt(f); 1612 function[1] = coeff * p[0]; 1613 for (int n = 2; n <= order; ++n) { 1614 1615 // update and evaluate polynomial P_n(x) 1616 double v = 0; 1617 p[n - 1] = (1 - n) * p[n - 2]; 1618 for (int k = n - 1; k >= 0; k -= 2) { 1619 v = v * x2 + p[k]; 1620 if (k > 2) { 1621 p[k - 2] = (k - 1) * p[k - 1] + (k - 2 * n) * p[k - 3]; 1622 } else if (k == 2) { 1623 p[0] = p[1]; 1624 } 1625 } 1626 if ((n & 0x1) == 0) { 1627 v *= x; 1628 } 1629 1630 coeff *= f; 1631 function[n] = coeff * v; 1632 1633 } 1634 } 1635 1636 // apply function composition 1637 compose(operand, operandOffset, function, result, resultOffset); 1638 1639 } 1640 1641 /** Compute inverse hyperbolic tangent of a derivative structure. 1642 * @param operand array holding the operand 1643 * @param operandOffset offset of the operand in its array 1644 * @param result array where result must be stored (for 1645 * inverse hyperbolic tangent the result array <em>cannot</em> be the input 1646 * array) 1647 * @param resultOffset offset of the result in its array 1648 */ 1649 public void atanh(final double[] operand, final int operandOffset, 1650 final double[] result, final int resultOffset) { 1651 1652 // create the function value and derivatives 1653 double[] function = new double[1 + order]; 1654 final double x = operand[operandOffset]; 1655 function[0] = FastMath.atanh(x); 1656 if (order > 0) { 1657 // the nth order derivative of atanh has the form: 1658 // dn(atanh(x)/dxn = Q_n(x) / (1 - x^2)^n 1659 // where Q_n(x) is a degree n-1 polynomial with same parity as n-1 1660 // Q_1(x) = 1, Q_2(x) = 2x, Q_3(x) = 6x^2 + 2 ... 1661 // the general recurrence relation for Q_n is: 1662 // Q_n(x) = (1-x^2) Q_(n-1)'(x) + 2(n-1) x Q_(n-1)(x) 1663 // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array 1664 final double[] q = new double[order]; 1665 q[0] = 1; 1666 final double x2 = x * x; 1667 final double f = 1.0 / (1 - x2); 1668 double coeff = f; 1669 function[1] = coeff * q[0]; 1670 for (int n = 2; n <= order; ++n) { 1671 1672 // update and evaluate polynomial Q_n(x) 1673 double v = 0; 1674 q[n - 1] = n * q[n - 2]; 1675 for (int k = n - 1; k >= 0; k -= 2) { 1676 v = v * x2 + q[k]; 1677 if (k > 2) { 1678 q[k - 2] = (k - 1) * q[k - 1] + (2 * n - k + 1) * q[k - 3]; 1679 } else if (k == 2) { 1680 q[0] = q[1]; 1681 } 1682 } 1683 if ((n & 0x1) == 0) { 1684 v *= x; 1685 } 1686 1687 coeff *= f; 1688 function[n] = coeff * v; 1689 1690 } 1691 } 1692 1693 // apply function composition 1694 compose(operand, operandOffset, function, result, resultOffset); 1695 1696 } 1697 1698 /** Compute composition of a derivative structure by a function. 1699 * @param operand array holding the operand 1700 * @param operandOffset offset of the operand in its array 1701 * @param f array of value and derivatives of the function at 1702 * the current point (i.e. at {@code operand[operandOffset]}). 1703 * @param result array where result must be stored (for 1704 * composition the result array <em>cannot</em> be the input 1705 * array) 1706 * @param resultOffset offset of the result in its array 1707 */ 1708 public void compose(final double[] operand, final int operandOffset, final double[] f, 1709 final double[] result, final int resultOffset) { 1710 for (int i = 0; i < compIndirection.length; ++i) { 1711 final int[][] mappingI = compIndirection[i]; 1712 double r = 0; 1713 for (int j = 0; j < mappingI.length; ++j) { 1714 final int[] mappingIJ = mappingI[j]; 1715 double product = mappingIJ[0] * f[mappingIJ[1]]; 1716 for (int k = 2; k < mappingIJ.length; ++k) { 1717 product *= operand[operandOffset + mappingIJ[k]]; 1718 } 1719 r += product; 1720 } 1721 result[resultOffset + i] = r; 1722 } 1723 } 1724 1725 /** Evaluate Taylor expansion of a derivative structure. 1726 * @param ds array holding the derivative structure 1727 * @param dsOffset offset of the derivative structure in its array 1728 * @param delta parameters offsets (Δx, Δy, ...) 1729 * @return value of the Taylor expansion at x + Δx, y + Δy, ... 1730 */ 1731 public double taylor(final double[] ds, final int dsOffset, final double ... delta) { 1732 double value = 0; 1733 for (int i = getSize() - 1; i >= 0; --i) { 1734 final int[] orders = getPartialDerivativeOrders(i); 1735 double term = ds[dsOffset + i]; 1736 for (int k = 0; k < orders.length; ++k) { 1737 if (orders[k] > 0) { 1738 term *= FastMath.pow(delta[k], orders[k]) / ArithmeticUtils.factorial(orders[k]); 1739 } 1740 } 1741 value += term; 1742 } 1743 return value; 1744 } 1745 1746 /** Check rules set compatibility. 1747 * @param compiler other compiler to check against instance 1748 * @exception DimensionMismatchException if number of free parameters or orders are inconsistent 1749 */ 1750 public void checkCompatibility(final DSCompiler compiler) 1751 throws DimensionMismatchException { 1752 if (parameters != compiler.parameters) { 1753 throw new DimensionMismatchException(parameters, compiler.parameters); 1754 } 1755 if (order != compiler.order) { 1756 throw new DimensionMismatchException(order, compiler.order); 1757 } 1758 } 1759 1760 }