org.apache.commons.math.ConvergenceException Java Examples

public void testPowell()
  throws CostException, ConvergenceException {

  CostFunction powell =
    new CostFunction() {
      public double cost(double[] x) {
        ++count;
        double a = x[0] + 10 * x[1];
        double b = x[2] - x[3];
        double c = x[1] - 2 * x[2];
        double d = x[0] - x[3];
        return a * a + 5 * b * b + c * c * c * c + 10 * d * d * d * d;
      }
    };

  count = 0;
  PointCostPair optimum =
    new MultiDirectional().minimize(powell, 1000, new ValueChecker(1.0e-3),
                                    new double[] {  3.0, -1.0, 0.0, 1.0 },
                                    new double[] {  4.0,  0.0, 1.0, 2.0 });
  assertTrue(count > 850);
  assertTrue(optimum.getCost() > 0.015);

}
 

public void testRosenbrock()
  throws CostException, ConvergenceException {

  CostFunction rosenbrock =
    new CostFunction() {
      public double cost(double[] x) {
        ++count;
        double a = x[1] - x[0] * x[0];
        double b = 1.0 - x[0];
        return 100 * a * a + b * b;
      }
    };

  count = 0;
  PointCostPair optimum =
    new MultiDirectional().minimize(rosenbrock, 100, new ValueChecker(1.0e-3),
                                    new double[][] {
                                      { -1.2,  1.0 }, { 0.9, 1.2 } , {  3.5, -2.3 }
                                    });

  assertTrue(count > 60);
  assertTrue(optimum.getCost() > 0.01);

}
 

/**
  * This method attempts to find two values a and b satisfying <ul>
  * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
  * <li> <code> f(a) * f(b) <= 0 </code> </li>
  * </ul>
  * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
  * and <code>b</code> bracket a root of f.
  * <p>
  * The algorithm starts by setting 
  * <code>a := initial -1; b := initial +1,</code> examines the value of the
  * function at <code>a</code> and <code>b</code> and keeps moving
  * the endpoints out by one unit each time through a loop that terminates 
  * when one of the following happens: <ul>
  * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
  * <li> <code> a = lower </code> and <code> b = upper</code> 
  * -- ConvergenceException </li>
  * <li> <code> maximumIterations</code> iterations elapse 
  * -- ConvergenceException </li></ul></p>
  * 
  * @param function the function
  * @param initial initial midpoint of interval being expanded to
  * bracket a root
  * @param lowerBound lower bound (a is never lower than this value)
  * @param upperBound upper bound (b never is greater than this
  * value)
  * @param maximumIterations maximum number of iterations to perform
  * @return a two element array holding {a, b}.
  * @throws ConvergenceException if the algorithm fails to find a and b
  * satisfying the desired conditions
  * @throws FunctionEvaluationException if an error occurs evaluating the 
  * function
  * @throws IllegalArgumentException if function is null, maximumIterations
  * is not positive, or initial is not between lowerBound and upperBound
  */
 public static double[] bracket(UnivariateRealFunction function,
         double initial, double lowerBound, double upperBound, 
         int maximumIterations) throws ConvergenceException, 
         FunctionEvaluationException {
     
     if (function == null) {
         throw MathRuntimeException.createIllegalArgumentException("function is null");
     }
     if (maximumIterations <= 0)  {
         throw MathRuntimeException.createIllegalArgumentException(
               "bad value for maximum iterations number: {0}", maximumIterations);
     }
     if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
         throw MathRuntimeException.createIllegalArgumentException(
               "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
               lowerBound, initial, upperBound);
     }
     double a = initial;
     double b = initial;
     double fa;
     double fb;
     int numIterations = 0 ;
 
     do {
         a = Math.max(a - 1.0, lowerBound);
         b = Math.min(b + 1.0, upperBound);
         fa = function.value(a);
         
         fb = function.value(b);
         numIterations++ ;
     } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
             ((a > lowerBound) || (b < upperBound)));

    if (maximumIterations <= 0) {
throw org.apache.commons.math.MathRuntimeException.createIllegalArgumentException("bad value for maximum iterations number: {0}", maximumIterations);
}
     
     return new double[]{a, b};
 }
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    
    return new double[]{a, b};
}
 

