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

@Test
public void testSmallError() {
    Random randomizer = new Random(53882150042l);
    double maxError = 0;
    for (int degree = 0; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
        for (double x = -1.0; x < 1.0; x += 0.01) {
            fitter.addObservedPoint(1.0, x,
                                    p.value(x) + 0.1 * randomizer.nextGaussian());
        }

        final double[] init = new double[degree + 1];
        PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init));

        for (double x = -1.0; x < 1.0; x += 0.01) {
            double error = FastMath.abs(p.value(x) - fitted.value(x)) /
                          (1.0 + FastMath.abs(p.value(x)));
            maxError = FastMath.max(maxError, error);
            Assert.assertTrue(FastMath.abs(error) < 0.1);
        }
    }
    Assert.assertTrue(maxError > 0.01);
}
 

/**
 * Test of solver for the quintic function.
 */
@Test
public void testQuinticFunction() {
    UnivariateFunction f = new QuinticFunction();
    UnivariateSolver solver = new RiddersSolver();
    double min, max, expected, result, tolerance;

    min = -0.4; max = 0.2; 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 = 0.75; max = 1.5; expected = 1.0;
    tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
                FastMath.abs(expected * solver.getRelativeAccuracy()));
    result = solver.solve(100, f, min, max);
    Assert.assertEquals(expected, result, tolerance);

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

