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}