public void testBracketCornerSolution() throws MathException {
    try {
        UnivariateRealSolverUtils.bracket(sin, 1.5, 0, 2.0); 
        fail("Expecting ConvergenceException");
    } catch (ConvergenceException ex) {
        // expected
    }
}
 

/** Minimizes a cost function.
 * <p>The simplex is built from all its vertices.</p>
 * <p>The optimization is performed in multi-start mode.</p>
 * @param f cost function
 * @param maxEvaluations maximal number of function calls for each
 * start (note that the number will be checked <em>after</em>
 * complete simplices have been evaluated, this means that in some
 * cases this number will be exceeded by a few units, depending on
 * the dimension of the problem)
 * @param checker object to use to check for convergence
 * @param vertices array containing all vertices of the simplex
 * @param starts number of starts to perform (including the
 * first one), multi-start is disabled if value is less than or
 * equal to 1
 * @param seed seed for the random vector generator
 * @return the point/cost pairs giving the minimal cost
 * @exception NotPositiveDefiniteMatrixException if the vertices
 * array is degenerated
 * @exception CostException if the cost function throws one during
 * the search
 * @exception ConvergenceException if none of the starts did
 * converge (it is not thrown if at least one start did converge)
 */
public PointCostPair minimize(CostFunction f, int maxEvaluations,
                              ConvergenceChecker checker,
                              double[][] vertices,
                              int starts, long seed)
throws NotPositiveDefiniteMatrixException,
CostException, ConvergenceException {

    try {
        // store the points into the simplex
        buildSimplex(vertices);

        // compute the statistical properties of the simplex points
        VectorialMean meanStat = new VectorialMean(vertices[0].length);
        VectorialCovariance covStat = new VectorialCovariance(vertices[0].length, true);
        for (int i = 0; i < vertices.length; ++i) {
            meanStat.increment(vertices[i]);
            covStat.increment(vertices[i]);
        }
        double[] mean = meanStat.getResult();
        RealMatrix covariance = covStat.getResult();
        

        RandomGenerator rg = new JDKRandomGenerator();
        rg.setSeed(seed);
        RandomVectorGenerator rvg =
            new CorrelatedRandomVectorGenerator(mean,
                                                covariance, 1.0e-12 * covariance.getNorm(),
                                                new UniformRandomGenerator(rg));
        setMultiStart(starts, rvg);

        // compute minimum
        return minimize(f, maxEvaluations, checker);

    } catch (DimensionMismatchException dme) {
        // this should not happen
        throw new RuntimeException("internal error");
    }

}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    
    return new double[]{a, b};
}
 

/**
  * This method attempts to find two values a and b satisfying <ul>
  * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
  * <li> <code> f(a) * f(b) <= 0 </code> </li>
  * </ul>
  * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
  * and <code>b</code> bracket a root of f.
  * <p>
  * The algorithm starts by setting 
  * <code>a := initial -1; b := initial +1,</code> examines the value of the
  * function at <code>a</code> and <code>b</code> and keeps moving
  * the endpoints out by one unit each time through a loop that terminates 
  * when one of the following happens: <ul>
  * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
  * <li> <code> a = lower </code> and <code> b = upper</code> 
  * -- ConvergenceException </li>
  * <li> <code> maximumIterations</code> iterations elapse 
  * -- ConvergenceException </li></ul></p>
  * 
  * @param function the function
  * @param initial initial midpoint of interval being expanded to
  * bracket a root
  * @param lowerBound lower bound (a is never lower than this value)
  * @param upperBound upper bound (b never is greater than this
  * value)
  * @param maximumIterations maximum number of iterations to perform
  * @return a two element array holding {a, b}.
  * @throws ConvergenceException if the algorithm fails to find a and b
  * satisfying the desired conditions
  * @throws FunctionEvaluationException if an error occurs evaluating the 
  * function
  * @throws IllegalArgumentException if function is null, maximumIterations
  * is not positive, or initial is not between lowerBound and upperBound
  */
 public static double[] bracket(UnivariateRealFunction function,
         double initial, double lowerBound, double upperBound, 
         int maximumIterations) throws ConvergenceException, 
         FunctionEvaluationException {
     
     if (function == null) {
         throw MathRuntimeException.createIllegalArgumentException("function is null");
     }
     if (maximumIterations <= 0)  {
         throw MathRuntimeException.createIllegalArgumentException(
               "bad value for maximum iterations number: {0}", maximumIterations);
     }
     if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
         throw MathRuntimeException.createIllegalArgumentException(
               "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
               lowerBound, initial, upperBound);
     }
     double a = initial;
     double b = initial;
     double fa;
     double fb;
     int numIterations = 0 ;
 
     do {
         a = Math.max(a - 1.0, lowerBound);
         b = Math.min(b + 1.0, upperBound);
         fa = function.value(a);
         
         fb = function.value(b);
         numIterations++ ;
     } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
             ((a > lowerBound) || (b < upperBound)));

    if (function == null) {
throw org.apache.commons.math.MathRuntimeException.createIllegalArgumentException("function is null");
}
     
     return new double[]{a, b};
 }
 

