Java Code Examples for org.apache.commons.math.util.FastMath#max()

/**
 * Test of solver for the exponential function.
 * <p>
 * It takes 10 to 15 iterations for the last two tests to converge.
 * In fact, if not for the bisection alternative, the solver would
 * exceed the default maximal iteration of 100.
 */
@Test
public void testExpm1Function() {
    UnivariateRealFunction f = new Expm1Function();
    UnivariateRealSolver solver = new MullerSolver();
    double min, max, expected, result, tolerance;

    min = -1.0; max = 2.0; expected = 0.0;
    tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
                FastMath.abs(expected * solver.getRelativeAccuracy()));
    result = solver.solve(100, f, min, max);
    Assert.assertEquals(expected, result, tolerance);

    min = -20.0; max = 10.0; expected = 0.0;
    tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
                FastMath.abs(expected * solver.getRelativeAccuracy()));
    result = solver.solve(100, f, min, max);
    Assert.assertEquals(expected, result, tolerance);

    min = -50.0; max = 100.0; expected = 0.0;
    tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
                FastMath.abs(expected * solver.getRelativeAccuracy()));
    result = solver.solve(100, f, min, max);
    Assert.assertEquals(expected, result, tolerance);
}
 

/** Perform some sanity checks on the integration parameters.
 * @param ode differential equations set
 * @param t0 start time
 * @param y0 state vector at t0
 * @param t target time for the integration
 * @param y placeholder where to put the state vector
 * @exception IntegratorException if some inconsistency is detected
 */
protected void sanityChecks(final FirstOrderDifferentialEquations ode,
                            final double t0, final double[] y0,
                            final double t, final double[] y)
    throws IntegratorException {

    if (ode.getDimension() != y0.length) {
        throw new IntegratorException(
                LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, ode.getDimension(), y0.length);
    }

    if (ode.getDimension() != y.length) {
        throw new IntegratorException(
                LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, ode.getDimension(), y.length);
    }

    if (FastMath.abs(t - t0) <= 1.0e-12 * FastMath.max(FastMath.abs(t0), FastMath.abs(t))) {
        throw new IntegratorException(
                LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
                FastMath.abs(t - t0));
    }

}
 

/** {@inheritDoc} */
@Override
public double getNorm() {
    final double[] colSums = new double[BLOCK_SIZE];
    double maxColSum = 0;
    for (int jBlock = 0; jBlock < blockColumns; jBlock++) {
        final int jWidth = blockWidth(jBlock);
        Arrays.fill(colSums, 0, jWidth, 0.0);
        for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
            final int iHeight = blockHeight(iBlock);
            final double[] block = blocks[iBlock * blockColumns + jBlock];
            for (int j = 0; j < jWidth; ++j) {
                double sum = 0;
                for (int i = 0; i < iHeight; ++i) {
                    sum += FastMath.abs(block[i * jWidth + j]);
                }
                colSums[j] += sum;
            }
        }
        for (int j = 0; j < jWidth; ++j) {
            maxColSum = FastMath.max(maxColSum, colSums[j]);
        }
    }
    return maxColSum;
}
 

/**
 * Check if a matrix is symmetric.
 *
 * @param matrix Matrix to check.
 * @param raiseException If {@code true}, the method will throw an
 * exception if {@code matrix} is not symmetric.
 * @return {@code true} if {@code matrix} is symmetric.
 * @throws NonSymmetricMatrixException if the matrix is not symmetric and
 * {@code raiseException} is {@code true}.
 */
private boolean isSymmetric(final RealMatrix matrix,
                            boolean raiseException) {
    final int rows = matrix.getRowDimension();
    final int columns = matrix.getColumnDimension();
    final double eps = 10 * rows * columns * MathUtils.EPSILON;
    for (int i = 0; i < rows; ++i) {
        for (int j = i + 1; j < columns; ++j) {
            final double mij = matrix.getEntry(i, j);
            final double mji = matrix.getEntry(j, i);
            if (FastMath.abs(mij - mji) >
                (FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) {
                if (raiseException) {
                    throw new NonSymmetricMatrixException(i, j, eps);
                }
                return false;
            }
        }
    }
    return true;
}
 

/**
 * Test of solver for the sine function.
 */
public void testSinFunction() throws MathException {
    UnivariateRealFunction f = new SinFunction();
    UnivariateRealSolver solver = new MullerSolver();
    double min, max, expected, result, tolerance;

    min = 3.0; max = 4.0; expected = FastMath.PI;
    tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
                FastMath.abs(expected * solver.getRelativeAccuracy()));
    result = solver.solve(f, min, max);
    assertEquals(expected, result, tolerance);

    min = -1.0; max = 1.5; expected = 0.0;
    tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
                FastMath.abs(expected * solver.getRelativeAccuracy()));
    result = solver.solve(f, min, max);
    assertEquals(expected, result, tolerance);
}
 