/** {@inheritDoc} */
@Override
protected double estimateError(final double[][] yDotK,
                               final double[] y0, final double[] y1,
                               final double h) {

  double error = 0;

  for (int j = 0; j < mainSetDimension; ++j) {
      final double errSum = E1 * yDotK[0][j] +  E3 * yDotK[2][j] +
                            E4 * yDotK[3][j] +  E5 * yDotK[4][j] +
                            E6 * yDotK[5][j] +  E7 * yDotK[6][j];

      final double yScale = FastMath.max(FastMath.abs(y0[j]), FastMath.abs(y1[j]));
      final double tol = (vecAbsoluteTolerance == null) ?
                         (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
                             (vecAbsoluteTolerance[j] + vecRelativeTolerance[j] * yScale);
      final double ratio  = h * errSum / tol;
      error += ratio * ratio;

  }

  return FastMath.sqrt(error / mainSetDimension);

}
 

/**
 * End visiting the Nordsieck vector.
 * <p>The correction is used to control stepsize. So its amplitude is
 * considered to be an error, which must be normalized according to
 * error control settings. If the normalized value is greater than 1,
 * the correction was too large and the step must be rejected.</p>
 * @return the normalized correction, if greater than 1, the step
 * must be rejected
 */
public double end() {

    double error = 0;
    for (int i = 0; i < after.length; ++i) {
        after[i] += previous[i] + scaled[i];
        if (i < mainSetDimension) {
            final double yScale = FastMath.max(FastMath.abs(previous[i]), FastMath.abs(after[i]));
            final double tol = (vecAbsoluteTolerance == null) ?
                               (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
                               (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);
            final double ratio  = (after[i] - before[i]) / tol;
            error += ratio * ratio;
        }
    }

    return FastMath.sqrt(error / mainSetDimension);

}
 

@Test
public void testLargeSample() {
    final Random randomizer = new Random(0x5551480dca5b369bl);
    double maxError = 0;
    for (int degree = 0; degree < 10; ++degree) {
        final PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
        final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(degree);
 
        final WeightedObservedPoints obs = new WeightedObservedPoints();
        for (int i = 0; i < 40000; ++i) {
            final double x = -1.0 + i / 20000.0;
            obs.add(1.0, x, p.value(x) + 0.1 * randomizer.nextGaussian());
        }

        final PolynomialFunction fitted = new PolynomialFunction(fitter.fit(obs.toList()));
        for (double x = -1.0; x < 1.0; x += 0.01) {
            final double error = FastMath.abs(p.value(x) - fitted.value(x)) /
                (1.0 + FastMath.abs(p.value(x)));
            maxError = FastMath.max(maxError, error);
            Assert.assertTrue(FastMath.abs(error) < 0.01);
        }
    }
    Assert.assertTrue(maxError > 0.001);
}
 

/**
 * End visiting the Nordsieck vector.
 * <p>The correction is used to control stepsize. So its amplitude is
 * considered to be an error, which must be normalized according to
 * error control settings. If the normalized value is greater than 1,
 * the correction was too large and the step must be rejected.</p>
 * @return the normalized correction, if greater than 1, the step
 * must be rejected
 */
public double end() {

    double error = 0;
    for (int i = 0; i < after.length; ++i) {
        after[i] += previous[i] + scaled[i];
        if (i < mainSetDimension) {
            final double yScale = FastMath.max(FastMath.abs(previous[i]), FastMath.abs(after[i]));
            final double tol = (vecAbsoluteTolerance == null) ?
                               (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
                               (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);
            final double ratio  = (after[i] - before[i]) / tol;
            error += ratio * ratio;
        }
    }

    return FastMath.sqrt(error / mainSetDimension);

}
 

@Test
public void testGetLInfNorm() {
    final double x = getPreferredEntryValue();
    final double[] data = new double[] { x, x, 1d, x, 2d, x, x, 3d, x };
    final RealVector v = create(data);
    final double actual = v.getLInfNorm();
    double expected = 0d;
    for (int i = 0; i < data.length; i++) {
        expected = FastMath.max(expected, FastMath.abs(data[i]));
    }
    Assert.assertEquals("", expected, actual, 0d);

}
 

/** {@inheritDoc} */
@Override
public double getLInfNorm() {
    double max = 0;
    for (double a : data) {
        max = FastMath.max(max, FastMath.abs(a));
    }
    return max;
}
 

/** {@inheritDoc} */
@Override
public double walkInOptimizedOrder(final RealMatrixChangingVisitor visitor,
                                   final int startRow, final int endRow,
                                   final int startColumn,
                                   final int endColumn)
    throws OutOfRangeException, NumberIsTooSmallException {
    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 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)
    throws OutOfRangeException, NumberIsTooSmallException {
    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();
}
 

/**
 * Checks whether a matrix is symmetric, within a given relative tolerance.
 *
 * @param matrix Matrix to check.
 * @param relativeTolerance Tolerance of the symmetry check.
 * @param raiseException If {@code true}, an exception will be raised if
 * the matrix is not symmetric.
 * @return {@code true} if {@code matrix} is symmetric.
 * @throws NonSquareMatrixException if the matrix is not square.
 * @throws NonSymmetricMatrixException if the matrix is not symmetric.
 */
private static boolean isSymmetricInternal(RealMatrix matrix,
                                           double relativeTolerance,
                                           boolean raiseException) {
    final int rows = matrix.getRowDimension();
    if (rows != matrix.getColumnDimension()) {
        if (raiseException) {
            throw new NonSquareMatrixException(rows, matrix.getColumnDimension());
        } else {
            return false;
        }
    }
    for (int i = 0; i < rows; i++) {
        for (int j = i + 1; j < rows; 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)) * relativeTolerance) {
                if (raiseException) {
                    throw new NonSymmetricMatrixException(i, j, relativeTolerance);
                } else {
                    return false;
                }
            }
        }
    }
    return true;
}
 