private void checkIntegerOverflow(double a) {
    try {
        new Fraction(a, 1.0e-12, 1000);
        fail("an exception should have been thrown");
    } catch (ConvergenceException ce) {
        // expected behavior
    } catch (Exception e) {
        fail("wrong exception caught");
    }
}
 

/** Minimizes a cost function.
 * <p>The initial simplex is built from two vertices that are
 * considered to represent two opposite vertices of a box parallel
 * to the canonical axes of the space. The simplex is the subset of
 * vertices encountered while going from vertexA to vertexB
 * traveling along the box edges only. This can be seen as a scaled
 * regular simplex using the projected separation between the given
 * points as the scaling factor along each coordinate axis.</p>
 * <p>The optimization is performed in multi-start mode.</p>
 * @param f cost function
 * @param maxEvaluations maximal number of function calls for each
 * start (note that the number will be checked <em>after</em>
 * complete simplices have been evaluated, this means that in some
 * cases this number will be exceeded by a few units, depending on
 * the dimension of the problem)
 * @param checker object to use to check for convergence
 * @param vertexA first vertex
 * @param vertexB last vertex
 * @param starts number of starts to perform (including the
 * first one), multi-start is disabled if value is less than or
 * equal to 1
 * @param seed seed for the random vector generator
 * @return the point/cost pairs giving the minimal cost
 * @exception CostException if the cost function throws one during
 * the search
 * @exception ConvergenceException if none of the starts did
 * converge (it is not thrown if at least one start did converge)
 */
public PointCostPair minimize(CostFunction f, int maxEvaluations,
                              ConvergenceChecker checker,
                              double[] vertexA, double[] vertexB,
                              int starts, long seed)
throws CostException, ConvergenceException {

    // set up the simplex traveling around the box
    buildSimplex(vertexA, vertexB);

    // we consider the simplex could have been produced by a generator
    // having its mean value at the center of the box, the standard
    // deviation along each axe being the corresponding half size
    double[] mean              = new double[vertexA.length];
    double[] standardDeviation = new double[vertexA.length];
    for (int i = 0; i < vertexA.length; ++i) {
        mean[i]              = 0.5 * (vertexA[i] + vertexB[i]);
        standardDeviation[i] = 0.5 * Math.abs(vertexA[i] - vertexB[i]);
    }

    RandomGenerator rg = new JDKRandomGenerator();
    rg.setSeed(seed);
    UniformRandomGenerator urg = new UniformRandomGenerator(rg);
    RandomVectorGenerator rvg =
        new UncorrelatedRandomVectorGenerator(mean, standardDeviation, urg);
    setMultiStart(starts, rvg);

    // compute minimum
    return minimize(f, maxEvaluations, checker);

}
 

/**
 * Find a real root in the given interval.
 * <p>
 * Despite the bracketing condition, the root returned by solve(Complex[],
 * Complex) may not be a real zero inside [min, max]. For example,
 * p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
 * another initial value, or, as we did here, call solveAll() to obtain
 * all roots and pick up the one that we're looking for.</p>
 *
 * @param min the lower bound for the interval
 * @param max the upper bound for the interval
 * @return the point at which the function value is zero
 * @throws ConvergenceException if the maximum iteration count is exceeded
 * or the solver detects convergence problems otherwise
 * @throws FunctionEvaluationException if an error occurs evaluating the
 * function 
 * @throws IllegalArgumentException if any parameters are invalid
 */
