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

/** Filter the integration step.
 * @param h signed step
 * @param forward forward integration indicator
 * @param acceptSmall if true, steps smaller than the minimal value
 * are silently increased up to this value, if false such small
 * steps generate an exception
 * @return a bounded integration step (h if no bound is reach, or a bounded value)
 * @exception IntegratorException if the step is too small and acceptSmall is false
 */
protected double filterStep(final double h, final boolean forward, final boolean acceptSmall)
  throws IntegratorException {

    double filteredH = h;
    if (FastMath.abs(h) < minStep) {
        if (acceptSmall) {
            filteredH = forward ? minStep : -minStep;
        } else {
            throw new IntegratorException(
                    LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
                    minStep, FastMath.abs(h));
        }
    }

    if (filteredH > maxStep) {
        filteredH = maxStep;
    } else if (filteredH < -maxStep) {
        filteredH = -maxStep;
    }

    return filteredH;

}
 

@Test
public void backward() throws MathUserException, IntegratorException {

    TestProblem5 pb = new TestProblem5();
    double range = FastMath.abs(pb.getFinalTime() - pb.getInitialTime());

    FirstOrderIntegrator integ = new AdamsBashforthIntegrator(4, 0, range, 1.0e-12, 1.0e-12);
    TestProblemHandler handler = new TestProblemHandler(pb, integ);
    integ.addStepHandler(handler);
    integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
                    pb.getFinalTime(), new double[pb.getDimension()]);

    Assert.assertTrue(handler.getLastError() < 1.5e-8);
    Assert.assertTrue(handler.getMaximalValueError() < 1.5e-8);
    Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
    Assert.assertEquals("Adams-Bashforth", integ.getName());
}
 

/** 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
 * @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 solveUpperTriangularSystem( 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 = rows-1 ; i >-1 ; i-- ){
        double diag = rm.getEntry(i, i);
        if( FastMath.abs(diag) < MathUtils.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)  );
        }
    }
}
 

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

    min = 0; max = FastMath.PI; expected = 2;
    tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
    result = integrator.integrate(f, min, max);
    Assert.assertEquals(expected, result, tolerance);

    min = -FastMath.PI/3; max = 0; expected = -0.5;
    tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
    result = integrator.integrate(f, min, max);
    Assert.assertEquals(expected, result, tolerance);
}
 

/** Filter the integration step.
 * @param h signed step
 * @param forward forward integration indicator
 * @param acceptSmall if true, steps smaller than the minimal value
 * are silently increased up to this value, if false such small
 * steps generate an exception
 * @return a bounded integration step (h if no bound is reach, or a bounded value)
 * @exception IntegratorException if the step is too small and acceptSmall is false
 */
protected double filterStep(final double h, final boolean forward, final boolean acceptSmall)
  throws IntegratorException {

    double filteredH = h;
    if (FastMath.abs(h) < minStep) {
        if (acceptSmall) {
            filteredH = forward ? minStep : -minStep;
        } else {
            throw new IntegratorException(
                    LocalizedFormats.MINIMAL_STEPSIZE_REACHED_DURING_INTEGRATION,
                    minStep, FastMath.abs(h));
        }
    }

    if (filteredH > maxStep) {
        filteredH = maxStep;
    } else if (filteredH < -maxStep) {
        filteredH = -maxStep;
    }

    return filteredH;

}
 

/** {@inheritDoc} */
public double getL1Norm() {
    double norm = 0;
    Iterator<Entry> it = sparseIterator();
    Entry e;
    while (it.hasNext() && (e = it.next()) != null) {
        norm += FastMath.abs(e.getValue());
    }
    return norm;
}
 

/** {@inheritDoc} */
public double getL1Distance(double[] v) {
    checkVectorDimensions(v.length);
    double d = 0;
    Iterator<Entry> it = iterator();
    Entry e;
    while (it.hasNext() && (e = it.next()) != null) {
        d += FastMath.abs(e.getValue() - v[e.getIndex()]);
    }
    return d;
}
 

/** {@inheritDoc} */
public double distanceInf(Vector<Euclidean3D> v) {
    final Vector3D v3 = (Vector3D) v;
    final double dx = FastMath.abs(v3.x - x);
    final double dy = FastMath.abs(v3.y - y);
    final double dz = FastMath.abs(v3.z - z);
    return FastMath.max(FastMath.max(dx, dy), dz);
}
 

