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.dfp; 018 019 020 import org.apache.commons.math3.analysis.solvers.AllowedSolution; 021 import org.apache.commons.math3.exception.MathInternalError; 022 import org.apache.commons.math3.exception.NoBracketingException; 023 import org.apache.commons.math3.exception.NullArgumentException; 024 import org.apache.commons.math3.exception.NumberIsTooSmallException; 025 import org.apache.commons.math3.util.Incrementor; 026 import org.apache.commons.math3.util.MathUtils; 027 028 /** 029 * This class implements a modification of the <a 030 * href="http://mathworld.wolfram.com/BrentsMethod.html"> Brent algorithm</a>. 031 * <p> 032 * The changes with respect to the original Brent algorithm are: 033 * <ul> 034 * <li>the returned value is chosen in the current interval according 035 * to user specified {@link AllowedSolution},</li> 036 * <li>the maximal order for the invert polynomial root search is 037 * user-specified instead of being invert quadratic only</li> 038 * </ul> 039 * </p> 040 * The given interval must bracket the root. 041 * 042 * @version $Id: BracketingNthOrderBrentSolverDFP.java 1416643 2012-12-03 19:37:14Z tn $ 043 */ 044 public class BracketingNthOrderBrentSolverDFP { 045 046 /** Maximal aging triggering an attempt to balance the bracketing interval. */ 047 private static final int MAXIMAL_AGING = 2; 048 049 /** Maximal order. */ 050 private final int maximalOrder; 051 052 /** Function value accuracy. */ 053 private final Dfp functionValueAccuracy; 054 055 /** Absolute accuracy. */ 056 private final Dfp absoluteAccuracy; 057 058 /** Relative accuracy. */ 059 private final Dfp relativeAccuracy; 060 061 /** Evaluations counter. */ 062 private final Incrementor evaluations = new Incrementor(); 063 064 /** 065 * Construct a solver. 066 * 067 * @param relativeAccuracy Relative accuracy. 068 * @param absoluteAccuracy Absolute accuracy. 069 * @param functionValueAccuracy Function value accuracy. 070 * @param maximalOrder maximal order. 071 * @exception NumberIsTooSmallException if maximal order is lower than 2 072 */ 073 public BracketingNthOrderBrentSolverDFP(final Dfp relativeAccuracy, 074 final Dfp absoluteAccuracy, 075 final Dfp functionValueAccuracy, 076 final int maximalOrder) 077 throws NumberIsTooSmallException { 078 if (maximalOrder < 2) { 079 throw new NumberIsTooSmallException(maximalOrder, 2, true); 080 } 081 this.maximalOrder = maximalOrder; 082 this.absoluteAccuracy = absoluteAccuracy; 083 this.relativeAccuracy = relativeAccuracy; 084 this.functionValueAccuracy = functionValueAccuracy; 085 } 086 087 /** Get the maximal order. 088 * @return maximal order 089 */ 090 public int getMaximalOrder() { 091 return maximalOrder; 092 } 093 094 /** 095 * Get the maximal number of function evaluations. 096 * 097 * @return the maximal number of function evaluations. 098 */ 099 public int getMaxEvaluations() { 100 return evaluations.getMaximalCount(); 101 } 102 103 /** 104 * Get the number of evaluations of the objective function. 105 * The number of evaluations corresponds to the last call to the 106 * {@code optimize} method. It is 0 if the method has not been 107 * called yet. 108 * 109 * @return the number of evaluations of the objective function. 110 */ 111 public int getEvaluations() { 112 return evaluations.getCount(); 113 } 114 115 /** 116 * Get the absolute accuracy. 117 * @return absolute accuracy 118 */ 119 public Dfp getAbsoluteAccuracy() { 120 return absoluteAccuracy; 121 } 122 123 /** 124 * Get the relative accuracy. 125 * @return relative accuracy 126 */ 127 public Dfp getRelativeAccuracy() { 128 return relativeAccuracy; 129 } 130 131 /** 132 * Get the function accuracy. 133 * @return function accuracy 134 */ 135 public Dfp getFunctionValueAccuracy() { 136 return functionValueAccuracy; 137 } 138 139 /** 140 * Solve for a zero in the given interval. 141 * A solver may require that the interval brackets a single zero root. 142 * Solvers that do require bracketing should be able to handle the case 143 * where one of the endpoints is itself a root. 144 * 145 * @param maxEval Maximum number of evaluations. 146 * @param f Function to solve. 147 * @param min Lower bound for the interval. 148 * @param max Upper bound for the interval. 149 * @param allowedSolution The kind of solutions that the root-finding algorithm may 150 * accept as solutions. 151 * @return a value where the function is zero. 152 * @exception NullArgumentException if f is null. 153 * @exception NoBracketingException if root cannot be bracketed 154 */ 155 public Dfp solve(final int maxEval, final UnivariateDfpFunction f, 156 final Dfp min, final Dfp max, final AllowedSolution allowedSolution) 157 throws NullArgumentException, NoBracketingException { 158 return solve(maxEval, f, min, max, min.add(max).divide(2), allowedSolution); 159 } 160 161 /** 162 * Solve for a zero in the given interval, start at {@code startValue}. 163 * A solver may require that the interval brackets a single zero root. 164 * Solvers that do require bracketing should be able to handle the case 165 * where one of the endpoints is itself a root. 166 * 167 * @param maxEval Maximum number of evaluations. 168 * @param f Function to solve. 169 * @param min Lower bound for the interval. 170 * @param max Upper bound for the interval. 171 * @param startValue Start value to use. 172 * @param allowedSolution The kind of solutions that the root-finding algorithm may 173 * accept as solutions. 174 * @return a value where the function is zero. 175 * @exception NullArgumentException if f is null. 176 * @exception NoBracketingException if root cannot be bracketed 177 */ 178 public Dfp solve(final int maxEval, final UnivariateDfpFunction f, 179 final Dfp min, final Dfp max, final Dfp startValue, 180 final AllowedSolution allowedSolution) 181 throws NullArgumentException, NoBracketingException { 182 183 // Checks. 184 MathUtils.checkNotNull(f); 185 186 // Reset. 187 evaluations.setMaximalCount(maxEval); 188 evaluations.resetCount(); 189 Dfp zero = startValue.getZero(); 190 Dfp nan = zero.newInstance((byte) 1, Dfp.QNAN); 191 192 // prepare arrays with the first points 193 final Dfp[] x = new Dfp[maximalOrder + 1]; 194 final Dfp[] y = new Dfp[maximalOrder + 1]; 195 x[0] = min; 196 x[1] = startValue; 197 x[2] = max; 198 199 // evaluate initial guess 200 evaluations.incrementCount(); 201 y[1] = f.value(x[1]); 202 if (y[1].isZero()) { 203 // return the initial guess if it is a perfect root. 204 return x[1]; 205 } 206 207 // evaluate first endpoint 208 evaluations.incrementCount(); 209 y[0] = f.value(x[0]); 210 if (y[0].isZero()) { 211 // return the first endpoint if it is a perfect root. 212 return x[0]; 213 } 214 215 int nbPoints; 216 int signChangeIndex; 217 if (y[0].multiply(y[1]).negativeOrNull()) { 218 219 // reduce interval if it brackets the root 220 nbPoints = 2; 221 signChangeIndex = 1; 222 223 } else { 224 225 // evaluate second endpoint 226 evaluations.incrementCount(); 227 y[2] = f.value(x[2]); 228 if (y[2].isZero()) { 229 // return the second endpoint if it is a perfect root. 230 return x[2]; 231 } 232 233 if (y[1].multiply(y[2]).negativeOrNull()) { 234 // use all computed point as a start sampling array for solving 235 nbPoints = 3; 236 signChangeIndex = 2; 237 } else { 238 throw new NoBracketingException(x[0].toDouble(), x[2].toDouble(), 239 y[0].toDouble(), y[2].toDouble()); 240 } 241 242 } 243 244 // prepare a work array for inverse polynomial interpolation 245 final Dfp[] tmpX = new Dfp[x.length]; 246 247 // current tightest bracketing of the root 248 Dfp xA = x[signChangeIndex - 1]; 249 Dfp yA = y[signChangeIndex - 1]; 250 Dfp absXA = xA.abs(); 251 Dfp absYA = yA.abs(); 252 int agingA = 0; 253 Dfp xB = x[signChangeIndex]; 254 Dfp yB = y[signChangeIndex]; 255 Dfp absXB = xB.abs(); 256 Dfp absYB = yB.abs(); 257 int agingB = 0; 258 259 // search loop 260 while (true) { 261 262 // check convergence of bracketing interval 263 Dfp maxX = absXA.lessThan(absXB) ? absXB : absXA; 264 Dfp maxY = absYA.lessThan(absYB) ? absYB : absYA; 265 final Dfp xTol = absoluteAccuracy.add(relativeAccuracy.multiply(maxX)); 266 if (xB.subtract(xA).subtract(xTol).negativeOrNull() || 267 maxY.lessThan(functionValueAccuracy)) { 268 switch (allowedSolution) { 269 case ANY_SIDE : 270 return absYA.lessThan(absYB) ? xA : xB; 271 case LEFT_SIDE : 272 return xA; 273 case RIGHT_SIDE : 274 return xB; 275 case BELOW_SIDE : 276 return yA.lessThan(zero) ? xA : xB; 277 case ABOVE_SIDE : 278 return yA.lessThan(zero) ? xB : xA; 279 default : 280 // this should never happen 281 throw new MathInternalError(null); 282 } 283 } 284 285 // target for the next evaluation point 286 Dfp targetY; 287 if (agingA >= MAXIMAL_AGING) { 288 // we keep updating the high bracket, try to compensate this 289 targetY = yB.divide(16).negate(); 290 } else if (agingB >= MAXIMAL_AGING) { 291 // we keep updating the low bracket, try to compensate this 292 targetY = yA.divide(16).negate(); 293 } else { 294 // bracketing is balanced, try to find the root itself 295 targetY = zero; 296 } 297 298 // make a few attempts to guess a root, 299 Dfp nextX; 300 int start = 0; 301 int end = nbPoints; 302 do { 303 304 // guess a value for current target, using inverse polynomial interpolation 305 System.arraycopy(x, start, tmpX, start, end - start); 306 nextX = guessX(targetY, tmpX, y, start, end); 307 308 if (!(nextX.greaterThan(xA) && nextX.lessThan(xB))) { 309 // the guessed root is not strictly inside of the tightest bracketing interval 310 311 // the guessed root is either not strictly inside the interval or it 312 // is a NaN (which occurs when some sampling points share the same y) 313 // we try again with a lower interpolation order 314 if (signChangeIndex - start >= end - signChangeIndex) { 315 // we have more points before the sign change, drop the lowest point 316 ++start; 317 } else { 318 // we have more points after sign change, drop the highest point 319 --end; 320 } 321 322 // we need to do one more attempt 323 nextX = nan; 324 325 } 326 327 } while (nextX.isNaN() && (end - start > 1)); 328 329 if (nextX.isNaN()) { 330 // fall back to bisection 331 nextX = xA.add(xB.subtract(xA).divide(2)); 332 start = signChangeIndex - 1; 333 end = signChangeIndex; 334 } 335 336 // evaluate the function at the guessed root 337 evaluations.incrementCount(); 338 final Dfp nextY = f.value(nextX); 339 if (nextY.isZero()) { 340 // we have found an exact root, since it is not an approximation 341 // we don't need to bother about the allowed solutions setting 342 return nextX; 343 } 344 345 if ((nbPoints > 2) && (end - start != nbPoints)) { 346 347 // we have been forced to ignore some points to keep bracketing, 348 // they are probably too far from the root, drop them from now on 349 nbPoints = end - start; 350 System.arraycopy(x, start, x, 0, nbPoints); 351 System.arraycopy(y, start, y, 0, nbPoints); 352 signChangeIndex -= start; 353 354 } else if (nbPoints == x.length) { 355 356 // we have to drop one point in order to insert the new one 357 nbPoints--; 358 359 // keep the tightest bracketing interval as centered as possible 360 if (signChangeIndex >= (x.length + 1) / 2) { 361 // we drop the lowest point, we have to shift the arrays and the index 362 System.arraycopy(x, 1, x, 0, nbPoints); 363 System.arraycopy(y, 1, y, 0, nbPoints); 364 --signChangeIndex; 365 } 366 367 } 368 369 // insert the last computed point 370 //(by construction, we know it lies inside the tightest bracketing interval) 371 System.arraycopy(x, signChangeIndex, x, signChangeIndex + 1, nbPoints - signChangeIndex); 372 x[signChangeIndex] = nextX; 373 System.arraycopy(y, signChangeIndex, y, signChangeIndex + 1, nbPoints - signChangeIndex); 374 y[signChangeIndex] = nextY; 375 ++nbPoints; 376 377 // update the bracketing interval 378 if (nextY.multiply(yA).negativeOrNull()) { 379 // the sign change occurs before the inserted point 380 xB = nextX; 381 yB = nextY; 382 absYB = yB.abs(); 383 ++agingA; 384 agingB = 0; 385 } else { 386 // the sign change occurs after the inserted point 387 xA = nextX; 388 yA = nextY; 389 absYA = yA.abs(); 390 agingA = 0; 391 ++agingB; 392 393 // update the sign change index 394 signChangeIndex++; 395 396 } 397 398 } 399 400 } 401 402 /** Guess an x value by n<sup>th</sup> order inverse polynomial interpolation. 403 * <p> 404 * The x value is guessed by evaluating polynomial Q(y) at y = targetY, where Q 405 * is built such that for all considered points (x<sub>i</sub>, y<sub>i</sub>), 406 * Q(y<sub>i</sub>) = x<sub>i</sub>. 407 * </p> 408 * @param targetY target value for y 409 * @param x reference points abscissas for interpolation, 410 * note that this array <em>is</em> modified during computation 411 * @param y reference points ordinates for interpolation 412 * @param start start index of the points to consider (inclusive) 413 * @param end end index of the points to consider (exclusive) 414 * @return guessed root (will be a NaN if two points share the same y) 415 */ 416 private Dfp guessX(final Dfp targetY, final Dfp[] x, final Dfp[] y, 417 final int start, final int end) { 418 419 // compute Q Newton coefficients by divided differences 420 for (int i = start; i < end - 1; ++i) { 421 final int delta = i + 1 - start; 422 for (int j = end - 1; j > i; --j) { 423 x[j] = x[j].subtract(x[j-1]).divide(y[j].subtract(y[j - delta])); 424 } 425 } 426 427 // evaluate Q(targetY) 428 Dfp x0 = targetY.getZero(); 429 for (int j = end - 1; j >= start; --j) { 430 x0 = x[j].add(x0.multiply(targetY.subtract(y[j]))); 431 } 432 433 return x0; 434 435 } 436 437 }