public double solve(double min, double max) throws ConvergenceException, 
    FunctionEvaluationException {

    // check for zeros before verifying bracketing
    if (p.value(min) == 0.0) { return min; }
    if (p.value(max) == 0.0) { return max; }
    verifyBracketing(min, max, p);

    double coefficients[] = p.getCoefficients();
    Complex c[] = new Complex[coefficients.length];
    for (int i = 0; i < coefficients.length; i++) {
        c[i] = new Complex(coefficients[i], 0.0);
    }
    Complex initial = new Complex(0.5 * (min + max), 0.0);
    Complex z = solve(c, initial);
    if (isRootOK(min, max, z)) {
        setResult(z.getReal(), iterationCount);
        return result;
    }

    // solve all roots and select the one we're seeking
    Complex[] root = solveAll(c, initial);
    for (int i = 0; i < root.length; i++) {
        if (isRootOK(min, max, root[i])) {
            setResult(root[i].getReal(), iterationCount);
            return result;
        }
    }

    // should never happen
    throw new ConvergenceException();
}
 

/**
     * This method attempts to find two values a and b satisfying <ul>
     * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
     * <li> <code> f(a) * f(b) <= 0 </code> </li>
     * </ul>
     * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
     * and <code>b</code> bracket a root of f.
     * <p>
     * The algorithm starts by setting 
     * <code>a := initial -1; b := initial +1,</code> examines the value of the
     * function at <code>a</code> and <code>b</code> and keeps moving
     * the endpoints out by one unit each time through a loop that terminates 
     * when one of the following happens: <ul>
     * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
     * <li> <code> a = lower </code> and <code> b = upper</code> 
     * -- ConvergenceException </li>
     * <li> <code> maximumIterations</code> iterations elapse 
     * -- ConvergenceException </li></ul></p>
     * 
     * @param function the function
     * @param initial initial midpoint of interval being expanded to
     * bracket a root
     * @param lowerBound lower bound (a is never lower than this value)
     * @param upperBound upper bound (b never is greater than this
     * value)
     * @param maximumIterations maximum number of iterations to perform
     * @return a two element array holding {a, b}.
     * @throws ConvergenceException if the algorithm fails to find a and b
     * satisfying the desired conditions
     * @throws FunctionEvaluationException if an error occurs evaluating the 
     * function
     * @throws IllegalArgumentException if function is null, maximumIterations
     * is not positive, or initial is not between lowerBound and upperBound
     */
    public static double[] bracket(UnivariateRealFunction function,
            double initial, double lowerBound, double upperBound, 
            int maximumIterations) throws ConvergenceException, 
            FunctionEvaluationException {
        
        if (function == null) {
            throw MathRuntimeException.createIllegalArgumentException("function is null");
        }
        if (maximumIterations <= 0)  {
            throw MathRuntimeException.createIllegalArgumentException(
                  "bad value for maximum iterations number: {0}", maximumIterations);
        }
        if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
            throw MathRuntimeException.createIllegalArgumentException(
                  "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
                  lowerBound, initial, upperBound);
        }
        double a = initial;
        double b = initial;
        double fa;
        double fb;
        int numIterations = 0 ;
    
        do {
            a = Math.max(a - 1.0, lowerBound);
            b = Math.min(b + 1.0, upperBound);
            fa = function.value(a);
            
            fb = function.value(b);
            numIterations++ ;
        } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
                ((a > lowerBound) || (b < upperBound)));
   
// start of generated patch
if(fa*upperBound>=0.0){
throw new ConvergenceException("number of iterations={0}, maximum iterations={1}, "+"initial={2}, lower bound={3}, upper bound={4}, final a value={5}, ",numIterations,maximumIterations,initial,lowerBound,upperBound,a,b,fa,fb);
}
// end of generated patch
/* start of original code
        if (fa * fb >= 0.0 ) {
            throw new ConvergenceException(
                      "number of iterations={0}, maximum iterations={1}, " +
                      "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                      "final b value={6}, f(a)={7}, f(b)={8}",
                      numIterations, maximumIterations, initial,
                      lowerBound, upperBound, a, b, fa, fb);
        }
 end of original code*/
        
        return new double[]{a, b};
    }
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if (fa * fb >= 0.0 ) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