/** {@inheritDoc} */
@Override
public double getNorm() {
    final double[] colSums = new double[BLOCK_SIZE];
    double maxColSum = 0;
    for (int jBlock = 0; jBlock < blockColumns; jBlock++) {
        final int jWidth = blockWidth(jBlock);
        Arrays.fill(colSums, 0, jWidth, 0.0);
        for (int iBlock = 0; iBlock < blockRows; ++iBlock) {
            final int iHeight = blockHeight(iBlock);
            final double[] block = blocks[iBlock * blockColumns + jBlock];
            for (int j = 0; j < jWidth; ++j) {
                double sum = 0;
                for (int i = 0; i < iHeight; ++i) {
                    sum += FastMath.abs(block[i * jWidth + j]);
                }
                colSums[j] += sum;
            }
        }
        for (int j = 0; j < jWidth; ++j) {
            maxColSum = FastMath.max(maxColSum, colSums[j]);
        }
    }
    return maxColSum;
}
 

/**
 * Verifies that probability computations are consistent
 */
@Test
public void testConsistency() throws Exception {
    for (int i=1; i < cumulativeTestPoints.length; i++) {

        // check that cdf(x, x) = 0
        TestUtils.assertEquals(0d,
           distribution.cumulativeProbability
             (cumulativeTestPoints[i], cumulativeTestPoints[i]), tolerance);

        // check that P(a < X < b) = P(X < b) - P(X < a)
        double upper = FastMath.max(cumulativeTestPoints[i], cumulativeTestPoints[i -1]);
        double lower = FastMath.min(cumulativeTestPoints[i], cumulativeTestPoints[i -1]);
        double diff = distribution.cumulativeProbability(upper) -
            distribution.cumulativeProbability(lower);
        double direct = distribution.cumulativeProbability(lower, upper);
        TestUtils.assertEquals("Inconsistent cumulative probabilities for ("
                + lower + "," + upper + ")", diff, direct, tolerance);
    }
}
 

/** {@inheritDoc} */
@Override
public double walkInOptimizedOrder(final RealMatrixPreservingVisitor visitor,
                                   final int startRow, final int endRow,
                                   final int startColumn, final int endColumn) {
    MatrixUtils.checkSubMatrixIndex(this, startRow, endRow, startColumn, endColumn);
    visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
    for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
        final int p0 = iBlock * BLOCK_SIZE;
        final int pStart = FastMath.max(startRow, p0);
        final int pEnd = FastMath.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
        for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
            final int jWidth = blockWidth(jBlock);
            final int q0 = jBlock * BLOCK_SIZE;
            final int qStart = FastMath.max(startColumn, q0);
            final int qEnd = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
            final double[] block = blocks[iBlock * blockColumns + jBlock];
            for (int p = pStart; p < pEnd; ++p) {
                int k = (p - p0) * jWidth + qStart - q0;
                for (int q = qStart; q < qEnd; ++q) {
                    visitor.visit(p, q, block[k]);
                    ++k;
                }
            }
        }
    }
    return visitor.end();
}
 

