Java Code Examples for org.apache.commons.math3.util.Precision#SAFE_MIN

/** test eigenvalues for a big matrix. */
@Test
public void testBigMatrix() {
    Random r = new Random(17748333525117l);
    double[] bigValues = new double[200];
    for (int i = 0; i < bigValues.length; ++i) {
        bigValues[i] = 2 * r.nextDouble() - 1;
    }
    Arrays.sort(bigValues);
    EigenDecomposition ed;
    ed = new EigenDecomposition(createTestMatrix(r, bigValues),
                                    Precision.SAFE_MIN);
    double[] eigenValues = ed.getRealEigenvalues();
    Assert.assertEquals(bigValues.length, eigenValues.length);
    for (int i = 0; i < bigValues.length; ++i) {
        Assert.assertEquals(bigValues[bigValues.length - i - 1], eigenValues[i], 2.0e-14);
    }
}
 

/** {@inheritDoc} */
@Override
protected void computeGeometricalProperties() {
    if (getTree(false).getCut() == null) {
        setBarycenter(Vector1D.NaN);
        setSize(((Boolean) getTree(false).getAttribute()) ? Double.POSITIVE_INFINITY : 0);
    } else {
        double size = 0.0;
        double sum = 0.0;
        for (final Interval interval : asList()) {
            size += interval.getSize();
            sum  += interval.getSize() * interval.getBarycenter();
        }
        setSize(size);
        if (Double.isInfinite(size)) {
            setBarycenter(Vector1D.NaN);
        } else if (size >= Precision.SAFE_MIN) {
            setBarycenter(new Vector1D(sum / size));
        } else {
            setBarycenter(((OrientedPoint) getTree(false).getCut().getHyperplane()).getLocation());
        }
    }
}
 

/**
 * @param xDegree Degree of the polynomial fitting functions along the
 * x-dimension.
 * @param yDegree Degree of the polynomial fitting functions along the
 * y-dimension.
 * @exception NotPositiveException if degrees are not positive
 */
public SmoothingPolynomialBicubicSplineInterpolator(int xDegree, int yDegree)
    throws NotPositiveException {
    if (xDegree < 0) {
        throw new NotPositiveException(xDegree);
    }
    if (yDegree < 0) {
        throw new NotPositiveException(yDegree);
    }
    this.xDegree = xDegree;
    this.yDegree = yDegree;

    final double safeFactor = 1e2;
    final SimpleVectorValueChecker checker
        = new SimpleVectorValueChecker(safeFactor * Precision.EPSILON,
                                       safeFactor * Precision.SAFE_MIN);
    xFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker));
    yFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker));
}
 

/** Compute the shortest distance between the instance and another line.
 * @param line line to check against the instance
 * @return shortest distance between the instance and the line
 */
public double distance(final Line line) {

    final Vector3D normal = Vector3D.crossProduct(direction, line.direction);
    final double n = normal.getNorm();
    if (n < Precision.SAFE_MIN) {
        // lines are parallel
        return distance(line.zero);
    }

    // signed separation of the two parallel planes that contains the lines
    final double offset = line.zero.subtract(zero).dotProduct(normal) / n;

    return FastMath.abs(offset);

}
 

/** Compute the shortest distance between the instance and another line.
 * @param line line to check against the instance
 * @return shortest distance between the instance and the line
 */
public double distance(final Line line) {

    final Vector3D normal = Vector3D.crossProduct(direction, line.direction);
    final double n = normal.getNorm();
    if (n < Precision.SAFE_MIN) {
        // lines are parallel
        return distance(line.zero);
    }

    // signed separation of the two parallel planes that contains the lines
    final double offset = line.zero.subtract(zero).dotProduct(normal) / n;

    return FastMath.abs(offset);

}
 

/** Compute the shortest distance between the instance and another line.
 * @param line line to check against the instance
 * @return shortest distance between the instance and the line
 */
