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.optim.univariate;

  19. import org.apache.commons.math3.util.FastMath;
  20. import org.apache.commons.math3.util.Incrementor;
  21. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  22. import org.apache.commons.math3.exception.TooManyEvaluationsException;
  23. import org.apache.commons.math3.exception.MaxCountExceededException;
  24. import org.apache.commons.math3.analysis.UnivariateFunction;
  25. import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;

  27. /**
  28.  * Provide an interval that brackets a local optimum of a function.
  29.  * This code is based on a Python implementation (from <em>SciPy</em>,
  30.  * module {@code optimize.py} v0.5).
  31.  *
  32.  * @since 2.2
  33.  */
  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;

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

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

  98.         this.growLimit = growLimit;
  99.         evaluations.setMaximalCount(maxEvaluations);
  100.     }

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

  119.         double fA = eval(func, xA);
  120.         double fB = eval(func, xB);
  121.         if (isMinim ?
  122.             fA < fB :
  123.             fA > fB) {

  125.             double tmp = xA;
  126.             xA = xB;
  127.             xB = tmp;

  129.             tmp = fA;
  130.             fA = fB;
  131.             fB = tmp;
  132.         }

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

  137.         while (isMinim ? fC < fB : fC > fB) {
  138.             double tmp1 = (xB - xA) * (fB - fC);
  139.             double tmp2 = (xB - xC) * (fB - fA);

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

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

  147.             double fW;
  148.             if ((w - xC) * (xB - w) > 0) {
  149.                 fW = eval(func, w);
  150.                 if (isMinim ?
  151.                     fW < fC :
  152.                     fW > fC) {
  153.                     xA = xB;
  154.                     xB = w;
  155.                     fA = fB;
  156.                     fB = fW;
  157.                     break;
  158.                 } else if (isMinim ?
  159.                            fW > fB :
  160.                            fW < fB) {
  161.                     xC = w;
  162.                     fC = fW;
  163.                     break;
  164.                 }
  165.                 w = xC + GOLD * (xC - xB);
  166.                 fW = eval(func, w);
  167.             } else if ((w - wLim) * (wLim - xC) >= 0) {
  168.                 w = wLim;
  169.                 fW = eval(func, w);
  170.             } else if ((w - wLim) * (xC - w) > 0) {
  171.                 fW = eval(func, w);
  172.                 if (isMinim ?
  173.                     fW < fC :
  174.                     fW > fC) {
  175.                     xB = xC;
  176.                     xC = w;
  177.                     w = xC + GOLD * (xC - xB);
  178.                     fB = fC;
  179.                     fC =fW;
  180.                     fW = eval(func, w);
  181.                 }
  182.             } else {
  183.                 w = xC + GOLD * (xC - xB);
  184.                 fW = eval(func, w);
  185.             }

  187.             xA = xB;
  188.             fA = fB;
  189.             xB = xC;
  190.             fB = fC;
  191.             xC = w;
  192.             fC = fW;
  193.         }

  195.         lo = xA;
  196.         fLo = fA;
  197.         mid = xB;
  198.         fMid = fB;
  199.         hi = xC;
  200.         fHi = fC;

  202.         if (lo > hi) {
  203.             double tmp = lo;
  204.             lo = hi;
  205.             hi = tmp;

  207.             tmp = fLo;
  208.             fLo = fHi;
  209.             fHi = tmp;
  210.         }
  211.     }

  213.     /**
  214.      * @return the number of evalutations.
  215.      */
  216.     public int getMaxEvaluations() {
  217.         return evaluations.getMaximalCount();
  218.     }

  220.     /**
  221.      * @return the number of evalutations.
  222.      */
  223.     public int getEvaluations() {
  224.         return evaluations.getCount();
  225.     }

  227.     /**
  228.      * @return the lower bound of the bracket.
  229.      * @see #getFLo()
  230.      */
  231.     public double getLo() {
  232.         return lo;
  233.     }

  235.     /**
  236.      * Get function value at {@link #getLo()}.
  237.      * @return function value at {@link #getLo()}
  238.      */
  239.     public double getFLo() {
  240.         return fLo;
  241.     }

  243.     /**
  244.      * @return the higher bound of the bracket.
  245.      * @see #getFHi()
  246.      */
  247.     public double getHi() {
  248.         return hi;
  249.     }

  251.     /**
  252.      * Get function value at {@link #getHi()}.
  253.      * @return function value at {@link #getHi()}
  254.      */
  255.     public double getFHi() {
  256.         return fHi;
  257.     }

  259.     /**
  260.      * @return a point in the middle of the bracket.
  261.      * @see #getFMid()
  262.      */
  263.     public double getMid() {
  264.         return mid;
  265.     }

  267.     /**
  268.      * Get function value at {@link #getMid()}.
  269.      * @return function value at {@link #getMid()}
  270.      */
  271.     public double getFMid() {
  272.         return fMid;
  273.     }

  275.     /**
  276.      * @param f Function.
  277.      * @param x Argument.
  278.      * @return {@code f(x)}
  279.      * @throws TooManyEvaluationsException if the maximal number of evaluations is
  280.      * exceeded.
  281.      */
  282.     private double eval(UnivariateFunction f, double x) {
  283.         try {
  284.             evaluations.incrementCount();
  285.         } catch (MaxCountExceededException e) {
  286.             throw new TooManyEvaluationsException(e.getMax());
  287.         }
  288.         return f.value(x);
  289.     }
  290. }
Created with JaCoCo 0.7.2.201409121644