/** {@inheritDoc} */
@Override
public double walkInOptimizedOrder(final RealMatrixChangingVisitor visitor,
                                   final int startRow, final int endRow,
                                   final int startColumn, final int endColumn) {
    MatrixUtils.checkSubMatrixIndex(this, startRow, endRow, startColumn, endColumn);
    visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
    for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
        final int p0 = iBlock * BLOCK_SIZE;
        final int pStart = FastMath.max(startRow, p0);
        final int pEnd = FastMath.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
        for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
            final int jWidth = blockWidth(jBlock);
            final int q0 = jBlock * BLOCK_SIZE;
            final int qStart = FastMath.max(startColumn, q0);
            final int qEnd = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
            final double[] block = blocks[iBlock * blockColumns + jBlock];
            for (int p = pStart; p < pEnd; ++p) {
                int k = (p - p0) * jWidth + qStart - q0;
                for (int q = qStart; q < qEnd; ++q) {
                    block[k] = visitor.visit(p, q, block[k]);
                    ++k;
                }
            }
        }
    }
    return visitor.end();
}
 

/** {@inheritDoc} */
@Override
public T walkInRowOrder(final FieldMatrixPreservingVisitor<T> visitor,
                        final int startRow, final int endRow,
                        final int startColumn, final int endColumn) {
    checkSubMatrixIndex(startRow, endRow, startColumn, endColumn);
    visitor.start(rows, columns, startRow, endRow, startColumn, endColumn);
    for (int iBlock = startRow / BLOCK_SIZE; iBlock < 1 + endRow / BLOCK_SIZE; ++iBlock) {
        final int p0     = iBlock * BLOCK_SIZE;
        final int pStart = FastMath.max(startRow, p0);
        final int pEnd   = FastMath.min((iBlock + 1) * BLOCK_SIZE, 1 + endRow);
        for (int p = pStart; p < pEnd; ++p) {
            for (int jBlock = startColumn / BLOCK_SIZE; jBlock < 1 + endColumn / BLOCK_SIZE; ++jBlock) {
                final int jWidth = blockWidth(jBlock);
                final int q0     = jBlock * BLOCK_SIZE;
                final int qStart = FastMath.max(startColumn, q0);
                final int qEnd   = FastMath.min((jBlock + 1) * BLOCK_SIZE, 1 + endColumn);
                final T[] block = blocks[iBlock * blockColumns + jBlock];
                int k = (p - p0) * jWidth + qStart - q0;
                for (int q = qStart; q < qEnd; ++q) {
                    visitor.visit(p, q, block[k]);
                    ++k;
                }
            }
         }
    }
    return visitor.end();
}
 

/** Reinitialize the beginning of the step.
 * @param interpolator valid for the current step
 * @exception EventException if the event handler
 * value cannot be evaluated at the beginning of the step
 */