@Override
public double[] computeTheoreticalState(double t) {

  // solve Kepler's equation
  double E = t;
  double d = 0;
  double corr = 999.0;
  for (int i = 0; (i < 50) && (FastMath.abs(corr) > 1.0e-12); ++i) {
    double f2  = e * FastMath.sin(E);
    double f0  = d - f2;
    double f1  = 1 - e * FastMath.cos(E);
    double f12 = f1 + f1;
    corr  = f0 * f12 / (f1 * f12 - f0 * f2);
    d -= corr;
    E = t + d;
  }

  double cosE = FastMath.cos(E);
  double sinE = FastMath.sin(E);

  y[0] = cosE - e;
  y[1] = FastMath.sqrt(1 - e * e) * sinE;
  y[2] = -sinE / (1 - e * cosE);
  y[3] = FastMath.sqrt(1 - e * e) * cosE / (1 - e * cosE);

  return y;
}
 

/** {@inheritDoc} */
@Override
public double getL1Distance(double[] v)
    throws IllegalArgumentException {
    checkVectorDimensions(v.length);
    double sum = 0;
    for (int i = 0; i < data.length; ++i) {
        final double delta = data[i] - v[i];
        sum += FastMath.abs(delta);
    }
    return sum;
}
 

/** {@inheritDoc} */
@Override
public double getLInfDistance(double[] v) {
    checkVectorDimensions(v.length);
    double max = 0;
    for (int i = 0; i < v.length; i++) {
        double delta = FastMath.abs(getEntry(i) - v[i]);
        if (delta > max) {
            max = delta;
        }
    }
    return max;
}
 

/** {@inheritDoc} */
public double getNorm() {
    return FastMath.abs(x);
}
 

/**
 * Create a fraction given the double value and either the maximum error
 * allowed or the maximum number of denominator digits.
 * <p>
 *
 * NOTE: This constructor is called with EITHER - a valid epsilon value and
 * the maxDenominator set to Integer.MAX_VALUE (that way the maxDenominator
 * has no effect). OR - a valid maxDenominator value and the epsilon value
 * set to zero (that way epsilon only has effect if there is an exact match
 * before the maxDenominator value is reached).
 * </p>
 * <p>
 *
 * It has been done this way so that the same code can be (re)used for both
 * scenarios. However this could be confusing to users if it were part of
 * the public API and this constructor should therefore remain PRIVATE.
 * </p>
 *
 * See JIRA issue ticket MATH-181 for more details:
 *
 * https://issues.apache.org/jira/browse/MATH-181
 *
 * @param value
 *            the double value to convert to a fraction.
 * @param epsilon
 *            maximum error allowed. The resulting fraction is within
 *            <code>epsilon</code> of <code>value</code>, in absolute terms.
 * @param maxDenominator
 *            maximum denominator value allowed.
 * @param maxIterations
 *            maximum number of convergents.
 * @throws FractionConversionException
 *             if the continued fraction failed to converge.
 */
private BigFraction(final double value, final double epsilon,
                    final int maxDenominator, int maxIterations)
    throws FractionConversionException {
    long overflow = Integer.MAX_VALUE;
    double r0 = value;
    long a0 = (long) FastMath.floor(r0);
    if (a0 > overflow) {
        throw new FractionConversionException(value, a0, 1l);
    }

    // check for (almost) integer arguments, which should not go
    // to iterations.
    if (FastMath.abs(a0 - value) < epsilon) {
        numerator = BigInteger.valueOf(a0);
        denominator = BigInteger.ONE;
        return;
    }

    long p0 = 1;
    long q0 = 0;
    long p1 = a0;
    long q1 = 1;

    long p2 = 0;
    long q2 = 1;

    int n = 0;
    boolean stop = false;
    do {
        ++n;
        final double r1 = 1.0 / (r0 - a0);
        final long a1 = (long) FastMath.floor(r1);
        p2 = (a1 * p1) + p0;
        q2 = (a1 * q1) + q0;
        if ((p2 > overflow) || (q2 > overflow)) {
            throw new FractionConversionException(value, p2, q2);
        }

        final double convergent = (double) p2 / (double) q2;
        if ((n < maxIterations) &&
            (FastMath.abs(convergent - value) > epsilon) &&
            (q2 < maxDenominator)) {
            p0 = p1;
            p1 = p2;
            q0 = q1;
            q1 = q2;
            a0 = a1;
            r0 = r1;
        } else {
            stop = true;
        }
    } while (!stop);

    if (n >= maxIterations) {
        throw new FractionConversionException(value, maxIterations);
    }

    if (q2 < maxDenominator) {
        numerator   = BigInteger.valueOf(p2);
        denominator = BigInteger.valueOf(q2);
    } else {
        numerator   = BigInteger.valueOf(p1);
        denominator = BigInteger.valueOf(q1);
    }
}
 