/** Minimizes a cost function.
 * @param f cost function
 * @param maxEvaluations maximal number of function calls for each
 * start (note that the number will be checked <em>after</em>
 * complete simplices have been evaluated, this means that in some
 * cases this number will be exceeded by a few units, depending on
 * the dimension of the problem)
 * @param checker object to use to check for convergence
 * @return the point/cost pairs giving the minimal cost
 * @exception CostException if the cost function throws one during
 * the search
 * @exception ConvergenceException if none of the starts did
 * converge (it is not thrown if at least one start did converge)
 */
private PointCostPair minimize(CostFunction f, int maxEvaluations,
                                ConvergenceChecker checker)
throws CostException, ConvergenceException {

    this.f = f;
    minima = new PointCostPair[starts];

    // multi-start loop
    for (int i = 0; i < starts; ++i) {

        evaluations = 0;
        evaluateSimplex();

        for (boolean loop = true; loop;) {
            if (checker.converged(simplex)) {
                // we have found a minimum
                minima[i] = simplex[0];
                loop = false;
            } else if (evaluations >= maxEvaluations) {
                // this start did not converge, try a new one
                minima[i] = null;
                loop = false;
            } else {
                iterateSimplex();
            }
        }

        if (i < (starts - 1)) {
            // restart
            buildSimplex(generator);
        }

    }

    // sort the minima from lowest cost to highest cost, followed by
    // null elements
    Arrays.sort(minima, pointCostPairComparator);

    // return the found point given the lowest cost
    if (minima[0] == null) {
        throw new ConvergenceException("none of the {0} start points" +
                                       " lead to convergence",
                                       new Object[] {
                                         Integer.toString(starts)
                                       });
    }
    return minima[0];

}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if (fa * fb >= 0.0 ) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

/**
     * This method attempts to find two values a and b satisfying <ul>
     * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
     * <li> <code> f(a) * f(b) <= 0 </code> </li>
     * </ul>
     * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
     * and <code>b</code> bracket a root of f.
     * <p>
     * The algorithm starts by setting 
     * <code>a := initial -1; b := initial +1,</code> examines the value of the
     * function at <code>a</code> and <code>b</code> and keeps moving
     * the endpoints out by one unit each time through a loop that terminates 
     * when one of the following happens: <ul>
     * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
     * <li> <code> a = lower </code> and <code> b = upper</code> 
     * -- ConvergenceException </li>
     * <li> <code> maximumIterations</code> iterations elapse 
     * -- ConvergenceException </li></ul></p>
     * 
     * @param function the function
     * @param initial initial midpoint of interval being expanded to
     * bracket a root
     * @param lowerBound lower bound (a is never lower than this value)
     * @param upperBound upper bound (b never is greater than this
     * value)
     * @param maximumIterations maximum number of iterations to perform
     * @return a two element array holding {a, b}.
     * @throws ConvergenceException if the algorithm fails to find a and b
     * satisfying the desired conditions
     * @throws FunctionEvaluationException if an error occurs evaluating the 
     * function
     * @throws IllegalArgumentException if function is null, maximumIterations
     * is not positive, or initial is not between lowerBound and upperBound
     */
    public static double[] bracket(UnivariateRealFunction function,
            double initial, double lowerBound, double upperBound, 
            int maximumIterations) throws ConvergenceException, 
            FunctionEvaluationException {
        
        if (function == null) {
            throw MathRuntimeException.createIllegalArgumentException("function is null");
        }
        if (maximumIterations <= 0)  {
            throw MathRuntimeException.createIllegalArgumentException(
                  "bad value for maximum iterations number: {0}", maximumIterations);
        }
        if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
            throw MathRuntimeException.createIllegalArgumentException(
                  "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
                  lowerBound, initial, upperBound);
        }
        double a = initial;
        double b = initial;
        double fa;
        double fb;
        int numIterations = 0 ;
    
        do {
            a = Math.max(a - 1.0, lowerBound);
            b = Math.min(b + 1.0, upperBound);
            fa = function.value(a);
            
            fb = function.value(b);
            numIterations++ ;
        } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
                ((a > lowerBound) || (b < upperBound)));
   