public void reinitializeBegin(final StepInterpolator interpolator)
    throws EventException {
    try {

        t0 = interpolator.getPreviousTime();
        interpolator.setInterpolatedTime(t0);
        g0 = handler.g(t0, interpolator.getInterpolatedState());
        if (g0 == 0) {
            // excerpt from MATH-421 issue:
            // If an ODE solver is setup with an EventHandler that return STOP
            // when the even is triggered, the integrator stops (which is exactly
            // the expected behavior). If however the user wants to restart the
            // solver from the final state reached at the event with the same
            // configuration (expecting the event to be triggered again at a
            // later time), then the integrator may fail to start. It can get stuck
            // at the previous event. The use case for the bug MATH-421 is fairly
            // general, so events occurring exactly at start in the first step should
            // be ignored.

            // extremely rare case: there is a zero EXACTLY at interval start
            // we will use the sign slightly after step beginning to force ignoring this zero
            final double epsilon = FastMath.max(solver.getAbsoluteAccuracy(),
                                                FastMath.abs(solver.getRelativeAccuracy() * t0));
            final double tStart = t0 + 0.5 * epsilon;
            interpolator.setInterpolatedTime(tStart);
            g0 = handler.g(tStart, interpolator.getInterpolatedState());
        }
        g0Positive = g0 >= 0;

    } catch (MathUserException mue) {
        throw new EventException(mue);
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve() {
    final double min = getMin();
    final double max = getMax();

    verifyInterval(min, max);

    final double relativeAccuracy = getRelativeAccuracy();
    final double absoluteAccuracy = getAbsoluteAccuracy();
    final double functionValueAccuracy = getFunctionValueAccuracy();

    // x2 is the last root approximation
    // x is the new approximation and new x2 for next round
    // x0 < x1 < x2 does not hold here

    double x0 = min;
    double y0 = computeObjectiveValue(x0);
    if (FastMath.abs(y0) < functionValueAccuracy) {
        return x0;
    }
    double x1 = max;
    double y1 = computeObjectiveValue(x1);
    if (FastMath.abs(y1) < functionValueAccuracy) {
        return x1;
    }

    if(y0 * y1 > 0) {
        throw new NoBracketingException(x0, x1, y0, y1);
    }

    double x2 = 0.5 * (x0 + x1);
    double y2 = computeObjectiveValue(x2);

    double oldx = Double.POSITIVE_INFINITY;
    while (true) {
        // quadratic interpolation through x0, x1, x2
        final double q = (x2 - x1) / (x1 - x0);
        final double a = q * (y2 - (1 + q) * y1 + q * y0);
        final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
        final double c = (1 + q) * y2;
        final double delta = b * b - 4 * a * c;
        double x;
        final double denominator;
        if (delta >= 0.0) {
            // choose a denominator larger in magnitude
            double dplus = b + FastMath.sqrt(delta);
            double dminus = b - FastMath.sqrt(delta);
            denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
        } else {
            // take the modulus of (B +/- FastMath.sqrt(delta))
            denominator = FastMath.sqrt(b * b - delta);
        }
        if (denominator != 0) {
            x = x2 - 2.0 * c * (x2 - x1) / denominator;
            // perturb x if it exactly coincides with x1 or x2
            // the equality tests here are intentional
            while (x == x1 || x == x2) {
                x += absoluteAccuracy;
            }
        } else {
            // extremely rare case, get a random number to skip it
            x = min + FastMath.random() * (max - min);
            oldx = Double.POSITIVE_INFINITY;
        }
        final double y = computeObjectiveValue(x);

        // check for convergence
        final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
        if (FastMath.abs(x - oldx) <= tolerance ||
            FastMath.abs(y) <= functionValueAccuracy) {
            return x;
        }

        // prepare the next iteration
        x0 = x1;
        y0 = y1;
        x1 = x2;
        y1 = y2;
        x2 = x;
        y2 = y;
        oldx = x;
    }
}
 

/** Initialize the integration step.
 * @param forward forward integration indicator
 * @param order order of the method
 * @param scale scaling vector for the state vector (can be shorter than state vector)
 * @param t0 start time
 * @param y0 state vector at t0
 * @param yDot0 first time derivative of y0
 * @param y1 work array for a state vector
 * @param yDot1 work array for the first time derivative of y1
 * @return first integration step
 */
public double initializeStep(final boolean forward, final int order, final double[] scale,
                             final double t0, final double[] y0, final double[] yDot0,
                             final double[] y1, final double[] yDot1) {

  if (initialStep > 0) {
    // use the user provided value
    return forward ? initialStep : -initialStep;
  }

  // very rough first guess : h = 0.01 * ||y/scale|| / ||y'/scale||
  // this guess will be used to perform an Euler step
  double ratio;
  double yOnScale2 = 0;
  double yDotOnScale2 = 0;
  for (int j = 0; j < scale.length; ++j) {
    ratio         = y0[j] / scale[j];
    yOnScale2    += ratio * ratio;
    ratio         = yDot0[j] / scale[j];
    yDotOnScale2 += ratio * ratio;
  }

  double h = ((yOnScale2 < 1.0e-10) || (yDotOnScale2 < 1.0e-10)) ?
             1.0e-6 : (0.01 * FastMath.sqrt(yOnScale2 / yDotOnScale2));
  if (! forward) {
    h = -h;
  }

  // perform an Euler step using the preceding rough guess
  for (int j = 0; j < y0.length; ++j) {
    y1[j] = y0[j] + h * yDot0[j];
  }
  computeDerivatives(t0 + h, y1, yDot1);

  // estimate the second derivative of the solution
  double yDDotOnScale = 0;
  for (int j = 0; j < scale.length; ++j) {
    ratio         = (yDot1[j] - yDot0[j]) / scale[j];
    yDDotOnScale += ratio * ratio;
  }
  yDDotOnScale = FastMath.sqrt(yDDotOnScale) / h;

  // step size is computed such that
  // h^order * max (||y'/tol||, ||y''/tol||) = 0.01
  final double maxInv2 = FastMath.max(FastMath.sqrt(yDotOnScale2), yDDotOnScale);
  final double h1 = (maxInv2 < 1.0e-15) ?
                    FastMath.max(1.0e-6, 0.001 * FastMath.abs(h)) :
                    FastMath.pow(0.01 / maxInv2, 1.0 / order);
  h = FastMath.min(100.0 * FastMath.abs(h), h1);
  h = FastMath.max(h, 1.0e-12 * FastMath.abs(t0));  // avoids cancellation when computing t1 - t0
  if (h < getMinStep()) {
    h = getMinStep();
  }
  if (h > getMaxStep()) {
    h = getMaxStep();
  }
  if (! forward) {
    h = -h;
  }

  return h;

}
 

/** Get the L<sub>&infin;</sub> norm for the vector.
 * @return L<sub>&infin;</sub> norm for the vector
 */
public double getNormInf() {
  return FastMath.max(FastMath.max(FastMath.abs(x), FastMath.abs(y)), FastMath.abs(z));
}
 

/**
 * 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 new NullArgumentException(LocalizedFormats.FUNCTION);
    }
    if (maximumIterations <= 0)  {
        throw MathRuntimeException.createIllegalArgumentException(
              LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
    }
    if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
        throw MathRuntimeException.createIllegalArgumentException(
              LocalizedFormats.INVALID_BRACKETING_PARAMETERS,
              lowerBound, initial, upperBound);
    }
    double a = initial;
    double b = initial;
    double fa;
    double fb;
    int numIterations = 0 ;

    do {
        a = FastMath.max(a - 1.0, lowerBound);
        b = FastMath.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(
                  LocalizedFormats.FAILED_BRACKETING,
                  numIterations, maximumIterations, initial,
                  lowerBound, upperBound, a, b, fa, fb);
    }

    return new double[]{a, b};
}
 

/** Get the L<sub>&infin;</sub> norm for the vector.
 * @return L<sub>&infin;</sub> norm for the vector
 */
public double getNormInf() {
  return FastMath.max(FastMath.max(FastMath.abs(x), FastMath.abs(y)), FastMath.abs(z));
}
 

/**
 * Find a zero in the given interval.
 * @param f the function to solve
 * @param min the lower bound for the interval.
 * @param max the upper bound for the interval.
 * @return the value where the function is zero
 * @throws MaxIterationsExceededException  if the maximum iteration count is exceeded
 * @throws FunctionEvaluationException if an error occurs evaluating the
 * function
 * @throws IllegalArgumentException if min is not less than max or the
 * signs of the values of the function at the endpoints are not opposites
 */
public double solve(final UnivariateRealFunction f,
                    final double min, final double max)
    throws MaxIterationsExceededException, FunctionEvaluationException {

    clearResult();
    verifyInterval(min, max);

    // Index 0 is the old approximation for the root.
    // Index 1 is the last calculated approximation  for the root.
    // Index 2 is a bracket for the root with respect to x0.
    // OldDelta is the length of the bracketing interval of the last
    // iteration.
    double x0 = min;
    double x1 = max;
    double y0 = f.value(x0);
    double y1 = f.value(x1);

    // Verify bracketing
    if (y0 * y1 >= 0) {
        throw MathRuntimeException.createIllegalArgumentException(
              LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, y0, y1);
    }

    double x2 = x0;
    double y2 = y0;
    double oldDelta = x2 - x1;
    int i = 0;
    while (i < maximalIterationCount) {
        if (FastMath.abs(y2) < FastMath.abs(y1)) {
            x0 = x1;
            x1 = x2;
            x2 = x0;
            y0 = y1;
            y1 = y2;
            y2 = y0;
        }
        if (FastMath.abs(y1) <= functionValueAccuracy) {
            setResult(x1, i);
            return result;
        }
        if (FastMath.abs(oldDelta) <
            FastMath.max(relativeAccuracy * FastMath.abs(x1), absoluteAccuracy)) {
            setResult(x1, i);
            return result;
        }
        double delta;
        if (FastMath.abs(y1) > FastMath.abs(y0)) {
            // Function value increased in last iteration. Force bisection.
            delta = 0.5 * oldDelta;
        } else {
            delta = (x0 - x1) / (1 - y0 / y1);
            if (delta / oldDelta > 1) {
                // New approximation falls outside bracket.
                // Fall back to bisection.
                delta = 0.5 * oldDelta;
            }
        }
        x0 = x1;
        y0 = y1;
        x1 = x1 + delta;
        y1 = f.value(x1);
        if ((y1 > 0) == (y2 > 0)) {
            // New bracket is (x0,x1).
            x2 = x0;
            y2 = y0;
        }
        oldDelta = x2 - x1;
        i++;
    }
    throw new MaxIterationsExceededException(maximalIterationCount);
}
 

/** {@inheritDoc} */
@Override
public void setSubMatrix(final T[][] subMatrix, final int row, final int column) {
    // safety checks
    final int refLength = subMatrix[0].length;
    if (refLength == 0) {
        throw new NoDataException(LocalizedFormats.AT_LEAST_ONE_COLUMN);
    }
    final int endRow    = row + subMatrix.length - 1;
    final int endColumn = column + refLength - 1;
    checkSubMatrixIndex(row, endRow, column, endColumn);
    for (final T[] subRow : subMatrix) {
        if (subRow.length != refLength) {
            throw new DimensionMismatchException(refLength, subRow.length);
        }
    }

    // compute blocks bounds
    final int blockStartRow    = row / BLOCK_SIZE;
    final int blockEndRow      = (endRow + BLOCK_SIZE) / BLOCK_SIZE;
    final int blockStartColumn = column / BLOCK_SIZE;
    final int blockEndColumn   = (endColumn + BLOCK_SIZE) / BLOCK_SIZE;

    // perform copy block-wise, to ensure good cache behavior
    for (int iBlock = blockStartRow; iBlock < blockEndRow; ++iBlock) {
        final int iHeight  = blockHeight(iBlock);
        final int firstRow = iBlock * BLOCK_SIZE;
        final int iStart   = FastMath.max(row,    firstRow);
        final int iEnd     = FastMath.min(endRow + 1, firstRow + iHeight);

        for (int jBlock = blockStartColumn; jBlock < blockEndColumn; ++jBlock) {
            final int jWidth      = blockWidth(jBlock);
            final int firstColumn = jBlock * BLOCK_SIZE;
            final int jStart      = FastMath.max(column,    firstColumn);
            final int jEnd        = FastMath.min(endColumn + 1, firstColumn + jWidth);
            final int jLength     = jEnd - jStart;

            // handle one block, row by row
            final T[] block = blocks[iBlock * blockColumns + jBlock];
            for (int i = iStart; i < iEnd; ++i) {
                System.arraycopy(subMatrix[i - row], jStart - column,
                                 block, (i - firstRow) * jWidth + (jStart - firstColumn),
                                 jLength);
            }

        }
    }
}
 

/**
 * Check if the optimization algorithm has converged considering the
 * last two points.
 * This method may be called several time from the same algorithm
 * iteration with different points. This can be detected by checking the
 * iteration number at each call if needed. Each time this method is
 * called, the previous and current point correspond to points with the
 * same role at each iteration, so they can be compared. As an example,
 * simplex-based algorithms call this method for all points of the simplex,
 * not only for the best or worst ones.
 *
 * @param iteration Index of current iteration
 * @param previous Best point in the previous iteration.
 * @param current Best point in the current iteration.
 * @return {@code true} if the algorithm has converged.
 */
@Override
public boolean converged(final int iteration,
                         final RealPointValuePair previous,
                         final RealPointValuePair current) {
    final double p = previous.getValue();
    final double c = current.getValue();
    final double difference = FastMath.abs(p - c);
    final double size = FastMath.max(FastMath.abs(p), FastMath.abs(c));
    return difference <= size * getRelativeThreshold() ||
        difference <= getAbsoluteThreshold();
}
 

/**
 * {@inheritDoc}
 *
 * For population size <code>N</code>,
 * number of successes <code>m</code>, and
 * sample size <code>n</code>,
 * the lower bound of the support is
 * <code>max(0, n + m - N)</code>
 *
 * @return lower bound of the support
 */
@Override
public int getSupportLowerBound() {
    return FastMath.max(0,
            getSampleSize() + getNumberOfSuccesses() - getPopulationSize());
}