public double distance(final Line line) {

    final Vector3D normal = Vector3D.crossProduct(direction, line.direction);
    final double n = normal.getNorm();
    if (n < Precision.SAFE_MIN) {
        // lines are parallel
        return distance(line.zero);
    }

    // signed separation of the two parallel planes that contains the lines
    final double offset = line.zero.subtract(zero).dotProduct(normal) / n;

    return FastMath.abs(offset);

}
 

/** Solver a  system composed  of an Upper Triangular Matrix
 * {@link RealMatrix}.
 * <p>
 * This method is called to solve systems of equations which are
 * of the lower triangular form. The matrix {@link RealMatrix}
 * is assumed, though not checked, to be in upper triangular form.
 * The vector {@link RealVector} is overwritten with the solution.
 * The matrix is checked that it is square and its dimensions match
 * the length of the vector.
 * </p>
 * @param rm RealMatrix which is upper triangular
 * @param b  RealVector this is overwritten
 * @throws DimensionMismatchException if the matrix and vector are not
 * conformable
 * @throws NonSquareMatrixException if the matrix {@code rm} is not
 * square
 * @throws MathArithmeticException if the absolute value of one of the diagonal
 * coefficient of {@code rm} is lower than {@link Precision#SAFE_MIN}
 */
public static void solveUpperTriangularSystem(RealMatrix rm, RealVector b)
    throws DimensionMismatchException, MathArithmeticException,
    NonSquareMatrixException {
    if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
        throw new DimensionMismatchException(
                (rm == null) ? 0 : rm.getRowDimension(),
                (b == null) ? 0 : b.getDimension());
    }
    if( rm.getColumnDimension() != rm.getRowDimension() ){
        throw new NonSquareMatrixException(rm.getRowDimension(),
                                           rm.getColumnDimension());
    }
    int rows = rm.getRowDimension();
    for( int i = rows-1 ; i >-1 ; i-- ){
        double diag = rm.getEntry(i, i);
        if( FastMath.abs(diag) < Precision.SAFE_MIN ){
            throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
        }
        double bi = b.getEntry(i)/diag;
        b.setEntry(i,  bi );
        for( int j = i-1; j>-1; j-- ){
            b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
        }
    }
}
 

/**Solve  a  system of composed of a Lower Triangular Matrix
 * {@link RealMatrix}.
 * <p>
 * This method is called to solve systems of equations which are
 * of the lower triangular form. The matrix {@link RealMatrix}
 * is assumed, though not checked, to be in lower triangular form.
 * The vector {@link RealVector} is overwritten with the solution.
 * The matrix is checked that it is square and its dimensions match
 * the length of the vector.
 * </p>
 * @param rm RealMatrix which is lower triangular
 * @param b  RealVector this is overwritten
 * @exception IllegalArgumentException if the matrix and vector are not conformable
 * @exception ArithmeticException there is a zero or near zero on the diagonal of rm
 */
public static void solveLowerTriangularSystem( RealMatrix rm, RealVector b){
    if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
        throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
                (rm == null) ? 0 : rm.getRowDimension(),
                (b == null) ? 0 : b.getDimension());
    }
    if( rm.getColumnDimension() != rm.getRowDimension() ){
        throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
                rm.getRowDimension(),rm.getRowDimension(),
                rm.getRowDimension(),rm.getColumnDimension());
    }
    int rows = rm.getRowDimension();
    for( int i = 0 ; i < rows ; i++ ){
        double diag = rm.getEntry(i, i);
        if( FastMath.abs(diag) < Precision.SAFE_MIN ){
            throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
        }
        double bi = b.getEntry(i)/diag;
        b.setEntry(i,  bi );
        for( int j = i+1; j< rows; j++ ){
            b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
        }
    }
}
 

/** Compute the shortest distance between the instance and another line.
 * @param line line to check against the instance
 * @return shortest distance between the instance and the line
 */