// start of generated patch
if(fa*upperBound>=0.0){
throw new ConvergenceException("number of iterations={0}, maximum iterations={1}, "+"initial={2}, lower bound={3}, upper bound={4}, final a value={5}, ",numIterations,maximumIterations,initial,lowerBound,upperBound,a,b,fa,fb);
}
// end of generated patch
/* start of original code
        if (fa * fb >= 0.0 ) {
            throw new ConvergenceException(
                      "number of iterations={0}, maximum iterations={1}, " +
                      "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                      "final b value={6}, f(a)={7}, f(b)={8}",
                      numIterations, maximumIterations, initial,
                      lowerBound, upperBound, a, b, fa, fb);
        }
 end of original code*/
        
        return new double[]{a, b};
    }
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if (fa * fb >= 0.0 ) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if (fa * fb >= 0.0 ) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if (fa * fb >= 0.0 ) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if (fa * fb >= 0.0 ) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if (fa * fb > 0.0 ) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if (fa * fb >= 0.0 ) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

/** {@inheritDoc} */
public double optimize(final UnivariateRealFunction f, final GoalType goalType,
                       final double min, final double max, final double startValue)
        throws ConvergenceException, FunctionEvaluationException {
    return optimize(f, goalType, min, max);
}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if ((((fa * fb) > 0.0) && (numIterations < maximumIterations)) && ((a > lowerBound) || (b < upperBound))) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

public void testDoubleConstructor() throws ConvergenceException  {
    assertFraction(1, 2, new Fraction((double)1 / (double)2));
    assertFraction(1, 3, new Fraction((double)1 / (double)3));
    assertFraction(2, 3, new Fraction((double)2 / (double)3));
    assertFraction(1, 4, new Fraction((double)1 / (double)4));
    assertFraction(3, 4, new Fraction((double)3 / (double)4));
    assertFraction(1, 5, new Fraction((double)1 / (double)5));
    assertFraction(2, 5, new Fraction((double)2 / (double)5));
    assertFraction(3, 5, new Fraction((double)3 / (double)5));
    assertFraction(4, 5, new Fraction((double)4 / (double)5));
    assertFraction(1, 6, new Fraction((double)1 / (double)6));
    assertFraction(5, 6, new Fraction((double)5 / (double)6));
    assertFraction(1, 7, new Fraction((double)1 / (double)7));
    assertFraction(2, 7, new Fraction((double)2 / (double)7));
    assertFraction(3, 7, new Fraction((double)3 / (double)7));
    assertFraction(4, 7, new Fraction((double)4 / (double)7));
    assertFraction(5, 7, new Fraction((double)5 / (double)7));
    assertFraction(6, 7, new Fraction((double)6 / (double)7));
    assertFraction(1, 8, new Fraction((double)1 / (double)8));
    assertFraction(3, 8, new Fraction((double)3 / (double)8));
    assertFraction(5, 8, new Fraction((double)5 / (double)8));
    assertFraction(7, 8, new Fraction((double)7 / (double)8));
    assertFraction(1, 9, new Fraction((double)1 / (double)9));
    assertFraction(2, 9, new Fraction((double)2 / (double)9));
    assertFraction(4, 9, new Fraction((double)4 / (double)9));
    assertFraction(5, 9, new Fraction((double)5 / (double)9));
    assertFraction(7, 9, new Fraction((double)7 / (double)9));
    assertFraction(8, 9, new Fraction((double)8 / (double)9));
    assertFraction(1, 10, new Fraction((double)1 / (double)10));
    assertFraction(3, 10, new Fraction((double)3 / (double)10));
    assertFraction(7, 10, new Fraction((double)7 / (double)10));
    assertFraction(9, 10, new Fraction((double)9 / (double)10));
    assertFraction(1, 11, new Fraction((double)1 / (double)11));
    assertFraction(2, 11, new Fraction((double)2 / (double)11));
    assertFraction(3, 11, new Fraction((double)3 / (double)11));
    assertFraction(4, 11, new Fraction((double)4 / (double)11));
    assertFraction(5, 11, new Fraction((double)5 / (double)11));
    assertFraction(6, 11, new Fraction((double)6 / (double)11));
    assertFraction(7, 11, new Fraction((double)7 / (double)11));
    assertFraction(8, 11, new Fraction((double)8 / (double)11));
    assertFraction(9, 11, new Fraction((double)9 / (double)11));
    assertFraction(10, 11, new Fraction((double)10 / (double)11));
}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if (fa * fb > 0.0 ) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