/** Perfect orthogonality on a 3X3 matrix.
 * @param m initial matrix (not exactly orthogonal)
 * @param threshold convergence threshold for the iterative
 * orthogonality correction (convergence is reached when the
 * difference between two steps of the Frobenius norm of the
 * correction is below this threshold)
 * @return an orthogonal matrix close to m
 * @exception NotARotationMatrixException if the matrix cannot be
 * orthogonalized with the given threshold after 10 iterations
 */
private double[][] orthogonalizeMatrix(double[][] m, double threshold)
  throws NotARotationMatrixException {
  double[] m0 = m[0];
  double[] m1 = m[1];
  double[] m2 = m[2];
  double x00 = m0[0];
  double x01 = m0[1];
  double x02 = m0[2];
  double x10 = m1[0];
  double x11 = m1[1];
  double x12 = m1[2];
  double x20 = m2[0];
  double x21 = m2[1];
  double x22 = m2[2];
  double fn = 0;
  double fn1;

  double[][] o = new double[3][3];
  double[] o0 = o[0];
  double[] o1 = o[1];
  double[] o2 = o[2];

  // iterative correction: Xn+1 = Xn - 0.5 * (Xn.Mt.Xn - M)
  int i = 0;
  while (++i < 11) {

    // Mt.Xn
    double mx00 = m0[0] * x00 + m1[0] * x10 + m2[0] * x20;
    double mx10 = m0[1] * x00 + m1[1] * x10 + m2[1] * x20;
    double mx20 = m0[2] * x00 + m1[2] * x10 + m2[2] * x20;
    double mx01 = m0[0] * x01 + m1[0] * x11 + m2[0] * x21;
    double mx11 = m0[1] * x01 + m1[1] * x11 + m2[1] * x21;
    double mx21 = m0[2] * x01 + m1[2] * x11 + m2[2] * x21;
    double mx02 = m0[0] * x02 + m1[0] * x12 + m2[0] * x22;
    double mx12 = m0[1] * x02 + m1[1] * x12 + m2[1] * x22;
    double mx22 = m0[2] * x02 + m1[2] * x12 + m2[2] * x22;

    // Xn+1
    o0[0] = x00 - 0.5 * (x00 * mx00 + x01 * mx10 + x02 * mx20 - m0[0]);
    o0[1] = x01 - 0.5 * (x00 * mx01 + x01 * mx11 + x02 * mx21 - m0[1]);
    o0[2] = x02 - 0.5 * (x00 * mx02 + x01 * mx12 + x02 * mx22 - m0[2]);
    o1[0] = x10 - 0.5 * (x10 * mx00 + x11 * mx10 + x12 * mx20 - m1[0]);
    o1[1] = x11 - 0.5 * (x10 * mx01 + x11 * mx11 + x12 * mx21 - m1[1]);
    o1[2] = x12 - 0.5 * (x10 * mx02 + x11 * mx12 + x12 * mx22 - m1[2]);
    o2[0] = x20 - 0.5 * (x20 * mx00 + x21 * mx10 + x22 * mx20 - m2[0]);
    o2[1] = x21 - 0.5 * (x20 * mx01 + x21 * mx11 + x22 * mx21 - m2[1]);
    o2[2] = x22 - 0.5 * (x20 * mx02 + x21 * mx12 + x22 * mx22 - m2[2]);

    // correction on each elements
    double corr00 = o0[0] - m0[0];
    double corr01 = o0[1] - m0[1];
    double corr02 = o0[2] - m0[2];
    double corr10 = o1[0] - m1[0];
    double corr11 = o1[1] - m1[1];
    double corr12 = o1[2] - m1[2];
    double corr20 = o2[0] - m2[0];
    double corr21 = o2[1] - m2[1];
    double corr22 = o2[2] - m2[2];

    // Frobenius norm of the correction
    fn1 = corr00 * corr00 + corr01 * corr01 + corr02 * corr02 +
          corr10 * corr10 + corr11 * corr11 + corr12 * corr12 +
          corr20 * corr20 + corr21 * corr21 + corr22 * corr22;

    // convergence test
    if (FastMath.abs(fn1 - fn) <= threshold) {
        return o;
    }

    // prepare next iteration
    x00 = o0[0];
    x01 = o0[1];
    x02 = o0[2];
    x10 = o1[0];
    x11 = o1[1];
    x12 = o1[2];
    x20 = o2[0];
    x21 = o2[1];
    x22 = o2[2];
    fn  = fn1;

  }

  // the algorithm did not converge after 10 iterations
  throw new NotARotationMatrixException(
          LocalizedFormats.UNABLE_TO_ORTHOGONOLIZE_MATRIX,
          i - 1);
}
 

