001/* 002* Minimization -- A collection of static methods for finding the minima or maxima of functions. 003* 004* Copyright (C) 2006-2025, by Joseph A. Huwaldt. All rights reserved. 005* 006* This library is free software; you can redistribute it and/or 007* modify it under the terms of the GNU Lesser General Public 008* License as published by the Free Software Foundation; either 009* version 2.1 of the License, or (at your option) any later version. 010* 011* This library is distributed in the hope that it will be useful, 012* but WITHOUT ANY WARRANTY; without even the implied warranty of 013* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 014* Lesser General Public License for more details. 015* 016* You should have received a copy of the GNU Lesser General Public License 017* along with this program; if not, write to the Free Software 018* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 019* Or visit: http://www.gnu.org/licenses/lgpl.html 020*/ 021package jahuwaldt.tools.math; 022 023import static java.lang.Math.*; 024 025/** 026 * A collection of static routines to find the minima or maxima of functions or sets of 027 * functions. 028 * 029 * <p> Modified by: Joseph A. Huwaldt </p> 030 * 031 * @author Joseph A. Huwaldt Date: July 22, 2006 032 * @version February 23, 2025 033 */ 034public class Minimization { 035 036 /** 037 * Machine floating point precision. 038 */ 039 private static final double ZEPS = 1e-10; 040 041 /** 042 * Maximum number of allowed iterations. 043 */ 044 private static final int MAXIT = 1000; 045 046 /** 047 * The default ratio by which successive intervals are magnified. 048 */ 049 private static final double GOLD = 1.618034; 050 051 /** 052 * The golden ratio. 053 */ 054 private static final double CGOLD = 0.3819660; 055 056 private static final double GLIMIT = 100; 057 private static final double TINY = 1.0e-20; 058 private static final double TINY2 = 1.0e-25; 059 private static final double TOL = 1.0e-8; 060 061 //----------------------------------------------------------------------------------- 062 /** 063 * <p> 064 * Method that isolates the minimum of a 1D function to a fractional precision of 065 * about "tol". A version of Brent's method is used with derivative information if it 066 * is available or a version without if it is not.</p> 067 * 068 * @param eval An evaluatable 1D function that returns the function value at x and 069 * optionally returns the function derivative value (d(fx)/dx) at x. It 070 * is guaranteed that "function()" will always be called before 071 * "derivative()". 072 * @param x1 The lower bracket surrounding the minimum (assumed that minimum lies 073 * between x1 and x2). 074 * @param x2 The upper bracket surrounding the root (assumed that minimum lies 075 * between x1 and x2). 076 * @param tol The root will be refined until it's accuracy is better than this. 077 * This should generally not be smaller than Math.sqrt(EPS)! 078 * @param output An optional 2-element array that will be filled in with the abscissa 079 * (x) and function minimum ordinate ( f(x) ) on output. 080 * output[xmin,f(xmin)]. Pass "null" if this is not required. 081 * @return The abscissa (x) of the minimum value of f(x). 082 * @exception RootException Unable to find a minimum of the function. 083 */ 084 public static double find(Evaluatable1D eval, double x1, double x2, double tol, double[] output) throws RootException { 085 086 // First find an initial triplet that brackets the minimum. 087 double[] triplet = new double[3]; 088 bracket1D(x1, x2, triplet, null, eval); 089 090 // Now find the minimum. 091 return find(eval, triplet[0], triplet[1], triplet[2], tol, output); 092 } 093 094 /** 095 * <p> 096 * Method that isolates the minimum of a 1D function to a fractional precision of 097 * about "tol". A version of Brent's method is used with derivative information if it 098 * is available or a version without if it is not.</p> 099 * 100 * @param eval An evaluatable 1D function that returns the function value at x and 101 * optionally returns the function derivative value (d(fx)/dx) at x. It 102 * is guaranteed that "function()" will always be called before 103 * "derivative()". The return value from the 1st call to "function()" is 104 * ignored. 105 * @param ax The lower bracket known to surround the minimum (assumed that minimum 106 * lies between ax and cx). 107 * @param bx Abscissa point that lies between ax and cx where f(bx) is less than 108 * both f(ax) and f(cx). 109 * @param cx The upper bracket known to surround the root (assumed that minimum 110 * lies between ax and cx). 111 * @param tol The root will be refined until it's accuracy is better than this. 112 * This should generally not be smaller than Math.sqrt(EPS)! 113 * @param output An optional 2-element array that will be filled in with the abscissa 114 * (x) and function minimum ordinate ( f(x) ) on output. 115 * output[xmin,f(xmin)]. Pass "null" if this is not required. 116 * @return The abscissa (x) of the minimum value of f(x). 117 * @exception RootException Unable to find a minimum of the function. 118 */ 119 public static double find(Evaluatable1D eval, double ax, double bx, double cx, double tol, double[] output) throws RootException { 120 if (tol < MathTools.EPS) 121 tol = MathTools.EPS; 122 eval.function(bx); 123 double d = eval.derivative(bx); 124 if (Double.isNaN(d)) 125 return brent(eval, ax, bx, cx, tol, output); 126 return dbrent(eval, ax, bx, cx, tol, output); 127 } 128 129 /** 130 * Method that searches in the downhill direction (defined by the function as 131 * evaluated at the initial points) and returns new points (in triple) that bracket a 132 * minimum of the function. 133 * 134 * @param ax The lower bracket surrounding the minimum (assumed that minimum lies 135 * between ax and bx). 136 * @param bx The upper bracket surrounding the root (assumed that minimum lies 137 * between ax and bx). 138 * @param triple A 3-element array that will be filled in with the abcissa of the 139 * three points that bracket the minium triple[ax,bx,cx]. 140 * @param values A 3-element array that contains the function values that correspond 141 * with the points ax, bx and cx in triple. Null may be passed for this 142 * if the function values are not required. 143 * @param eval An evaluatable 1D function that returns the function value at x and 144 * optionally returns the function derivative value (d(fx)/dx) at x. It 145 * is guaranteed that "function()" will always be called before 146 * "derivative()". 147 * @throws RootException if there is any problem finding the bracket. 148 */ 149 public static void bracket1D(double ax, double bx, double[] triple, double[] values, Evaluatable1D eval) 150 throws RootException { 151 double fa = eval.function(ax); 152 double fb = eval.function(bx); 153 if (fb > fa) { 154 // Switch roles of a and b so that we can go downhill in the direction from a to b. 155 double dum = ax; 156 ax = bx; 157 bx = dum; 158 dum = fb; 159 fb = fa; 160 fa = dum; 161 } 162 163 double cx = bx + GOLD * (bx - ax); 164 double fc = eval.function(cx); 165 double fu; 166 while (fb > fc) { 167 // Keep returning here until we bracket. 168 double r = (bx - ax) * (fb - fc); // Compute u by parabolic interpolation from a,b,c. 169 double q = (bx - cx) * (fb - fa); 170 double u = bx - ((bx - cx) * q - (bx - ax) * r) 171 / (2.0 * MathTools.sign(max(abs(q - r), TINY), q - r)); 172 double ulim = bx + GLIMIT * (cx - bx); 173 174 if ((bx - u) * (u - cx) > 0.0) { 175 // Parabolic u is between b and c. 176 fu = eval.function(u); 177 if (fu < fc) { 178 // Got a minimum between b and c. 179 ax = bx; 180 bx = u; 181 fa = fb; 182 fb = fu; 183 triple[0] = ax; 184 triple[1] = bx; 185 triple[2] = cx; 186 if (values != null) { 187 values[0] = fa; 188 values[1] = fb; 189 values[2] = fc; 190 } 191 return; 192 193 } else if (fu > fb) { 194 // Got a minimum between a and u. 195 cx = u; 196 fc = fu; 197 triple[0] = ax; 198 triple[1] = bx; 199 triple[2] = cx; 200 if (values != null) { 201 values[0] = fa; 202 values[1] = fb; 203 values[2] = fc; 204 } 205 return; 206 } 207 u = cx + GOLD * (cx - bx); // Parabolic fit was no use. Use default magnification. 208 fu = eval.function(u); 209 210 } else if ((cx - u) * (u - ulim) > 0.0) { 211 // Parabolic fit is between c and it's allowed limit. 212 fu = eval.function(u); 213 if (fu < fc) { 214 bx = cx; 215 cx = u; 216 u = u + GOLD * (u - bx); 217 fb = fc; 218 fc = fu; 219 fu = eval.function(u); 220 } 221 222 } else if ((u - ulim) * (ulim - cx) >= 0.0) { 223 // Limit parabolic u to maximum allowed value. 224 u = ulim; 225 fu = eval.function(u); 226 227 } else { 228 // Reject parabolic u, use default magnification. 229 u = cx + GOLD * (cx - bx); 230 fu = eval.function(u); 231 } 232 233 // Eliminate oldest point and continue. 234 ax = bx; 235 bx = cx; 236 cx = u; 237 fa = fb; 238 fb = fc; 239 fc = fu; 240 } 241 242 triple[0] = ax; 243 triple[1] = bx; 244 triple[2] = cx; 245 if (values != null) { 246 values[0] = fa; 247 values[1] = fb; 248 values[2] = fc; 249 } 250 251 } 252 253 /** 254 * <p> 255 * Method that isolates the minimum of a 1D function to a fractional precision of 256 * about "tol" using Brent's method.</p> 257 * 258 * <p> 259 * Ported from brent() code found in 260 * _Numerical_Recipes_in_C:_The_Art_of_Scientific_Computing_, 2nd Edition, 1992. </p> 261 * 262 * @param eval An evaluatable 1D function that returns the function value at x and 263 * optionally returns the function derivative value (d(fx)/dx) at x. It 264 * is guaranteed that "function()" will always be called before 265 * "derivative()". 266 * @param ax The lower bracket surrounding the minimum (assumed that minimum lies 267 * between ax and cx). 268 * @param bx Abscissa point that lies between ax and cx where f(bx) is less than 269 * both f(ax) and f(cx). 270 * @param cx The upper bracket surrounding the root (assumed that minimum lies 271 * between ax and cx). 272 * @param tol The minimum will be refined until it's accuracy is better than this. 273 * This should generally not be smaller than Math.sqrt(EPS)! 274 * @param output An optional 2-element array that will be filled in with the abscissa 275 * (x) and function minimum ordinate ( f(x) ) on output. 276 * output[xmin,f(xmin)]. Pass <code>null</code> if this is not required. 277 * @return The abscissa (x) of the minimum value of f(x). 278 * @exception RootException Unable to find a minimum of the function. 279 */ 280 private static double brent(Evaluatable1D eval, double ax, double bx, double cx, double tol, double[] output) throws RootException { 281 282 double e = 0; // This is the distance moved on the step before last. 283 284 // a and b must be in ascending order, but the input abscissas need not be. 285 double a = (ax < cx ? ax : cx); 286 double b = (ax > cx ? ax : cx); 287 double x = bx, w = bx, v = bx; 288 double fx = eval.function(x); 289 double fw = fx, fv = fx; 290 double u, fu, d = 0; 291 292 // Main program loop. 293 for (int iter = 0; iter < MAXIT; ++iter) { 294 double xm = 0.5 * (a + b); 295 double tol1 = tol * abs(x) + ZEPS; 296 double tol2 = 2.0 * tol1; 297 298 // Test for done here. 299 if (abs(x - xm) <= (tol2 - 0.5 * (b - a))) { 300 if (output != null) { 301 output[0] = x; 302 output[1] = fx; 303 } 304 return x; 305 } 306 307 if (abs(e) > tol1) { 308 // Construct a trial parabolic fit. 309 double r = (x - w) * (fx - fv); 310 double q = (x - v) * (fx - fw); 311 double p = (x - v) * q - (x - w) * r; 312 q = 2.0 * (q - r); 313 if (q > 0.0) 314 p = -p; 315 q = abs(q); 316 double etemp = e; 317 e = d; 318 if (abs(p) >= abs(0.5 * q * etemp) || p <= q * (a - x) || p >= q * (b - x)) { 319 // Take a golden ratio step. 320 e = (x >= xm ? a - x : b - x); 321 d = CGOLD * e; 322 323 } else { 324 // Take a parabolic step. 325 d = p / q; 326 u = x + d; 327 if (u - a < tol2 || b - u < tol2) 328 d = MathTools.sign(tol1, xm - x); 329 } 330 331 } else { 332 // Take a golden ratio step. 333 e = (x >= xm ? a - x : b - x); 334 d = CGOLD * e; 335 } 336 337 u = (abs(d) >= tol1 ? x + d : x + MathTools.sign(tol1, d)); 338 fu = eval.function(u); // This is the one function eval per iteration. 339 340 if (fu <= fx) { 341 if (u >= x) 342 a = x; 343 else 344 b = x; 345 v = w; 346 w = x; 347 x = u; 348 fv = fw; 349 fw = fx; 350 fx = fu; 351 352 } else { 353 if (u < x) 354 a = u; 355 else 356 b = u; 357 if (fu <= fw || w == x) { 358 v = w; 359 w = u; 360 fv = fw; 361 fw = fu; 362 363 } else if (fu <= fv || v == x || v == w) { 364 v = u; 365 fv = fu; 366 } 367 } 368 369 } // Next iter 370 371 // Max iterations exceeded. 372 throw new RootException("Maximum number of iterations exceeded"); 373 374 } 375 376 /** 377 * <p> 378 * Method that isolates the minimum of a 1D function to a fractional precision of 379 * about "tol" using Brent's method plus derivatives. </p> 380 * 381 * <p> 382 * Ported from dbrent() code found in 383 * _Numerical_Recipes_in_C:_The_Art_of_Scientific_Computing_, 2nd Edition, 1992, pg 384 * 406. </p> 385 * 386 * @param eval An evaluatable 1D function that returns the function value at x and 387 * returns the function derivative value (d(fx)/dx) at x. It is 388 * guaranteed that "function()" will always be called before 389 * "derivative()". 390 * @param ax The lower bracket surrounding the minimum (assumed that minimum lies 391 * between ax and cx). 392 * @param bx Abscissa point that lies between ax and cx where f(bx) is less than 393 * both f(ax) and f(cx). 394 * @param cx The upper bracket surrounding the root (assumed that minimum lies 395 * between ax and cx). 396 * @param tol The minimum will be refined until it's accuracy is better than this. 397 * This should generally not be smaller than Math.sqrt(EPS)! 398 * @param output An optional 2-element array that will be filled in with the abscissa 399 * (x) and function minimum ordinate ( f(x) ) on output. 400 * output[xmin,f(xmin)]. Pass <code>null</code> if this is not required. 401 * @return The abscissa (x) of the minimum value of f(x). 402 * @exception RootException Unable to find a minimum of the function. 403 */ 404 private static double dbrent(Evaluatable1D eval, double ax, double bx, double cx, double tol, double[] output) 405 throws RootException { 406 double a = (ax < cx ? ax : cx); 407 double b = (ax > cx ? ax : cx); 408 double x = bx; 409 double w = bx; 410 double v = bx; 411 double fw = eval.function(x); 412 double fv = fw; 413 double fx = fw; 414 double fu; 415 double dw = eval.derivative(x); 416 double dv = dw; 417 double dx = dw; 418 double d = 0.0; 419 double e = 0.0; 420 double u; 421 422 for (int iter = 0; iter < MAXIT; iter++) { 423 double xm = 0.5 * (a + b); 424 double tol1 = tol * abs(x) + ZEPS; 425 double tol2 = 2.0 * tol1; 426 if (abs(x - xm) <= (tol2 - 0.5 * (b - a))) { 427 if (output != null) { 428 output[0] = x; 429 output[1] = fx; 430 } 431 return x; 432 } 433 if (abs(e) > tol1) { 434 double d1 = 2.0 * (b - a); 435 double d2 = d1; 436 if (dw != dx) 437 d1 = (w - x) * dx / (dx - dw); 438 if (dv != dx) 439 d2 = (v - x) * dx / (dx - dv); 440 double u1 = x + d1; 441 double u2 = x + d2; 442 boolean ok1 = (a - u1) * (u1 - b) > 0.0 && dx * d1 <= 0.0; 443 boolean ok2 = (a - u2) * (u2 - b) > 0.0 && dx * d2 <= 0.0; 444 double olde = e; 445 e = d; 446 if (ok1 || ok2) { 447 if (ok1 && ok2) 448 d = (abs(d1) < abs(d2) ? d1 : d2); 449 else if (ok1) 450 d = d1; 451 else 452 d = d2; 453 if (abs(d) <= abs(0.5 * olde)) { 454 u = x + d; 455 if (u - a < tol2 || b - u < tol2) 456 d = MathTools.sign(tol1, xm - x); 457 } else { 458 d = 0.5 * (e = (dx >= 0.0 ? a - x : b - x)); 459 } 460 } else { 461 d = 0.5 * (e = (dx >= 0.0 ? a - x : b - x)); 462 } 463 } else { 464 d = 0.5 * (e = (dx >= 0.0 ? a - x : b - x)); 465 } 466 if (abs(d) >= tol1) { 467 u = x + d; 468 fu = eval.function(u); 469 } else { 470 u = x + MathTools.sign(tol1, d); 471 fu = eval.function(u); 472 if (fu > fx) { 473 if (output != null) { 474 output[0] = x; 475 output[1] = fx; 476 } 477 return x; 478 } 479 } 480 double du = eval.derivative(u); 481 if (fu <= fx) { 482 if (u >= x) 483 a = x; 484 else 485 b = x; 486 v = w; 487 fv = fw; 488 dv = dw; 489 w = x; 490 fw = fx; 491 dw = dx; 492 x = u; 493 fx = fu; 494 dx = du; 495 } else { 496 if (u < x) 497 a = u; 498 else 499 b = u; 500 if (fu <= fw || w == x) { 501 v = w; 502 fv = fw; 503 dv = dw; 504 w = u; 505 fw = fu; 506 dw = du; 507 } else if (fu < fv || v == x || v == w) { 508 v = u; 509 fv = fu; 510 dv = du; 511 } 512 } 513 } 514 515 // Max iterations exceeded. 516 throw new RootException("Maximum number of iterations exceeded"); 517 } 518 519 /** 520 * <p> 521 * Method that isolates the minimum of an n-Dimensional scalar function to a 522 * fractional precision of about "tol". If no derivatives are available then Powell's 523 * method is used. If derivatives are available, then a Polak-Reibiere minimization is 524 * used.</p> 525 * 526 * @param func An evaluatable n-D scalar function that returns the function value at a 527 * point, p[1..n], and optionally returns the function derivatives 528 * (d(fx[1..n])/dx) at p[1..n]. It is guaranteed that "function()" will 529 * always be called before "derivatives()". 530 * @param p On input, this is the initial guess at the minimum point. On output, 531 * this is filled in with the values that minimize the function. 532 * @param n The number of variables in the array p. 533 * @param tol The convergence tolerance on the function value. 534 * @return The minimum value of f(p[1..n]). 535 * @exception RootException if unable to find a minimum of the function. 536 */ 537 public static double findND(ScalarFunctionND func, double[] p, int n, double tol) 538 throws RootException { 539 if (tol < MathTools.EPS) 540 tol = MathTools.EPS; 541 double[] df = new double[n]; 542 double f = func.function(p); 543 boolean hasDerivatives = func.derivatives(p, df); 544 if (hasDerivatives) 545 // We have derivatives, so use Polak-Reibiere method. 546 return frprmn(p, n, tol, func); 547 //df = null; 548 549 // No derivatives, so use Powell's method. 550 double[][] xi = new double[n][n]; 551 for (int i = 0; i < n; i++) 552 for (int j = 0; j < n; j++) 553 xi[i][j] = (i == j ? 1.0 : 0.0); 554 555 return powell(p, n, xi, tol, func); 556 } 557 558 /** 559 * Given a starting point, p[1..n], a Powell's method minimization is performed on a 560 * function, func. 561 * 562 * Reference: Numerical Recipes in C, 2nd Edition, pg 417. 563 * 564 * @param p On input, this is the initial guess at the minimum point. On output, 565 * this is filled in with the values that minimize the function. 566 * @param n The number of variables in the array p. 567 * @param xi An existing nxn matrix who's columns contain the initial set of 568 * directions (usually the n unit vectors). 569 * @param ftol The convergence tolerance on the function value. 570 * @param func The function being minimized and it's derivative. 571 * @return The minimum function value. 572 */ 573 private static double powell(double[] p, int n, double[][] xi, double ftol, ScalarFunctionND func) 574 throws RootException { 575 double fptt; 576 double[] pt = new double[n]; 577 double[] ptt = new double[n]; 578 double[] xit = new double[n]; 579 double fret = func.function(p); 580 System.arraycopy(p, 0, pt, 0, n); // for (int j=0;j < n;j++) pt[j]=p[j]; // Save the initial point. 581 for (int iter = 0;; ++iter) { 582 double fp = fret; 583 int ibig = 0; 584 double del = 0.0; // Will be the biggest function decrease. 585 586 // Loop over all the directions in the set. 587 for (int i = 0; i < n; i++) { 588 for (int j = 0; j < n; j++) 589 xit[j] = xi[j][i]; // Copy the direction. 590 fptt = fret; 591 fret = linmin(n, p, xit, func); // Minimize along direction. 592 593 // Record if the largest decrease so far. 594 if (abs(fptt - fret) > del) { 595 del = abs(fptt - fret); 596 ibig = i; 597 } 598 } 599 600 // Termination criteria. 601 double ferr = 2.0 * abs(fp - fret); 602 double fcrit = ftol * (abs(fp) + abs(fret)) + TINY2; 603 if (ferr <= fcrit) { 604 return fret; 605 } 606 if (iter == MAXIT) 607 throw new RootException("Exceeded maximum iterations."); 608 609 // Construct extrapolated point and average direction moved and save old starting point. 610 for (int j = 0; j < n; j++) { 611 double pj = p[j]; 612 double ptj = pt[j]; 613 ptt[j] = 2.0 * pj - ptj; 614 xit[j] = pj - ptj; 615 pt[j] = pj; 616 } 617 618 fptt = func.function(ptt); // Function value at extrapolated point. 619 if (fptt < fp) { 620 double tmp1 = fp - fret - del; 621 double tmp2 = fp - fptt; 622 double t = 2.0 * (fp - 2.0 * fret + fptt) * tmp1 * tmp1 - del * tmp2 * tmp2; 623 if (t < 0.0) { 624 // Move to the minimum of the new direction, and save the new direciton. 625 fret = linmin(n, p, xit, func); 626 for (int j = 0; j < n; j++) { 627 xi[j][ibig] = xi[j][n - 1]; 628 xi[j][n - 1] = xit[j]; 629 } 630 } 631 } 632 } // Next iteration 633 } 634 635 /** 636 * Given a starting point, p[1..n], a Polak-Reibiere minimization is performed on a 637 * function, func, using it's gradient. 638 * 639 * Reference: Numerical Recipes in C, 2nd Edition, pg 423. 640 * 641 * @param p On input, this is the initial guess at the minimum point. On output, 642 * this is filled in with the variables that minimize the function. 643 * @param n The number of variables in the array p. 644 * @param ftol The convergence tolerance on the function value. 645 * @param func The function being minimized and it's derivative. 646 * @return The minimum function value. 647 */ 648 private static double frprmn(double[] p, int n, double ftol, ScalarFunctionND func) 649 throws RootException { 650 651 double[] g = new double[n]; 652 double[] h = new double[n]; 653 double[] xi = new double[n]; 654 655 // Initializations 656 double fp = func.function(p); 657 func.derivatives(p, xi); 658 659 for (int j = 0; j < n; j++) { 660 g[j] = -xi[j]; 661 } 662 System.arraycopy(g, 0, h, 0, n); 663 System.arraycopy(g, 0, xi, 0, n); 664 665 // Loop over the iterations. 666 for (int its = 0; its < MAXIT; its++) { 667 double fret = linmin(n, p, xi, func); 668 if (2.0 * abs(fret - fp) <= ftol * (abs(fret) + abs(fp) + TINY)) 669 return fret; // Normal return. 670 fp = fret; 671 func.derivatives(p, xi); 672 double dgg = 0.0; 673 double gg = 0.0; 674 for (int j = 0; j < n; j++) { 675 double gj = g[j]; 676 gg += gj * gj; 677 double xij = xi[j]; 678 // dgg += xij*xij]; // Statement for Fletcher-Reeves. 679 dgg += (xij + gj) * xij; // Statement for Polak-Ribiere. 680 } 681 if (gg == 0.0) // Unlikely, if gradient is exactly 0, then we are already done. 682 return fret; 683 684 double gam = dgg / gg; 685 for (int j = 0; j < n; j++) { 686 g[j] = -xi[j]; 687 xi[j] = h[j] = g[j] + gam * h[j]; 688 } 689 } 690 throw new RootException("Too many iterations required for minimization."); 691 } 692 693 /** 694 * Given an n-dimensional point, p[1..n], and an n-dimensional direction, xi[1..n], 695 * moves and resets p to where the function, func(p), takes on a minimum along the 696 * direction xi from p, and replaces xi by the actual vector displacement that p was 697 * moved. Also returns the function value at p. 698 * 699 * Reference: Numerical Recipes in C, 2nd Edition, pg 419. 700 */ 701 private static double linmin(int n, double[] p, double[] xi, ScalarFunctionND func) 702 throws RootException { 703 704 Evaluatable1D f1dim = new F1DimFunction(p, xi, n, func); 705 double[] triple = new double[3]; 706 double[] values = new double[3]; 707 bracket1D(0, 1, triple, values, f1dim); 708 double ax = triple[0]; 709 double xx = triple[1]; 710 double bx = triple[2]; 711 712 double[] outputs = new double[2]; 713 double xmin = find(f1dim, ax, xx, bx, TOL, outputs); 714 double fret = outputs[1]; 715 for (int j = 0; j < n; j++) { 716 xi[j] *= xmin; 717 p[j] += xi[j]; 718 } 719 720 return fret; 721 } 722 723 /** 724 * This is the value of the function, nrfunc, along the line going through the point p 725 * in the direction of xi. This is used by linmin(). 726 */ 727 private static class F1DimFunction extends AbstractEvaluatable1D { 728 729 private final double[] _pcom; 730 private final double[] _xicom; 731 private final int _ncom; 732 private final double[] _xt; 733 private final double[] _df; 734 private final ScalarFunctionND _nrfunc; 735 736 public F1DimFunction(double[] p, double[] xi, int n, ScalarFunctionND nrfunc) { 737 _pcom = p; 738 _xicom = xi; 739 _ncom = n; 740 _xt = new double[n]; 741 _df = new double[n]; 742 _nrfunc = nrfunc; 743 } 744 745 @Override 746 public double function(double x) throws RootException { 747 for (int j = 0; j < _ncom; j++) 748 _xt[j] = _pcom[j] + x * _xicom[j]; 749 return _nrfunc.function(_xt); 750 } 751 752 @Override 753 public double derivative(double x) throws RootException { 754 boolean hasDerivatives = _nrfunc.derivatives(_xt, _df); 755 if (hasDerivatives) { 756 double df1 = 0; 757 for (int j = 0; j < _ncom; j++) 758 df1 += _df[j] * _xicom[j]; 759 return df1; 760 } 761 return Double.NaN; 762 } 763 } 764 765 /** 766 * Used to test out the methods in this class. 767 * 768 * @param args the command-line arguments 769 */ 770 public static void main(String args[]) { 771 772 System.out.println(); 773 System.out.println("Testing Minimization..."); 774 775 try { 776 // Find the minimum of f(x) = x^2 - 3*x + 2. 777 778 // First create an instance of our evaluatable function. 779 Evaluatable1D function = new DemoMinFunction(); 780 781 // Create some temporary storage. 782 double[] output = new double[2]; 783 784 // Then call the minimizer with brackets on the minima and a tolerance. 785 // The brackets can be used to isolate which minima you want (for instance). 786 double xmin = find(function, -10, 10, 0.0001, output); 787 788 System.out.println(" Minimum of f(x) = x^2 - 3*x + 2 is x = " + (float) xmin 789 + ", f(x) = " + (float) output[1]); 790 791 } catch (Exception e) { 792 e.printStackTrace(); 793 } 794 795 } 796 797 /** 798 * <p> 799 * A simple demonstration class that shows how to use the Evaluatable1D class to find 800 * the minimum of a 1D equation.</p> 801 */ 802 private static class DemoMinFunction extends AbstractEvaluatable1D { 803 804 /** 805 * User supplied method that calculates the function y = fn(x). This demonstration 806 * returns the value of y = x^2 - 3*x + 2. 807 * 808 * @param x Independent parameter to the function. 809 * @return The x^2 - 3*x + 2 evaluated at x. 810 */ 811 @Override 812 public double function(double x) { 813 return (x * x - 3. * x + 2.); 814 } 815 816 /** 817 * Calculates the derivative of the function <code>dy/dx = d fn(x)/dx</code>. This 818 * demonstration returns <code>dy/dx = 2*x - 3</code>. 819 * 820 * @param x Independent parameter to the function. 821 * @return The 2*x - 3 evaluated at x. 822 */ 823 @Override 824 public double derivative(double x) { 825 return (2. * x - 3.); 826 } 827 828 } 829 830}