BracketFinder.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.optimization.univariate;

  18. import org.apache.commons.math3.util.FastMath;
  19. import org.apache.commons.math3.util.Incrementor;
  20. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  21. import org.apache.commons.math3.exception.TooManyEvaluationsException;
  22. import org.apache.commons.math3.exception.MaxCountExceededException;
  23. import org.apache.commons.math3.analysis.UnivariateFunction;
  24. import org.apache.commons.math3.optimization.GoalType;

  25. /**
  26.  * Provide an interval that brackets a local optimum of a function.
  27.  * This code is based on a Python implementation (from <em>SciPy</em>,
  28.  * module {@code optimize.py} v0.5).
  29.  *
  30.  * @deprecated As of 3.1 (to be removed in 4.0).
  31.  * @since 2.2
  32.  */
  33. @Deprecated
  34. public class BracketFinder {
  35.     /** Tolerance to avoid division by zero. */
  36.     private static final double EPS_MIN = 1e-21;
  37.     /**
  38.      * Golden section.
  39.      */
  40.     private static final double GOLD = 1.618034;
  41.     /**
  42.      * Factor for expanding the interval.
  43.      */
  44.     private final double growLimit;
  45.     /**
  46.      * Counter for function evaluations.
  47.      */
  48.     private final Incrementor evaluations = new Incrementor();
  49.     /**
  50.      * Lower bound of the bracket.
  51.      */
  52.     private double lo;
  53.     /**
  54.      * Higher bound of the bracket.
  55.      */
  56.     private double hi;
  57.     /**
  58.      * Point inside the bracket.
  59.      */
  60.     private double mid;
  61.     /**
  62.      * Function value at {@link #lo}.
  63.      */
  64.     private double fLo;
  65.     /**
  66.      * Function value at {@link #hi}.
  67.      */
  68.     private double fHi;
  69.     /**
  70.      * Function value at {@link #mid}.
  71.      */
  72.     private double fMid;

  73.     /**
  74.      * Constructor with default values {@code 100, 50} (see the
  75.      * {@link #BracketFinder(double,int) other constructor}).
  76.      */
  77.     public BracketFinder() {
  78.         this(100, 50);
  79.     }

  80.     /**
  81.      * Create a bracketing interval finder.
  82.      *
  83.      * @param growLimit Expanding factor.
  84.      * @param maxEvaluations Maximum number of evaluations allowed for finding
  85.      * a bracketing interval.
  86.      */
  87.     public BracketFinder(double growLimit,
  88.                          int maxEvaluations) {
  89.         if (growLimit <= 0) {
  90.             throw new NotStrictlyPositiveException(growLimit);
  91.         }
  92.         if (maxEvaluations <= 0) {
  93.             throw new NotStrictlyPositiveException(maxEvaluations);
  94.         }

  95.         this.growLimit = growLimit;
  96.         evaluations.setMaximalCount(maxEvaluations);
  97.     }

  98.     /**
  99.      * Search new points that bracket a local optimum of the function.
  100.      *
  101.      * @param func Function whose optimum should be bracketed.
  102.      * @param goal {@link GoalType Goal type}.
  103.      * @param xA Initial point.
  104.      * @param xB Initial point.
  105.      * @throws TooManyEvaluationsException if the maximum number of evaluations
  106.      * is exceeded.
  107.      */
  108.     public void search(UnivariateFunction func, GoalType goal, double xA, double xB) {
  109.         evaluations.resetCount();
  110.         final boolean isMinim = goal == GoalType.MINIMIZE;

  111.         double fA = eval(func, xA);
  112.         double fB = eval(func, xB);
  113.         if (isMinim ?
  114.             fA < fB :
  115.             fA > fB) {

  116.             double tmp = xA;
  117.             xA = xB;
  118.             xB = tmp;

  119.             tmp = fA;
  120.             fA = fB;
  121.             fB = tmp;
  122.         }

  123.         double xC = xB + GOLD * (xB - xA);
  124.         double fC = eval(func, xC);

  125.         while (isMinim ? fC < fB : fC > fB) {
  126.             double tmp1 = (xB - xA) * (fB - fC);
  127.             double tmp2 = (xB - xC) * (fB - fA);

  128.             double val = tmp2 - tmp1;
  129.             double denom = FastMath.abs(val) < EPS_MIN ? 2 * EPS_MIN : 2 * val;

  130.             double w = xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom;
  131.             double wLim = xB + growLimit * (xC - xB);

  132.             double fW;
  133.             if ((w - xC) * (xB - w) > 0) {
  134.                 fW = eval(func, w);
  135.                 if (isMinim ?
  136.                     fW < fC :
  137.                     fW > fC) {
  138.                     xA = xB;
  139.                     xB = w;
  140.                     fA = fB;
  141.                     fB = fW;
  142.                     break;
  143.                 } else if (isMinim ?
  144.                            fW > fB :
  145.                            fW < fB) {
  146.                     xC = w;
  147.                     fC = fW;
  148.                     break;
  149.                 }
  150.                 w = xC + GOLD * (xC - xB);
  151.                 fW = eval(func, w);
  152.             } else if ((w - wLim) * (wLim - xC) >= 0) {
  153.                 w = wLim;
  154.                 fW = eval(func, w);
  155.             } else if ((w - wLim) * (xC - w) > 0) {
  156.                 fW = eval(func, w);
  157.                 if (isMinim ?
  158.                     fW < fC :
  159.                     fW > fC) {
  160.                     xB = xC;
  161.                     xC = w;
  162.                     w = xC + GOLD * (xC - xB);
  163.                     fB = fC;
  164.                     fC =fW;
  165.                     fW = eval(func, w);
  166.                 }
  167.             } else {
  168.                 w = xC + GOLD * (xC - xB);
  169.                 fW = eval(func, w);
  170.             }

  171.             xA = xB;
  172.             fA = fB;
  173.             xB = xC;
  174.             fB = fC;
  175.             xC = w;
  176.             fC = fW;
  177.         }

  178.         lo = xA;
  179.         fLo = fA;
  180.         mid = xB;
  181.         fMid = fB;
  182.         hi = xC;
  183.         fHi = fC;

  184.         if (lo > hi) {
  185.             double tmp = lo;
  186.             lo = hi;
  187.             hi = tmp;

  188.             tmp = fLo;
  189.             fLo = fHi;
  190.             fHi = tmp;
  191.         }
  192.     }

  193.     /**
  194.      * @return the number of evalutations.
  195.      */
  196.     public int getMaxEvaluations() {
  197.         return evaluations.getMaximalCount();
  198.     }

  199.     /**
  200.      * @return the number of evalutations.
  201.      */
  202.     public int getEvaluations() {
  203.         return evaluations.getCount();
  204.     }

  205.     /**
  206.      * @return the lower bound of the bracket.
  207.      * @see #getFLo()
  208.      */
  209.     public double getLo() {
  210.         return lo;
  211.     }

  212.     /**
  213.      * Get function value at {@link #getLo()}.
  214.      * @return function value at {@link #getLo()}
  215.      */
  216.     public double getFLo() {
  217.         return fLo;
  218.     }

  219.     /**
  220.      * @return the higher bound of the bracket.
  221.      * @see #getFHi()
  222.      */
  223.     public double getHi() {
  224.         return hi;
  225.     }

  226.     /**
  227.      * Get function value at {@link #getHi()}.
  228.      * @return function value at {@link #getHi()}
  229.      */
  230.     public double getFHi() {
  231.         return fHi;
  232.     }

  233.     /**
  234.      * @return a point in the middle of the bracket.
  235.      * @see #getFMid()
  236.      */
  237.     public double getMid() {
  238.         return mid;
  239.     }

  240.     /**
  241.      * Get function value at {@link #getMid()}.
  242.      * @return function value at {@link #getMid()}
  243.      */
  244.     public double getFMid() {
  245.         return fMid;
  246.     }

  247.     /**
  248.      * @param f Function.
  249.      * @param x Argument.
  250.      * @return {@code f(x)}
  251.      * @throws TooManyEvaluationsException if the maximal number of evaluations is
  252.      * exceeded.
  253.      */
  254.     private double eval(UnivariateFunction f, double x) {
  255.         try {
  256.             evaluations.incrementCount();
  257.         } catch (MaxCountExceededException e) {
  258.             throw new TooManyEvaluationsException(e.getMax());
  259.         }
  260.         return f.value(x);
  261.     }
  262. }