@Test
public void testStepSizeUnstability() {
    UnivariateDifferentiableFunction quintic = new QuinticFunction();
    UnivariateDifferentiableFunction goodStep =
            new FiniteDifferencesDifferentiator(7, 0.25).differentiate(quintic);
    UnivariateDifferentiableFunction badStep =
            new FiniteDifferencesDifferentiator(7, 1.0e-6).differentiate(quintic);
    double[] maxErrorGood = new double[7];
    double[] maxErrorBad  = new double[7];
    for (double x = -10; x < 10; x += 0.1) {
        DerivativeStructure dsX  = new DerivativeStructure(1, 6, 0, x);
        DerivativeStructure yRef  = quintic.value(dsX);
        DerivativeStructure yGood = goodStep.value(dsX);
        DerivativeStructure yBad  = badStep.value(dsX);
        for (int order = 0; order <= 6; ++order) {
            maxErrorGood[order] = FastMath.max(maxErrorGood[order],
                                               FastMath.abs(yRef.getPartialDerivative(order) -
                                                            yGood.getPartialDerivative(order)));
            maxErrorBad[order]  = FastMath.max(maxErrorBad[order],
                                               FastMath.abs(yRef.getPartialDerivative(order) -
                                                            yBad.getPartialDerivative(order)));
        }
    }

    // the 0.25 step size is good for finite differences in the quintic on this abscissa range for 7 points
    // the errors are fair
    final double[] expectedGood = new double[] {
        7.276e-12, 7.276e-11, 9.968e-10, 3.092e-9, 5.432e-8, 8.196e-8, 1.818e-6
    };

    // the 1.0e-6 step size is far too small for finite differences in the quintic on this abscissa range for 7 points
    // the errors are huge!
    final double[] expectedBad = new double[] {
        1.792e-22, 6.926e-5, 56.25, 1.783e8, 2.468e14, 3.056e20, 5.857e26            
    };

    for (int i = 0; i < maxErrorGood.length; ++i) {
        Assert.assertEquals(expectedGood[i], maxErrorGood[i], 0.01 * expectedGood[i]);
        Assert.assertEquals(expectedBad[i],  maxErrorBad[i],  0.01 * expectedBad[i]);
    }

}
 

/**
 * 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 UnivariatePointValuePair previous,
                         final UnivariatePointValuePair current) {
    if (maxIterationCount != ITERATION_CHECK_DISABLED && iteration >= maxIterationCount) {
        return true;
    }

    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();
}
 

/** Force a root found by a non-bracketing solver to lie on a specified side,
 * as if the solver was a bracketing one.
 * @param maxEval maximal number of new evaluations of the function
 * (evaluations already done for finding the root should have already been subtracted
 * from this number)
 * @param f function to solve
 * @param bracketing bracketing solver to use for shifting the root
 * @param baseRoot original root found by a previous non-bracketing solver
 * @param min minimal bound of the search interval
 * @param max maximal bound of the search interval
 * @param allowedSolution the kind of solutions that the root-finding algorithm may
 * accept as solutions.
 * @return a root approximation, on the specified side of the exact root
 * @throws NoBracketingException if the function has the same sign at the
 * endpoints.
 */
public static double forceSide(final int maxEval, final UnivariateFunction f,
                               final BracketedUnivariateSolver<UnivariateFunction> bracketing,
                               final double baseRoot, final double min, final double max,
                               final AllowedSolution allowedSolution)
    throws NoBracketingException {

    if (allowedSolution == AllowedSolution.ANY_SIDE) {
        // no further bracketing required
        return baseRoot;
    }

    // find a very small interval bracketing the root
    final double step = FastMath.max(bracketing.getAbsoluteAccuracy(),
                                     FastMath.abs(baseRoot * bracketing.getRelativeAccuracy()));
    double xLo        = FastMath.max(min, baseRoot - step);
    double fLo        = f.value(xLo);
    double xHi        = FastMath.min(max, baseRoot + step);
    double fHi        = f.value(xHi);
    int remainingEval = maxEval - 2;
    while (remainingEval > 0) {

        if ((fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0)) {
            // compute the root on the selected side
            return bracketing.solve(remainingEval, f, xLo, xHi, baseRoot, allowedSolution);
        }

        // try increasing the interval
        boolean changeLo = false;
        boolean changeHi = false;
        if (fLo < fHi) {
            // increasing function
            if (fLo >= 0) {
                changeLo = true;
            } else {
                changeHi = true;
            }
        } else if (fLo > fHi) {
            // decreasing function
            if (fLo <= 0) {
                changeLo = true;
            } else {
                changeHi = true;
            }
        } else {
            // unknown variation
            changeLo = true;
            changeHi = true;
        }

        // update the lower bound
        if (changeLo) {
            xLo = FastMath.max(min, xLo - step);
            fLo  = f.value(xLo);
            remainingEval--;
        }

        // update the higher bound
        if (changeHi) {
            xHi = FastMath.min(max, xHi + step);
            fHi  = f.value(xHi);
            remainingEval--;
        }

    }

    throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
                                    xLo, xHi, fLo, fHi,
                                    maxEval - remainingEval, maxEval, baseRoot,
                                    min, max);

}
 

