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

  19. import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
  20. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  21. import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
  22. import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException;
  23. import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
  24. import org.apache.commons.math4.core.jdkmath.JdkMath;
  25. import org.apache.commons.math4.legacy.core.IntegerSequence.Incrementor;

  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.      * 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;

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

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

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

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

  123.         double fA = eval.value(xA);
  124.         double fB = eval.value(xB);
  125.         if (isMinim ?
  126.             fA < fB :
  127.             fA > fB) {

  129.             double tmp = xA;
  130.             xA = xB;
  131.             xB = tmp;

  133.             tmp = fA;
  134.             fA = fB;
  135.             fB = tmp;
  136.         }

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

  141.         while (isMinim ? fC < fB : fC > fB) {
  142.             double tmp1 = (xB - xA) * (fB - fC);
  143.             double tmp2 = (xB - xC) * (fB - fA);

  145.             double val = tmp2 - tmp1;
  146.             double denom = JdkMath.abs(val) < EPS_MIN ? 2 * EPS_MIN : 2 * val;

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

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

  191.             xA = xB;
  192.             fA = fB;
  193.             xB = xC;
  194.             fB = fC;
  195.             xC = w;
  196.             fC = fW;
  197.         }

  199.         lo = xA;
  200.         fLo = fA;
  201.         mid = xB;
  202.         fMid = fB;
  203.         hi = xC;
  204.         fHi = fC;

  206.         if (lo > hi) {
  207.             double tmp = lo;
  208.             lo = hi;
  209.             hi = tmp;

  211.             tmp = fLo;
  212.             fLo = fHi;
  213.             fHi = tmp;
  214.         }
  215.     }

  217.     /**
  218.      * @return the number of evaluations.
  219.      */
  220.     public int getMaxEvaluations() {
  221.         return maxEvaluations;
  222.     }

  224.     /**
  225.      * @return the number of evaluations.
  226.      */
  227.     public int getEvaluations() {
  228.         return evaluations;
  229.     }

  231.     /**
  232.      * @return the lower bound of the bracket.
  233.      * @see #getFLo()
  234.      */
  235.     public double getLo() {
  236.         return lo;
  237.     }

  239.     /**
  240.      * Get function value at {@link #getLo()}.
  241.      * @return function value at {@link #getLo()}
  242.      */
  243.     public double getFLo() {
  244.         return fLo;
  245.     }

  247.     /**
  248.      * @return the higher bound of the bracket.
  249.      * @see #getFHi()
  250.      */
  251.     public double getHi() {
  252.         return hi;
  253.     }

  255.     /**
  256.      * Get function value at {@link #getHi()}.
  257.      * @return function value at {@link #getHi()}
  258.      */
  259.     public double getFHi() {
  260.         return fHi;
  261.     }

  263.     /**
  264.      * @return a point in the middle of the bracket.
  265.      * @see #getFMid()
  266.      */
  267.     public double getMid() {
  268.         return mid;
  269.     }

  271.     /**
  272.      * Get function value at {@link #getMid()}.
  273.      * @return function value at {@link #getMid()}
  274.      */
  275.     public double getFMid() {
  276.         return fMid;
  277.     }

  279.     /**
  280.      * Utility for incrementing a counter at each function evaluation.
  281.      */
  282.     private final class FunctionEvaluator {
  283.         /** Function. */
  284.         private final UnivariateFunction func;
  285.         /** Counter. */
  286.         private final Incrementor inc;

  288.         /**
  289.          * @param func Function.
  290.          */
  291.         FunctionEvaluator(UnivariateFunction func) {
  292.             this.func = func;
  293.             inc = Incrementor.create().withMaximalCount(maxEvaluations);
  294.             evaluations = 0;
  295.         }

  297.         /**
  298.          * @param x Argument.
  299.          * @return {@code f(x)}
  300.          * @throws TooManyEvaluationsException if the maximal number of evaluations is
  301.          * exceeded.
  302.          */
  303.         double value(double x) {
  304.             try {
  305.                 inc.increment();
  306.                 evaluations = inc.getCount();
  307.             } catch (MaxCountExceededException e) {
  308.                 throw new TooManyEvaluationsException(e.getMax());
  309.             }

  311.             return func.value(x);
  312.         }
  313.     }
  314. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT