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    }