public void testPdf() throws FunctionEvaluationException  {

        final int numberOfTests = 100;
        double totErr = 0;
        double ptotErr = 0; int np = 0;
        double qtotErr = 0;

        Random random = new Random(37);

        for(int i = 0; i < numberOfTests; i++){
            final double mean = .01 + (3-0.01) * random.nextDouble();
            final double var = .01 + (3-0.01) * random.nextDouble();

            final double scale = var / mean;
            final double shape = mean / scale;

            final GammaDistribution gamma = new GammaDistribution(shape,scale);

            final double value = gamma.nextGamma();

            final double mypdf = mypdf(value, shape, scale);
            final double pdf = gamma.pdf(value);
            if ( Double.isInfinite(mypdf) && Double.isInfinite(pdf)) {
                continue;
            }

            assertFalse(Double.isNaN(mypdf));
            assertFalse(Double.isNaN(pdf));

            totErr +=  mypdf != 0 ? Math.abs((pdf - mypdf)/mypdf) : pdf;

            assertFalse("nan", Double.isNaN(totErr));
            //assertEquals("" + shape + "," + scale + "," + value, mypdf,gamma.pdf(value),1e-10);

            final double cdf = gamma.cdf(value);
            UnivariateRealFunction f = new UnivariateRealFunction() {
                public double value(double v) throws FunctionEvaluationException {
                    return mypdf(v, shape, scale);
                }
            };
            final UnivariateRealIntegrator integrator = new RombergIntegrator();
            integrator.setAbsoluteAccuracy(1e-14);
            integrator.setMaximalIterationCount(16);  // fail if it takes too much time

            double x;
            try {
                x = integrator.integrate(f, 0, value);
                ptotErr += cdf != 0.0 ? Math.abs(x-cdf)/cdf : x;
                np += 1;
                //assertTrue("" + shape + "," + scale + "," + value + " " + Math.abs(x-cdf)/x + "> 1e-6", Math.abs(1-cdf/x) < 1e-6);

                //System.out.println(shape + ","  + scale + " " + value);
            } catch( ConvergenceException e ) {
                 // can't integrate , skip test
              //  System.out.println(shape + ","  + scale + " skipped");
            }

            final double q = gamma.quantile(cdf);
            qtotErr += q != 0 ? Math.abs(q-value)/q : value;
           // assertEquals("" + shape + "," + scale + "," + value + " " + Math.abs(q-value)/value, q, value, 1e-6);
        }
        //System.out.println( !Double.isNaN(totErr) );
       // System.out.println(totErr);
        // bad test, but I can't find a good threshold that works for all individual cases 
        assertTrue("failed " + totErr/numberOfTests, totErr/numberOfTests < 1e-7);
        assertTrue("failed " + ptotErr/np, np > 0 ? (ptotErr/np < 1e-5) : true);
        assertTrue("failed " + qtotErr/numberOfTests , qtotErr/numberOfTests < 1e-7);
	}
 