/**
 * Create a fraction given the double value and either the maximum error
 * allowed or the maximum number of denominator digits.
 * <p>
 *
 * NOTE: This constructor is called with EITHER - a valid epsilon value and
 * the maxDenominator set to Integer.MAX_VALUE (that way the maxDenominator
 * has no effect). OR - a valid maxDenominator value and the epsilon value
 * set to zero (that way epsilon only has effect if there is an exact match
 * before the maxDenominator value is reached).
 * </p>
 * <p>
 *
 * It has been done this way so that the same code can be (re)used for both
 * scenarios. However this could be confusing to users if it were part of
 * the public API and this constructor should therefore remain PRIVATE.
 * </p>
 *
 * See JIRA issue ticket MATH-181 for more details:
 *
 * https://issues.apache.org/jira/browse/MATH-181
 *
 * @param value
 *            the double value to convert to a fraction.
 * @param epsilon
 *            maximum error allowed. The resulting fraction is within
 *            <code>epsilon</code> of <code>value</code>, in absolute terms.
 * @param maxDenominator
 *            maximum denominator value allowed.
 * @param maxIterations
 *            maximum number of convergents.
 * @throws FractionConversionException
 *             if the continued fraction failed to converge.
 */
private BigFraction(final double value, final double epsilon,
                    final int maxDenominator, int maxIterations)
    throws FractionConversionException {
    long overflow = Integer.MAX_VALUE;
    double r0 = value;
    long a0 = (long) FastMath.floor(r0);
    if (a0 > overflow) {
        throw new FractionConversionException(value, a0, 1l);
    }

    // check for (almost) integer arguments, which should not go
    // to iterations.
    if (FastMath.abs(a0 - value) < epsilon) {
        numerator = BigInteger.valueOf(a0);
        denominator = BigInteger.ONE;
        return;
    }

    long p0 = 1;
    long q0 = 0;
    long p1 = a0;
    long q1 = 1;

    long p2 = 0;
    long q2 = 1;

    int n = 0;
    boolean stop = false;
    do {
        ++n;
        final double r1 = 1.0 / (r0 - a0);
        final long a1 = (long) FastMath.floor(r1);
        p2 = (a1 * p1) + p0;
        q2 = (a1 * q1) + q0;
        if ((p2 > overflow) || (q2 > overflow)) {
            throw new FractionConversionException(value, p2, q2);
        }

        final double convergent = (double) p2 / (double) q2;
        if ((n < maxIterations) &&
            (FastMath.abs(convergent - value) > epsilon) &&
            (q2 < maxDenominator)) {
            p0 = p1;
            p1 = p2;
            q0 = q1;
            q1 = q2;
            a0 = a1;
            r0 = r1;
        } else {
            stop = true;
        }
    } while (!stop);

    if (n >= maxIterations) {
        throw new FractionConversionException(value, maxIterations);
    }

    if (q2 < maxDenominator) {
        numerator   = BigInteger.valueOf(p2);
        denominator = BigInteger.valueOf(q2);
    } else {
        numerator   = BigInteger.valueOf(p1);
        denominator = BigInteger.valueOf(q1);
    }
}
 

