commons-statistics-inference/src/main/java/org/apache/commons/statistics/inference/BracketFinder.java - commons-statistics - Git at Google

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

 import java.util.function.DoubleUnaryOperator;

 /**
  * Provide an interval that brackets a local minimum of a function.
  * This code is based on a Python implementation (from <em>SciPy</em>,
  * module {@code optimize.py} v0.5).
  *
  * <p>This class has been extracted from {@code o.a.c.math4.optim.univariate}
  * and modified to: remove support for bracketing a maximum; support bounds
  * on the bracket; correct the sign of the denominator when the magnitude is small;
  * and return true/false if there is a minimum strictly inside the bounds.
  *
  * @since 1.1
  */
 class BracketFinder {
     /** Tolerance to avoid division by zero. */
     private static final double EPS_MIN = 1e-21;
     /** Golden section. */
     private static final double GOLD = 1.6180339887498948482;
     /** Factor for expanding the interval. */
     private final double growLimit;
     /**  Number of allowed function evaluations. */
     private final int maxEvaluations;
     /** Number of function evaluations performed in the last search. */
     private int evaluations;
     /** Lower bound of the bracket. */
     private double lo;
     /** Higher bound of the bracket. */
     private double hi;
     /** Point inside the bracket. */
     private double mid;
     /** Function value at {@link #lo}. */
     private double fLo;
     /** Function value at {@link #hi}. */
     private double fHi;
     /** Function value at {@link #mid}. */
     private double fMid;

     /**
      * Constructor with default values {@code 100, 100000} (see the
      * {@link #BracketFinder(double,int) other constructor}).
      */
     BracketFinder() {
         this(100, 100000);
     }

     /**
      * Create a bracketing interval finder.
      *
      * @param growLimit Expanding factor.
      * @param maxEvaluations Maximum number of evaluations allowed for finding
      * a bracketing interval.
      * @throws IllegalArgumentException if the {@code growLimit} or {@code maxEvalutations}
      * are not strictly positive.
      */
     BracketFinder(double growLimit, int maxEvaluations) {
         Arguments.checkStrictlyPositive(growLimit);
         Arguments.checkStrictlyPositive(maxEvaluations);
         this.growLimit = growLimit;
         this.maxEvaluations = maxEvaluations;
     }

     /**
      * Search downhill from the initial points to obtain new points that bracket a local
      * minimum of the function. Note that the initial points do not have to bracket a minimum.
      * An exception is raised if a minimum cannot be found within the configured number
      * of function evaluations.
      *
      * <p>The bracket is limited to the provided bounds if they create a positive interval
      * {@code min < max}. It is possible that the middle of the bracket is at the bounds as
      * the final bracket is {@code f(mid) <= min(f(lo), f(hi))} and {@code lo <= mid <= hi}.
      *
      * <p>No exception is raised if the initial points are not within the bounds; the points
      * are updated to be within the bounds.
      *
      * <p>No exception is raised if the initial points are equal; the bracket will be returned
      * as a single point {@code lo == mid == hi}.
      *
      * @param func Function whose optimum should be bracketed.
      * @param a Initial point.
      * @param b Initial point.
      * @param min Minimum bound of the bracket (inclusive).
      * @param max Maximum bound of the bracket (inclusive).
      * @return true if the mid-point is strictly within the final bracket {@code [lo, hi]};
      * false if there is no local minima.
      * @throws IllegalStateException if the maximum number of evaluations is exceeded.
      */
     boolean search(DoubleUnaryOperator func,
                    double a, double b,
                    double min, double max) {
         evaluations = 0;

         // Limit the range of x
         DoubleUnaryOperator range;
         if (min < max) {
             // Limit: min <= x <= max
             range = x -> {
                 if (x > min) {
                     return x < max ? x : max;
                 }
                 return min;
             };
         } else {
             range = DoubleUnaryOperator.identity();
         }

         double xA = range.applyAsDouble(a);
         double xB = range.applyAsDouble(b);
         double fA = value(func, xA);
         double fB = value(func, xB);
         // Ensure fB <= fA
         if (fA < fB) {
             double tmp = xA;
             xA = xB;
             xB = tmp;
             tmp = fA;
             fA = fB;
             fB = tmp;
         }

         double xC = range.applyAsDouble(xB + GOLD * (xB - xA));
         double fC = value(func, xC);

         // Note: When a [min, max] interval is provided and there is no minima then this
         // loop will terminate when B == C and both are at the min/max bound.
         while (fC < fB) {
             final double tmp1 = (xB - xA) * (fB - fC);
             final double tmp2 = (xB - xC) * (fB - fA);

             final double val = tmp2 - tmp1;
             // limit magnitude of val to a small value
             final double denom = 2 * Math.copySign(Math.max(Math.abs(val), EPS_MIN), val);

             double w = range.applyAsDouble(xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom);
             final double wLim = range.applyAsDouble(xB + growLimit * (xC - xB));

             double fW;
             if ((w - xC) * (xB - w) > 0) {
                 // xB < w < xC
                 fW = value(func, w);
                 if (fW < fC) {
                     // minimum in [xB, xC]
                     xA = xB;
                     xB = w;
                     fA = fB;
                     fB = fW;
                     break;
                 } else if (fW > fB) {
                     // minimum in [xA, w]
                     xC = w;
                     fC = fW;
                     break;
                 }
                 // continue downhill
                 w = range.applyAsDouble(xC + GOLD * (xC - xB));
                 fW = value(func, w);
             } else if ((w - wLim) * (xC - w) > 0) {
                 // xC < w < limit
                 fW = value(func, w);
                 if (fW < fC) {
                     // continue downhill
                     xB = xC;
                     xC = w;
                     w = range.applyAsDouble(xC + GOLD * (xC - xB));
                     fB = fC;
                     fC = fW;
                     fW = value(func, w);
                 }
             } else if ((w - wLim) * (wLim - xC) >= 0) {
                 // xC <= limit <= w
                 w = wLim;
                 fW = value(func, w);
             } else {
                 // possibly w == xC; reject w and take a default step
                 w = range.applyAsDouble(xC + GOLD * (xC - xB));
                 fW = value(func, w);
             }

             xA = xB;
             fA = fB;
             xB = xC;
             fB = fC;
             xC = w;
             fC = fW;
         }

         mid = xB;
         fMid = fB;

         // Store the bracket: lo <= mid <= hi
         if (xC < xA) {
             lo = xC;
             fLo = fC;
             hi = xA;
             fHi = fA;
         } else {
             lo = xA;
             fLo = fA;
             hi = xC;
             fHi = fC;
         }

         return lo < mid && mid < hi;
     }

     /**
      * @return the number of evaluations.
      */
     int getEvaluations() {
         return evaluations;
     }

     /**
      * @return the lower bound of the bracket.
      * @see #getFLo()
      */
     double getLo() {
         return lo;
     }

     /**
      * Get function value at {@link #getLo()}.
      * @return function value at {@link #getLo()}
      */
     double getFLo() {
         return fLo;
     }

     /**
      * @return the higher bound of the bracket.
      * @see #getFHi()
      */
     double getHi() {
         return hi;
     }

     /**
      * Get function value at {@link #getHi()}.
      * @return function value at {@link #getHi()}
      */
     double getFHi() {
         return fHi;
     }

     /**
      * @return a point in the middle of the bracket.
      * @see #getFMid()
      */
     double getMid() {
         return mid;
     }

     /**
      * Get function value at {@link #getMid()}.
      * @return function value at {@link #getMid()}
      */
     double getFMid() {
         return fMid;
     }

     /**
      * Get the value of the function.
      *
      * @param func Function.
      * @param x Point.
      * @return the value
      * @throws IllegalStateException if the maximal number of evaluations is exceeded.
      */
     private double value(DoubleUnaryOperator func, double x) {
         if (evaluations >= maxEvaluations) {
             throw new IllegalStateException("Too many evaluations: " + evaluations);
         }
         evaluations++;
         return func.applyAsDouble(x);
     }
 }
	/*
	* Licensed to the Apache Software Foundation (ASF) under one or more
	* contributor license agreements. See the NOTICE file distributed with
	* this work for additional information regarding copyright ownership.
	* The ASF licenses this file to You under the Apache License, Version 2.0
	* (the "License"); you may not use this file except in compliance with
	* the License. You may obtain a copy of the License at
	*
	* http://www.apache.org/licenses/LICENSE-2.0
	*
	* Unless required by applicable law or agreed to in writing, software
	* distributed under the License is distributed on an "AS IS" BASIS,
	* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
	* See the License for the specific language governing permissions and
	* limitations under the License.
	*/
	package org.apache.commons.statistics.inference;

	import java.util.function.DoubleUnaryOperator;

	/**
	* Provide an interval that brackets a local minimum of a function.
	* This code is based on a Python implementation (from <em>SciPy</em>,
	* module {@code optimize.py} v0.5).
	*
	* <p>This class has been extracted from {@code o.a.c.math4.optim.univariate}
	* and modified to: remove support for bracketing a maximum; support bounds
	* on the bracket; correct the sign of the denominator when the magnitude is small;
	* and return true/false if there is a minimum strictly inside the bounds.
	*
	* @since 1.1
	*/
	class BracketFinder {
	/** Tolerance to avoid division by zero. */
	private static final double EPS_MIN = 1e-21;
	/** Golden section. */
	private static final double GOLD = 1.6180339887498948482;
	/** Factor for expanding the interval. */
	private final double growLimit;
	/** Number of allowed function evaluations. */
	private final int maxEvaluations;
	/** Number of function evaluations performed in the last search. */
	private int evaluations;
	/** Lower bound of the bracket. */
	private double lo;
	/** Higher bound of the bracket. */
	private double hi;
	/** Point inside the bracket. */
	private double mid;
	/** Function value at {@link #lo}. */
	private double fLo;
	/** Function value at {@link #hi}. */
	private double fHi;
	/** Function value at {@link #mid}. */
	private double fMid;

	/**
	* Constructor with default values {@code 100, 100000} (see the
	* {@link #BracketFinder(double,int) other constructor}).
	*/
	BracketFinder() {
	this(100, 100000);
	}

	/**
	* Create a bracketing interval finder.
	*
	* @param growLimit Expanding factor.
	* @param maxEvaluations Maximum number of evaluations allowed for finding
	* a bracketing interval.
	* @throws IllegalArgumentException if the {@code growLimit} or {@code maxEvalutations}
	* are not strictly positive.
	*/
	BracketFinder(double growLimit, int maxEvaluations) {
	Arguments.checkStrictlyPositive(growLimit);
	Arguments.checkStrictlyPositive(maxEvaluations);
	this.growLimit = growLimit;
	this.maxEvaluations = maxEvaluations;
	}

	/**
	* Search downhill from the initial points to obtain new points that bracket a local
	* minimum of the function. Note that the initial points do not have to bracket a minimum.
	* An exception is raised if a minimum cannot be found within the configured number
	* of function evaluations.
	*
	* <p>The bracket is limited to the provided bounds if they create a positive interval
	* {@code min < max}. It is possible that the middle of the bracket is at the bounds as
	* the final bracket is {@code f(mid) <= min(f(lo), f(hi))} and {@code lo <= mid <= hi}.
	*
	* <p>No exception is raised if the initial points are not within the bounds; the points
	* are updated to be within the bounds.
	*
	* <p>No exception is raised if the initial points are equal; the bracket will be returned
	* as a single point {@code lo == mid == hi}.
	*
	* @param func Function whose optimum should be bracketed.
	* @param a Initial point.
	* @param b Initial point.
	* @param min Minimum bound of the bracket (inclusive).
	* @param max Maximum bound of the bracket (inclusive).
	* @return true if the mid-point is strictly within the final bracket {@code [lo, hi]};
	* false if there is no local minima.
	* @throws IllegalStateException if the maximum number of evaluations is exceeded.
	*/
	boolean search(DoubleUnaryOperator func,
	double a, double b,
	double min, double max) {
	evaluations = 0;

	// Limit the range of x
	DoubleUnaryOperator range;
	if (min < max) {
	// Limit: min <= x <= max
	range = x -> {
	if (x > min) {
	return x < max ? x : max;
	}
	return min;
	};
	} else {
	range = DoubleUnaryOperator.identity();
	}

	double xA = range.applyAsDouble(a);
	double xB = range.applyAsDouble(b);
	double fA = value(func, xA);
	double fB = value(func, xB);
	// Ensure fB <= fA
	if (fA < fB) {
	double tmp = xA;
	xA = xB;
	xB = tmp;
	tmp = fA;
	fA = fB;
	fB = tmp;
	}

	double xC = range.applyAsDouble(xB + GOLD * (xB - xA));
	double fC = value(func, xC);

	// Note: When a [min, max] interval is provided and there is no minima then this
	// loop will terminate when B == C and both are at the min/max bound.
	while (fC < fB) {
	final double tmp1 = (xB - xA) * (fB - fC);
	final double tmp2 = (xB - xC) * (fB - fA);

	final double val = tmp2 - tmp1;
	// limit magnitude of val to a small value
	final double denom = 2 * Math.copySign(Math.max(Math.abs(val), EPS_MIN), val);

	double w = range.applyAsDouble(xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom);
	final double wLim = range.applyAsDouble(xB + growLimit * (xC - xB));

	double fW;
	if ((w - xC) * (xB - w) > 0) {
	// xB < w < xC
	fW = value(func, w);
	if (fW < fC) {
	// minimum in [xB, xC]
	xA = xB;
	xB = w;
	fA = fB;
	fB = fW;
	break;
	} else if (fW > fB) {
	// minimum in [xA, w]
	xC = w;
	fC = fW;
	break;
	}
	// continue downhill
	w = range.applyAsDouble(xC + GOLD * (xC - xB));
	fW = value(func, w);
	} else if ((w - wLim) * (xC - w) > 0) {
	// xC < w < limit
	fW = value(func, w);
	if (fW < fC) {
	// continue downhill
	xB = xC;
	xC = w;
	w = range.applyAsDouble(xC + GOLD * (xC - xB));
	fB = fC;
	fC = fW;
	fW = value(func, w);
	}
	} else if ((w - wLim) * (wLim - xC) >= 0) {
	// xC <= limit <= w
	w = wLim;
	fW = value(func, w);
	} else {
	// possibly w == xC; reject w and take a default step
	w = range.applyAsDouble(xC + GOLD * (xC - xB));
	fW = value(func, w);
	}

	xA = xB;
	fA = fB;
	xB = xC;
	fB = fC;
	xC = w;
	fC = fW;
	}

	mid = xB;
	fMid = fB;

	// Store the bracket: lo <= mid <= hi
	if (xC < xA) {
	lo = xC;
	fLo = fC;
	hi = xA;
	fHi = fA;
	} else {
	lo = xA;
	fLo = fA;
	hi = xC;
	fHi = fC;
	}

	return lo < mid && mid < hi;
	}

	/**
	* @return the number of evaluations.
	*/
	int getEvaluations() {
	return evaluations;
	}

	/**
	* @return the lower bound of the bracket.
	* @see #getFLo()
	*/
	double getLo() {
	return lo;
	}

	/**
	* Get function value at {@link #getLo()}.
	* @return function value at {@link #getLo()}
	*/
	double getFLo() {
	return fLo;
	}

	/**
	* @return the higher bound of the bracket.
	* @see #getFHi()
	*/
	double getHi() {
	return hi;
	}

	/**
	* Get function value at {@link #getHi()}.
	* @return function value at {@link #getHi()}
	*/
	double getFHi() {
	return fHi;
	}

	/**
	* @return a point in the middle of the bracket.
	* @see #getFMid()
	*/
	double getMid() {
	return mid;
	}

	/**
	* Get function value at {@link #getMid()}.
	* @return function value at {@link #getMid()}
	*/
	double getFMid() {
	return fMid;
	}

	/**
	* Get the value of the function.
	*
	* @param func Function.
	* @param x Point.
	* @return the value
	* @throws IllegalStateException if the maximal number of evaluations is exceeded.
	*/
	private double value(DoubleUnaryOperator func, double x) {
	if (evaluations >= maxEvaluations) {
	throw new IllegalStateException("Too many evaluations: " + evaluations);
	}
	evaluations++;
	return func.applyAsDouble(x);
	}
	}