public double distance(final Line line) {

    final Vector3D normal = Vector3D.crossProduct(direction, line.direction);
    final double n = normal.getNorm();
    if (n < Precision.SAFE_MIN) {
        // lines are parallel
        return distance(line.zero);
    }

    // signed separation of the two parallel planes that contains the lines
    final double offset = line.zero.subtract(zero).dotProduct(normal) / n;

    return FastMath.abs(offset);

}
 

/** Compute the shortest distance between the instance and another line.
 * @param line line to check against the instance
 * @return shortest distance between the instance and the line
 */
public double distance(final Line line) {

    final Vector3D normal = Vector3D.crossProduct(direction, line.direction);
    final double n = normal.getNorm();
    if (n < Precision.SAFE_MIN) {
        // lines are parallel
        return distance(line.zero);
    }

    // signed separation of the two parallel planes that contains the lines
    final double offset = line.zero.subtract(zero).dotProduct(normal) / n;

    return FastMath.abs(offset);

}
 

/** Compute the shortest distance between the instance and another line.
 * @param line line to check against the instance
 * @return shortest distance between the instance and the line
 */
public double distance(final Line line) {

    final Vector3D normal = Vector3D.crossProduct(direction, line.direction);
    final double n = normal.getNorm();
    if (n < Precision.SAFE_MIN) {
        // lines are parallel
        return distance(line.zero);
    }

    // signed separation of the two parallel planes that contains the lines
    final double offset = line.zero.subtract(zero).dotProduct(normal) / n;

    return FastMath.abs(offset);

}
 

/**
 * Computes the normalized quaternion (the versor of the instance).
 * The norm of the quaternion must not be zero.
 *
 * @return a normalized quaternion.
 * @throws ZeroException if the norm of the quaternion is zero.
 */
public Quaternion normalize() {
    final double norm = getNorm();

    if (norm < Precision.SAFE_MIN) {
        throw new ZeroException(LocalizedFormats.NORM, norm);
    }

    return new Quaternion(q0 / norm,
                          q1 / norm,
                          q2 / norm,
                          q3 / norm);
}
 

/**
 * Matrix with eigenvalues {8, -1, -1}
 */
@Test
public void testRepeatedEigenvalue() {
    RealMatrix repeated = MatrixUtils.createRealMatrix(new double[][] {
            {3,  2,  4},
            {2,  0,  2},
            {4,  2,  3}
    });
    EigenDecomposition ed;
    ed = new EigenDecomposition(repeated, Precision.SAFE_MIN);
    checkEigenValues((new double[] {8, -1, -1}), ed, 1E-12);
    checkEigenVector((new double[] {2, 1, 2}), ed, 1E-12);
}
 

/** test diagonal matrix */
@Test
public void testDiagonal() {
    double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 };
    RealMatrix m = createDiagonalMatrix(diagonal, diagonal.length, diagonal.length);
    EigenDecomposition ed;
    ed = new EigenDecomposition(m, Precision.SAFE_MIN);
    Assert.assertEquals(diagonal[0], ed.getRealEigenvalue(3), 2.0e-15);
    Assert.assertEquals(diagonal[1], ed.getRealEigenvalue(2), 2.0e-15);
    Assert.assertEquals(diagonal[2], ed.getRealEigenvalue(1), 2.0e-15);
    Assert.assertEquals(diagonal[3], ed.getRealEigenvalue(0), 2.0e-15);
}
 

/**
 * Returns the inverse of this instance.
 * The norm of the quaternion must not be zero.
 *
 * @return the inverse.
 * @throws ZeroException if the norm (squared) of the quaternion is zero.
 */
public Quaternion getInverse() {
    final double squareNorm = q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3;
    if (squareNorm < Precision.SAFE_MIN) {
        throw new ZeroException(LocalizedFormats.NORM, squareNorm);
    }

    return new Quaternion(q0 / squareNorm,
                          -q1 / squareNorm,
                          -q2 / squareNorm,
                          -q3 / squareNorm);
}
 

