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.  *      https://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. /*
  18.  * This is not the original file distributed by the Apache Software Foundation
  19.  * It has been modified by the Hipparchus project
  20.  */
  21. package org.hipparchus.optim.univariate;

  22. import org.hipparchus.analysis.UnivariateFunction;
  23. import org.hipparchus.exception.LocalizedCoreFormats;
  24. import org.hipparchus.exception.MathIllegalArgumentException;
  25. import org.hipparchus.optim.nonlinear.scalar.GoalType;
  26. import org.hipparchus.util.FastMath;
  27. import org.hipparchus.util.Incrementor;

  28. /**
  29.  * Provide an interval that brackets a local optimum of a function.
  30.  * This code is based on a Python implementation (from <em>SciPy</em>,
  31.  * module {@code optimize.py} v0.5).
  32.  *
  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.      * Number of allowed function evaluations.
  47.      */
  48.     private final int maxEvaluations;
  49.     /**
  50.      * Number of function evaluations performed in the last search.
  51.      */
  52.     private int evaluations;
  53.     /**
  54.      * Lower bound of the bracket.
  55.      */
  56.     private double lo;
  57.     /**
  58.      * Higher bound of the bracket.
  59.      */
  60.     private double hi;
  61.     /**
  62.      * Point inside the bracket.
  63.      */
  64.     private double mid;
  65.     /**
  66.      * Function value at {@link #lo}.
  67.      */
  68.     private double fLo;
  69.     /**
  70.      * Function value at {@link #hi}.
  71.      */
  72.     private double fHi;
  73.     /**
  74.      * Function value at {@link #mid}.
  75.      */
  76.     private double fMid;

  77.     /**
  78.      * Constructor with default values {@code 100, 500} (see the
  79.      * {@link #BracketFinder(double,int) other constructor}).
  80.      */
  81.     public BracketFinder() {
  82.         this(100, 500);
  83.     }

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

  101.         this.growLimit = growLimit;
  102.         this.maxEvaluations = maxEvaluations;
  103.     }

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

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

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

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

  132.         double xC = xB + GOLD * (xB - xA);
  133.         double fC = eval.value(xC);

  134.         while (isMinim ? fC < fB : fC > fB) {
  135.             double tmp1 = (xB - xA) * (fB - fC);
  136.             double tmp2 = (xB - xC) * (fB - fA);

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

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

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

  180.             xA = xB;
  181.             fA = fB;
  182.             xB = xC;
  183.             fB = fC;
  184.             xC = w;
  185.             fC = fW;
  186.         }

  187.         lo = xA;
  188.         fLo = fA;
  189.         mid = xB;
  190.         fMid = fB;
  191.         hi = xC;
  192.         fHi = fC;

  193.         if (lo > hi) {
  194.             double tmp = lo;
  195.             lo = hi;
  196.             hi = tmp;

  197.             tmp = fLo;
  198.             fLo = fHi;
  199.             fHi = tmp;
  200.         }
  201.     }

  202.     /** Get maximum number of evaluations.
  203.      * @return the maximum number of evaluations
  204.      */
  205.     public int getMaxEvaluations() {
  206.         return maxEvaluations;
  207.     }

  208.     /** Get number of evaluations.
  209.      * @return the number of evaluations
  210.      */
  211.     public int getEvaluations() {
  212.         return evaluations;
  213.     }

  214.     /** Get lower bound of the bracket.
  215.      * @return the lower bound of the bracket
  216.      * @see #getFLo()
  217.      */
  218.     public double getLo() {
  219.         return lo;
  220.     }

  221.     /** Get function value at {@link #getLo()}.
  222.      * @return function value at {@link #getLo()}
  223.      */
  224.     public double getFLo() {
  225.         return fLo;
  226.     }

  227.     /** Get higher bound of the bracket.
  228.      * @return the higher bound of the bracket
  229.      * @see #getFHi()
  230.      */
  231.     public double getHi() {
  232.         return hi;
  233.     }

  234.     /**
  235.      * Get function value at {@link #getHi()}.
  236.      * @return function value at {@link #getHi()}
  237.      */
  238.     public double getFHi() {
  239.         return fHi;
  240.     }

  241.     /** Get a point in the middle of the bracket.
  242.      * @return a point in the middle of the bracket
  243.      * @see #getFMid()
  244.      */
  245.     public double getMid() {
  246.         return mid;
  247.     }

  248.     /** Get function value at {@link #getMid()}.
  249.      * @return function value at {@link #getMid()}
  250.      */
  251.     public double getFMid() {
  252.         return fMid;
  253.     }

  254.     /**
  255.      * Utility for incrementing a counter at each function evaluation.
  256.      */
  257.     private class FunctionEvaluator {
  258.         /** Function. */
  259.         private final UnivariateFunction func;
  260.         /** Counter. */
  261.         private final Incrementor inc;

  262.         /** Simple constructor.
  263.          * @param func Function.
  264.          */
  265.         FunctionEvaluator(UnivariateFunction func) {
  266.             this.func = func;
  267.             inc = new Incrementor(maxEvaluations);
  268.             evaluations = 0;
  269.         }

  270.         /** Evaluate function.
  271.          * @param x Argument.
  272.          * @return {@code f(x)}
  273.          * @throws org.hipparchus.exception.MathIllegalStateException if the maximal number of evaluations is
  274.          * exceeded.
  275.          */
  276.         double value(double x) {
  277.             inc.increment();
  278.             evaluations = inc.getCount();
  279.             return func.value(x);
  280.         }
  281.     }
  282. }