/**
 * This method attempts to find two values a and b satisfying <ul>
 * <li> <code> lowerBound <= a < initial < b <= upperBound</code> </li>
 * <li> <code> f(a) * f(b) <= 0 </code> </li>
 * </ul>
 * If f is continuous on <code>[a,b],</code> this means that <code>a</code>
 * and <code>b</code> bracket a root of f.
 * <p>
 * The algorithm starts by setting 
 * <code>a := initial -1; b := initial +1,</code> examines the value of the
 * function at <code>a</code> and <code>b</code> and keeps moving
 * the endpoints out by one unit each time through a loop that terminates 
 * when one of the following happens: <ul>
 * <li> <code> f(a) * f(b) <= 0 </code> --  success!</li>
 * <li> <code> a = lower </code> and <code> b = upper</code> 
 * -- ConvergenceException </li>
 * <li> <code> maximumIterations</code> iterations elapse 
 * -- ConvergenceException </li></ul></p>
 * 
 * @param function the function
 * @param initial initial midpoint of interval being expanded to
 * bracket a root
 * @param lowerBound lower bound (a is never lower than this value)
 * @param upperBound upper bound (b never is greater than this
 * value)
 * @param maximumIterations maximum number of iterations to perform
 * @return a two element array holding {a, b}.
 * @throws ConvergenceException if the algorithm fails to find a and b
 * satisfying the desired conditions
 * @throws FunctionEvaluationException if an error occurs evaluating the 
 * function
 * @throws IllegalArgumentException if function is null, maximumIterations
 * is not positive, or initial is not between lowerBound and upperBound
 */
public static double[] bracket(UnivariateRealFunction function,
        double initial, double lowerBound, double upperBound, 
        int maximumIterations) throws ConvergenceException, 
        FunctionEvaluationException {
    
    if (function == null) {
        throw MathRuntimeException.createIllegalArgumentException("function is null");
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              "bad value for maximum iterations number: {0}", maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              "invalid bracketing parameters:  lower bound={0},  initial={1}, upper bound={2}",
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = Math.max(a - 1.0, lowerBound);
        b = Math.min(b + 1.0, upperBound);
        fa = function.value(a);
        
        fb = function.value(b);
        numIterations++ ;
    } while ((fa * fb > 0.0) && (numIterations < maximumIterations) && 
            ((a > lowerBound) || (b < upperBound)));
   
    if (fa * fb >= 0.0 ) {
        throw new ConvergenceException(
                  "number of iterations={0}, maximum iterations={1}, " +
                  "initial={2}, lower bound={3}, upper bound={4}, final a value={5}, " +
                  "final b value={6}, f(a)={7}, f(b)={8}",
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }
    
    return new double[]{a, b};
}
 

/**
 * Convenience method to find a zero of a univariate real function.  A default
 * solver is used. 
 * 
 * @param f the function
 * @param x0 the lower bound for the interval
 * @param x1 the upper bound for the interval
 * @param absoluteAccuracy the accuracy to be used by the solver
 * @return a value where the function is zero
 * @throws ConvergenceException if the iteration count is exceeded
 * @throws FunctionEvaluationException if an error occurs evaluating the
 * function
 * @throws IllegalArgumentException if f is null, the endpoints do not 
 * specify a valid interval, or the absoluteAccuracy is not valid for the
 * default solver
 */
public static double solve(UnivariateRealFunction f, double x0, double x1,
        double absoluteAccuracy) throws ConvergenceException, 
        FunctionEvaluationException {    
   
    setup(f);
    UnivariateRealSolver solver = LazyHolder.FACTORY.newDefaultSolver();
    solver.setAbsoluteAccuracy(absoluteAccuracy);
    return solver.solve(f, x0, x1);
}
 

/**
 * Convenience method to find a zero of a univariate real function.  A default
 * solver is used. 
 * 
 * @param f the function
 * @param x0 the lower bound for the interval
 * @param x1 the upper bound for the interval
 * @param absoluteAccuracy the accuracy to be used by the solver
 * @return a value where the function is zero
 * @throws ConvergenceException if the iteration count is exceeded
 * @throws FunctionEvaluationException if an error occurs evaluating the
 * function
 * @throws IllegalArgumentException if f is null, the endpoints do not 
 * specify a valid interval, or the absoluteAccuracy is not valid for the
 * default solver
 */
public static double solve(UnivariateRealFunction f, double x0, double x1,
        double absoluteAccuracy) throws ConvergenceException, 
        FunctionEvaluationException {    
   
    setup(f);
    UnivariateRealSolver solver = LazyHolder.FACTORY.newDefaultSolver();
    solver.setAbsoluteAccuracy(absoluteAccuracy);
    return solver.solve(f, x0, x1);
}