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