/** {@inheritDoc} */
@Override
protected final double doSolve()
    throws TooManyEvaluationsException,
           NoBracketingException {
    // Get initial solution
    double x0 = getMin();
    double x1 = getMax();
    double f0 = computeObjectiveValue(x0);
    double f1 = computeObjectiveValue(x1);

    // If one of the bounds is the exact root, return it. Since these are
    // not under-approximations or over-approximations, we can return them
    // regardless of the allowed solutions.
    if (f0 == 0.0) {
        return x0;
    }
    if (f1 == 0.0) {
        return x1;
    }

    // Verify bracketing of initial solution.
    verifyBracketing(x0, x1);

    // Get accuracies.
    final double ftol = getFunctionValueAccuracy();
    final double atol = getAbsoluteAccuracy();
    final double rtol = getRelativeAccuracy();

    // Keep finding better approximations.
    while (true) {
        // Calculate the next approximation.
        final double x = x1 - ((f1 * (x1 - x0)) / (f1 - f0));
        final double fx = computeObjectiveValue(x);

        // If the new approximation is the exact root, return it. Since
        // this is not an under-approximation or an over-approximation,
        // we can return it regardless of the allowed solutions.
        if (fx == 0.0) {
            return x;
        }

        // Update the bounds with the new approximation.
        x0 = x1;
        f0 = f1;
        x1 = x;
        f1 = fx;

        // If the function value of the last approximation is too small,
        // given the function value accuracy, then we can't get closer to
        // the root than we already are.
        if (FastMath.abs(f1) <= ftol) {
            return x1;
        }

        // If the current interval is within the given accuracies, we
        // are satisfied with the current approximation.
        if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1), atol)) {
            return x1;
        }
    }
}
 

public static Vector2D max(Vector2D a, Vector2D b) {
	return new Vector2D(FastMath.max(a.getX(), b.getX()), FastMath.max(a.getY(), b.getY()));
}
 

/** {@inheritDoc} */
public double getNormInf() {
    return FastMath.max(FastMath.max(FastMath.abs(x), FastMath.abs(y)), FastMath.abs(z));
}
 

/**
 * Distance between two vectors.
 * <p>This method computes the distance consistent with
 * L<sub>&infin;</sub> norm, i.e. the max of the absolute values of
 * element differences.</p>
 *
 * @param v Vector to which distance is requested.
 * @return the distance between two vectors.
 * @throws DimensionMismatchException if {@code v} is not the same size as
 * {@code this} vector.
 * @see #getDistance(RealVector)
 * @see #getL1Distance(RealVector)
 * @see #getLInfNorm()
 */
public double getLInfDistance(RealVector v)
    throws DimensionMismatchException {
    checkVectorDimensions(v);
    double d = 0;
    Iterator<Entry> it = iterator();
    while (it.hasNext()) {
        final Entry e = it.next();
        d = FastMath.max(FastMath.abs(e.getValue() - v.getEntry(e.getIndex())), d);
    }
    return d;
}
 

/**
 * Check whether the given complex root is actually a real zero
 * in the given interval, within the solver tolerance level.
 *
 * @param min Lower bound for the interval.
 * @param max Upper bound for the interval.
 * @param z Complex root.
 * @return {@code true} if z is a real zero.
 */
public boolean isRoot(double min, double max, Complex z) {
    if (isSequence(min, z.getReal(), max)) {
        double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy());
        return (FastMath.abs(z.getImaginary()) <= tolerance) ||
             (z.abs() <= getFunctionValueAccuracy());
    }
    return false;
}