/**
 * Build an optimizer for least squares problems with default values
 * for some of the tuning parameters (see the {@link
 * #LevenbergMarquardtOptimizer(double,double,double,double,double)
 * other contructor}.
 * The default values for the algorithm settings are:
 * <ul>
 *  <li>Initial step bound factor}: 100</li>
 *  <li>QR ranking threshold}: {@link Precision#SAFE_MIN}</li>
 * </ul>
 *
 * @param costRelativeTolerance Desired relative error in the sum of
 * squares.
 * @param parRelativeTolerance Desired relative error in the approximate
 * solution parameters.
 * @param orthoTolerance Desired max cosine on the orthogonality between
 * the function vector and the columns of the Jacobian.
 */
public LevenbergMarquardtOptimizer(double costRelativeTolerance,
                                   double parRelativeTolerance,
                                   double orthoTolerance) {
    this(100,
         costRelativeTolerance, parRelativeTolerance, orthoTolerance,
         Precision.SAFE_MIN);
}
 

/**
 * Constructor that allows the specification of a custom convergence
 * checker.
 * Note that all the usual convergence checks will be <em>disabled</em>.
 * The default values for the algorithm settings are:
 * <ul>
 *  <li>Initial step bound factor: 100</li>
 *  <li>Cost relative tolerance: 1e-10</li>
 *  <li>Parameters relative tolerance: 1e-10</li>
 *  <li>Orthogonality tolerance: 1e-10</li>
 *  <li>QR ranking threshold: {@link Precision#SAFE_MIN}</li>
 * </ul>
 *
 * @param checker Convergence checker.
 */
public LevenbergMarquardtOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
    this(100, checker, 1e-10, 1e-10, 1e-10, Precision.SAFE_MIN);
}
 

/**
 * Constructor that allows the specification of a custom convergence
 * checker.
 * Note that all the usual convergence checks will be <em>disabled</em>.
 * The default values for the algorithm settings are:
 * <ul>
 *  <li>Initial step bound factor: 100</li>
 *  <li>Cost relative tolerance: 1e-10</li>
 *  <li>Parameters relative tolerance: 1e-10</li>
 *  <li>Orthogonality tolerance: 1e-10</li>
 *  <li>QR ranking threshold: {@link Precision#SAFE_MIN}</li>
 * </ul>
 *
 * @param checker Convergence checker.
 */
public LevenbergMarquardtOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
    this(100, checker, 1e-10, 1e-10, 1e-10, Precision.SAFE_MIN);
}
 

/**
 * Constructor that allows the specification of a custom convergence
 * checker.
 * Note that all the usual convergence checks will be <em>disabled</em>.
 * The default values for the algorithm settings are:
 * <ul>
 *  <li>Initial step bound factor: 100</li>
 *  <li>Cost relative tolerance: 1e-10</li>
 *  <li>Parameters relative tolerance: 1e-10</li>
 *  <li>Orthogonality tolerance: 1e-10</li>
 *  <li>QR ranking threshold: {@link Precision#SAFE_MIN}</li>
 * </ul>
 *
 * @param checker Convergence checker.
 */
public LevenbergMarquardtOptimizer(ConvergenceChecker<PointVectorValuePair> checker) {
    this(100, checker, 1e-10, 1e-10, 1e-10, Precision.SAFE_MIN);
}
 

/**
 * Build an optimizer for least squares problems with default values
 * for all the tuning parameters (see the {@link
 * #LevenbergMarquardtOptimizer(double,double,double,double,double)
 * other contructor}.
 * The default values for the algorithm settings are:
 * <ul>
 *  <li>Initial step bound factor: 100</li>
 *  <li>Cost relative tolerance: 1e-10</li>
 *  <li>Parameters relative tolerance: 1e-10</li>
 *  <li>Orthogonality tolerance: 1e-10</li>
 *  <li>QR ranking threshold: {@link Precision#SAFE_MIN}</li>
 * </ul>
 */
public LevenbergMarquardtOptimizer() {
    this(100, 1e-10, 1e-10, 1e-10, Precision.SAFE_MIN);
}