/**
 * 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 points Points used for checking convergence. The list must
 * contain two elements:
 * <ul>
 *  <li>the previous best point,</li>
 *  <li>the current best point.</li>
 * </ul>
 * @return {@code true} if the algorithm has converged.
 * @throws DimensionMismatchException if the length of the {@code points}
 * list is not equal to 2.
 */
public boolean converged(final int iteration,
                         final RealPointValuePair ... points) {
    if (points.length != 2) {
        throw new DimensionMismatchException(points.length, 2);
    }

    final double[] p = points[0].getPoint();
    final double[] c = points[1].getPoint();
    for (int i = 0; i < p.length; ++i) {
        final double difference = FastMath.abs(p[i] - c[i]);
        final double size = FastMath.max(FastMath.abs(p[i]), FastMath.abs(c[i]));
        if (difference > size * getRelativeThreshold() &&
            difference > getAbsoluteThreshold()) {
            return false;
        }
    }
    return true;
}
 

/** Compute the distance between two vectors according to the L<sub>1</sub> norm.
 * <p>Calling this method is equivalent to calling:
 * <code>v1.subtract(v2).getNorm1()</code> except that no intermediate
 * vector is built</p>
 * @param v1 first vector
 * @param v2 second vector
 * @return the distance between v1 and v2 according to the L<sub>1</sub> norm
 */
public static double distance1(Vector3D v1, Vector3D v2) {
  final double dx = FastMath.abs(v2.x - v1.x);
  final double dy = FastMath.abs(v2.y - v1.y);
  final double dz = FastMath.abs(v2.z - v1.z);
  return dx + dy + dz;
}
 

/**
 * For this distribution, {@code X}, this method returns {@code P(X < x)}.
 * If {@code x}is more than 40 standard deviations from the mean, 0 or 1 is returned,
 * as in these cases the actual value is within {@code Double.MIN_VALUE} of 0 or 1.
 *
 * @param x Value at which the CDF is evaluated.
 * @return CDF evaluated at {@code x}.
 * @throws MathException if the algorithm fails to converge
 */
public double cumulativeProbability(double x) throws MathException {
    final double dev = x - mean;
    if (FastMath.abs(dev) > 40 * standardDeviation) {
        return dev < 0 ? 0.0d : 1.0d;
    }
    return 0.5 * (1 + Erf.erf(dev / (standardDeviation * FastMath.sqrt(2))));
}
 

/**
 * Returns the error function.
 *
 * <p>erf(x) = 2/&radic;&pi; <sub>0</sub>&int;<sup>x</sup> e<sup>-t<sup>2</sup></sup>dt </p>
 *
 * <p>This implementation computes erf(x) using the
 * {@link Gamma#regularizedGammaP(double, double, double, int) regularized gamma function},
 * following <a href="http://mathworld.wolfram.com/Erf.html"> Erf</a>, equation (3)</p>
 *
 * <p>The value returned is always between -1 and 1 (inclusive).
 * If {@code abs(x) > 40}, then {@code erf(x)} is indistinguishable from
 * either 1 or -1 as a double, so the appropriate extreme value is returned.
 * </p>
 *
 * @param x the value.
 * @return the error function erf(x)
 * @throws org.apache.commons.math.exception.MaxCountExceededException
 * if the algorithm fails to converge.
 * @see Gamma#regularizedGammaP(double, double, double, int)
 */
public static double erf(double x) {
    if (FastMath.abs(x) > 40) {
        return x > 0 ? 1 : -1;
    }
    double ret = Gamma.regularizedGammaP(0.5, x * x, 1.0e-15, 10000);
    if (x < 0) {
        ret = -ret;
    }
    return ret;
}
 

/**
 * Determine if this value is within epsilon of zero.
 *
 * @param value Value to test
 * @return {@code true} if this value is within epsilon to zero,
 * {@code false} otherwise.
 * @since 2.1
 */
protected boolean isDefaultValue(double value) {
    return FastMath.abs(value) < epsilon;
}