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