org.apache.commons.math.util.FastMath Java Examples

/**
 * Test of integrator for the sine function.
 */
@Test
public void testSinFunction() throws MathException {
    UnivariateRealFunction f = new SinFunction();
    UnivariateRealIntegrator integrator = new RombergIntegrator();
    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);
}
 

public void handleStep(StepInterpolator interpolator,
                       boolean isLast) {

  double step = FastMath.abs(interpolator.getCurrentTime()
                         - interpolator.getPreviousTime());
  if (firstTime) {
    minStep   = FastMath.abs(step);
    maxStep   = minStep;
    firstTime = false;
  } else {
    if (step < minStep) {
      minStep = step;
    }
    if (step > maxStep) {
      maxStep = step;
    }
  }

  if (isLast) {
    Assert.assertTrue(minStep < (1.0 / 100.0));
    Assert.assertTrue(maxStep > (1.0 / 2.0));
  }
}
 

@Test
public void testMath226() {
    double[] mean = { 1, 1, 10, 1 };
    double[][] cov = {
            { 1, 3, 2, 6 },
            { 3, 13, 16, 2 },
            { 2, 16, 38, -1 },
            { 6, 2, -1, 197 }
    };
    RealMatrix covRM = MatrixUtils.createRealMatrix(cov);
    JDKRandomGenerator jg = new JDKRandomGenerator();
    jg.setSeed(5322145245211l);
    NormalizedRandomGenerator rg = new GaussianRandomGenerator(jg);
    CorrelatedRandomVectorGenerator sg =
        new CorrelatedRandomVectorGenerator(mean, covRM, 0.00001, rg);

    for (int i = 0; i < 10; i++) {
        double[] generated = sg.nextVector();
        Assert.assertTrue(FastMath.abs(generated[0] - 1) > 0.1);
    }

}
 

@Test
public void testNonInvertible() {

    LinearProblem problem = new LinearProblem(new double[][] {
            {  1, 2, -3 },
            {  2, 1,  3 },
            { -3, 0, -9 }
    }, new double[] { 1, 1, 1 });

    LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
    optimizer.optimize(100, problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 });
    Assert.assertTrue(FastMath.sqrt(problem.target.length) * optimizer.getRMS() > 0.6);
    try {
        optimizer.getCovariances();
        Assert.fail("an exception should have been thrown");
    } catch (SingularMatrixException ee) {
        // expected behavior
    }
}
 

/** Compute the shortest distance between the instance and another line.
 * @param line line to check agains 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);
    if (normal.getNorm() < 1.0e-10) {
        // lines are parallel
        return distance(line.zero);
    }

    // separating middle plane
    final Plane middle = new Plane(new Vector3D(0.5, zero, 0.5, line.zero), normal);

    // the lines are at the same distance on either side of the plane
    return 2 * FastMath.abs(middle.getOffset(zero));

}
 

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

  TestProblemAbstract pb = new TestProblem5();
  double minStep = 1.25;
  double maxStep = FastMath.abs(pb.getFinalTime() - pb.getInitialTime());
  double scalAbsoluteTolerance = 6.0e-4;
  double scalRelativeTolerance = 6.0e-4;

  AdaptiveStepsizeIntegrator integ =
    new DormandPrince54Integrator(minStep, maxStep,
                                  scalAbsoluteTolerance,
                                  scalRelativeTolerance);

  DP54SmallLastHandler handler = new DP54SmallLastHandler(minStep);
  integ.addStepHandler(handler);
  integ.setInitialStepSize(1.7);
  integ.integrate(pb,
                  pb.getInitialTime(), pb.getInitialState(),
                  pb.getFinalTime(), new double[pb.getDimension()]);
  Assert.assertTrue(handler.wasLastSeen());
  Assert.assertEquals("Dormand-Prince 5(4)", integ.getName());

}
 

public void testDensity() {
    ExponentialDistribution d1 = new ExponentialDistributionImpl(1);
    assertEquals(0.0, d1.density(-1e-9));
    assertEquals(1.0, d1.density(0.0));
    assertEquals(0.0, d1.density(1000.0));
    assertEquals(FastMath.exp(-1), d1.density(1.0));
    assertEquals(FastMath.exp(-2), d1.density(2.0));

    ExponentialDistribution d2 = new ExponentialDistributionImpl(3);
    assertEquals(1/3.0, d2.density(0.0));
    // computed using  print(dexp(1, rate=1/3), digits=10) in R 2.5
    assertEquals(0.2388437702, d2.density(1.0), 1e-8);

    // computed using  print(dexp(2, rate=1/3), digits=10) in R 2.5
    assertEquals(0.1711390397, d2.density(2.0), 1e-8);
}
 

/**
 * 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
public void testSomeValues() {
    final double k = 4;
    final double m = 5;
    final double b = 2;
    final double q = 3;
    final double a = -1;
    final double n = 2;

    final UnivariateRealFunction f = new Logistic(k, m, b, q, a, n);

    double x;
    x = m;
    Assert.assertEquals("x=" + x, a + (k - a) / FastMath.sqrt(1 + q), f.value(x), EPS);

    x = Double.NEGATIVE_INFINITY;
    Assert.assertEquals("x=" + x, a, f.value(x), EPS);

    x = Double.POSITIVE_INFINITY;
    Assert.assertEquals("x=" + x, k, f.value(x), EPS);
}
 

/**
 * tests the firstDerivative function by comparison
 *
 * <p>This will test the functions
 * <tt>f(x) = x^3 - 2x^2 + 6x + 3, g(x) = 3x^2 - 4x + 6</tt>
 * and <tt>h(x) = 6x - 4</tt>
 */
@Test
public void testfirstDerivativeComparison() {
    double[] f_coeff = { 3, 6, -2, 1 };
    double[] g_coeff = { 6, -4, 3 };
    double[] h_coeff = { -4, 6 };

    PolynomialFunction f = new PolynomialFunction(f_coeff);
    PolynomialFunction g = new PolynomialFunction(g_coeff);
    PolynomialFunction h = new PolynomialFunction(h_coeff);

    // compare f' = g
    Assert.assertEquals(f.derivative().value(0), g.value(0), tolerance);
    Assert.assertEquals(f.derivative().value(1), g.value(1), tolerance);
    Assert.assertEquals(f.derivative().value(100), g.value(100), tolerance);
    Assert.assertEquals(f.derivative().value(4.1), g.value(4.1), tolerance);
    Assert.assertEquals(f.derivative().value(-3.25), g.value(-3.25), tolerance);

    // compare g' = h
    Assert.assertEquals(g.derivative().value(FastMath.PI), h.value(FastMath.PI), tolerance);
    Assert.assertEquals(g.derivative().value(FastMath.E),  h.value(FastMath.E),  tolerance);
}
 

/**
 * Returns the natural logarithm of the gamma function Γ(x).
 *
 * The implementation of this method is based on:
 * <ul>
 * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">
 * Gamma Function</a>, equation (28).</li>
 * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
 * Lanczos Approximation</a>, equations (1) through (5).</li>
 * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
 * the computation of the convergent Lanczos complex Gamma approximation
 * </a></li>
 * </ul>
 *
 * @param x the value.
 * @return log(&#915;(x))
 */
public static double logGamma(double x) {
    double ret;

    if (Double.isNaN(x) || (x <= 0.0)) {
        ret = Double.NaN;
    } else {
        double g = 607.0 / 128.0;

        double sum = 0.0;
        for (int i = LANCZOS.length - 1; i > 0; --i) {
            sum = sum + (LANCZOS[i] / (x + i));
        }
        sum = sum + LANCZOS[0];

        double tmp = x + g + .5;
        ret = ((x + .5) * FastMath.log(tmp)) - tmp +
            HALF_LOG_2_PI + FastMath.log(sum / x);
    }

    return ret;
}
 

/** Reset the instance as if built from two points.
 * <p>The line is oriented from p1 to p2</p>
 * @param p1 first point
 * @param p2 second point
 */
public void reset(final Vector2D p1, final Vector2D p2) {
    final double dx = p2.getX() - p1.getX();
    final double dy = p2.getY() - p1.getY();
    final double d = FastMath.hypot(dx, dy);
    if (d == 0.0) {
        angle        = 0.0;
        cos          = 1.0;
        sin          = 0.0;
        originOffset = p1.getY();
    } else {
        angle        = FastMath.PI + FastMath.atan2(-dy, -dx);
        cos          = FastMath.cos(angle);
        sin          = FastMath.sin(angle);
        originOffset = (p2.getX() * p1.getY() - p1.getX() * p2.getY()) / d;
    }
}
 

/**
 * Test of integrator for the quintic function.
 */
public void testQuinticFunction() throws MathException {
    UnivariateRealFunction f = new QuinticFunction();
    UnivariateRealIntegrator integrator = new SimpsonIntegrator();
    double min, max, expected, result, tolerance;

    min = 0; max = 1; expected = -1.0/48;
    tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
    result = integrator.integrate(f, min, max);
    assertEquals(expected, result, tolerance);

    min = 0; max = 0.5; expected = 11.0/768;
    tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
    result = integrator.integrate(f, min, max);
    assertEquals(expected, result, tolerance);

    min = -1; max = 4; expected = 2048/3.0 - 78 + 1.0/48;
    tolerance = FastMath.abs(expected * integrator.getRelativeAccuracy());
    result = integrator.integrate(f, min, max);
    assertEquals(expected, result, tolerance);
}
 

/** {@inheritDoc} */
@Override
public double walkInRowOrder(final RealMatrixPreservingVisitor visitor,
                             final int startRow, final int endRow,
                             final int startColumn, final int endColumn)
    throws MatrixIndexException, MatrixVisitorException {
    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 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 double[] 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();
}
 

/**
 * {@inheritDoc}
 */
@Override
public double getResult() {
    if (sumOfLogs.getN() > 0) {
        return FastMath.exp(sumOfLogs.getResult() / sumOfLogs.getN());
    } else {
        return Double.NaN;
    }
}
 

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

/**
 * Computes the n-th roots of this complex number.
 * The nth roots are defined by the formula:
 * <pre>
 *  <code>
 *   z<sub>k</sub> = abs<sup>1/n</sup> (cos(phi + 2&pi;k/n) + i (sin(phi + 2&pi;k/n))
 *  </code>
 * </pre>
 * for <i>{@code k=0, 1, ..., n-1}</i>, where {@code abs} and {@code phi}
 * are respectively the {@link #abs() modulus} and
 * {@link #getArgument() argument} of this complex number.
 * <br/>
 * If one or both parts of this complex number is NaN, a list with just
 * one element, {@link #NaN} is returned.
 * if neither part is NaN, but at least one part is infinite, the result
 * is a one-element list containing {@link #INF}.
 *
 * @param n Degree of root.
 * @return a List<Complex> of all {@code n}-th roots of {@code this}.
 * @throws NotPositiveException if {@code n <= 0}.
 * @since 2.0
 */
public List<Complex> nthRoot(int n) {

    if (n <= 0) {
        throw new NotPositiveException(LocalizedFormats.CANNOT_COMPUTE_NTH_ROOT_FOR_NEGATIVE_N,
                                       n);
    }

    final List<Complex> result = new ArrayList<Complex>();

    if (isNaN) {
        result.add(NaN);
        return result;
    }
    if (isInfinite()) {
        result.add(INF);
        return result;
    }

    // nth root of abs -- faster / more accurate to use a solver here?
    final double nthRootOfAbs = FastMath.pow(abs(), 1.0 / n);

    // Compute nth roots of complex number with k = 0, 1, ... n-1
    final double nthPhi = getArgument() / n;
    final double slice = 2 * FastMath.PI / n;
    double innerPart = nthPhi;
    for (int k = 0; k < n ; k++) {
        // inner part
        final double realPart = nthRootOfAbs *  FastMath.cos(innerPart);
        final double imaginaryPart = nthRootOfAbs *  FastMath.sin(innerPart);
        result.add(createComplex(realPart, imaginaryPart));
        innerPart += slice;
    }

    return result;
}
 

/**
 * Computes the n-th roots of this complex number.
 * The nth roots are defined by the formula:
 * <pre>
 *  <code>
 *   z<sub>k</sub> = abs<sup>1/n</sup> (cos(phi + 2&pi;k/n) + i (sin(phi + 2&pi;k/n))
 *  </code>
 * </pre>
 * for <i>{@code k=0, 1, ..., n-1}</i>, where {@code abs} and {@code phi}
 * are respectively the {@link #abs() modulus} and
 * {@link #getArgument() argument} of this complex number.
 * <br/>
 * If one or both parts of this complex number is NaN, a list with just
 * one element, {@link #NaN} is returned.
 * if neither part is NaN, but at least one part is infinite, the result
 * is a one-element list containing {@link #INF}.
 *
 * @param n Degree of root.
 * @return a List<Complex> of all {@code n}-th roots of {@code this}.
 * @throws NotPositiveException if {@code n <= 0}.
 * @since 2.0
 */
public List<Complex> nthRoot(int n) {

    if (n <= 0) {
        throw new NotPositiveException(LocalizedFormats.CANNOT_COMPUTE_NTH_ROOT_FOR_NEGATIVE_N,
                                       n);
    }

    final List<Complex> result = new ArrayList<Complex>();

    if (isNaN) {
        result.add(NaN);
        return result;
    }
    if (isInfinite()) {
        result.add(INF);
        return result;
    }

    // nth root of abs -- faster / more accurate to use a solver here?
    final double nthRootOfAbs = FastMath.pow(abs(), 1.0 / n);

    // Compute nth roots of complex number with k = 0, 1, ... n-1
    final double nthPhi = getArgument() / n;
    final double slice = 2 * FastMath.PI / n;
    double innerPart = nthPhi;
    for (int k = 0; k < n ; k++) {
        // inner part
        final double realPart = nthRootOfAbs *  FastMath.cos(innerPart);
        final double imaginaryPart = nthRootOfAbs *  FastMath.sin(innerPart);
        result.add(createComplex(realPart, imaginaryPart));
        innerPart += slice;
    }

    return result;
}
 

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

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

/**
 * {@inheritDoc}
 */
@Override
public double density(double x) {
    if (x < 0) {
        return 0;
    }
    return FastMath.exp(-x / mean) / mean;
}
 

/** {@inheritDoc} */
@Override
public double getDistance(double[] v) {
    checkVectorDimensions(v.length);
    double res = 0;
    for (int i = 0; i < v.length; i++) {
        double delta = entries.get(i) - v[i];
        res += delta * delta;
    }
    return FastMath.sqrt(res);
}
 

public void testSqrtPolar() {
    double r = 1;
    for (int i = 0; i < 5; i++) {
        r += i;
        double theta = 0;
        for (int j =0; j < 11; j++) {
            theta += pi /12;
            Complex z = ComplexUtils.polar2Complex(r, theta);
            Complex sqrtz = ComplexUtils.polar2Complex(FastMath.sqrt(r), theta / 2);
            TestUtils.assertEquals(sqrtz, z.sqrt(), 10e-12);
        }
    }
}
 

/** {@inheritDoc} */
public double solve(final UnivariateRealFunction f, double min, double max)
    throws MaxIterationsExceededException, FunctionEvaluationException {

    clearResult();
    verifyInterval(min,max);
    double m;
    double fm;
    double fmin;

    int i = 0;
    while (i < maximalIterationCount) {
        m = UnivariateRealSolverUtils.midpoint(min, max);
       fmin = f.value(min);
       fm = f.value(m);

        if (fm * fmin > 0.0) {
            // max and m bracket the root.
            min = m;
        } else {
            // min and m bracket the root.
            max = m;
        }

        if (FastMath.abs(max - min) <= absoluteAccuracy) {
            m = UnivariateRealSolverUtils.midpoint(min, max);
            setResult(m, i);
            return m;
        }
        ++i;
    }

    throw new MaxIterationsExceededException(maximalIterationCount);
}
 

@Test
public void testIncreasingRobustnessItersIncreasesSmoothnessWithOutliers() {
    int numPoints = 100;
    double[] xval = new double[numPoints];
    double[] yval = new double[numPoints];
    double xnoise = 0.1;
    double ynoise = 0.1;

    generateSineData(xval, yval, xnoise, ynoise);

    // Introduce a couple of outliers
    yval[numPoints/3] *= 100;
    yval[2 * numPoints/3] *= -100;

    // Check that variance decreases as the number of robustness
    // iterations increases

    double[] variances = new double[4];
    for (int i = 0; i < 4; i++) {
        LoessInterpolator li = new LoessInterpolator(0.3, i, 1e-12);

        double[] res = li.smooth(xval, yval);

        for (int j = 1; j < res.length; ++j) {
            variances[i] += FastMath.abs(res[j] - res[j-1]);
        }
    }

    for(int i = 1; i < variances.length; ++i) {
        Assert.assertTrue(variances[i] < variances[i-1]);
    }
}
 

/** {@inheritDoc} */
@Override
public double getDistance(double[] v) {
    checkVectorDimensions(v.length);
    double res = 0;
    for (int i = 0; i < v.length; i++) {
        double delta = entries.get(i) - v[i];
        res += delta * delta;
    }
    return FastMath.sqrt(res);
}
 

public double[] exactY(double t) {
    double cos = FastMath.cos(omega * t);
    double sin = FastMath.sin(omega * t);
    double dx0 = y0[0] - cx;
    double dy0 = y0[1] - cy;
    return new double[] {
        cx + cos * dx0 - sin * dy0,
        cy + sin * dx0 + cos * dy0
    };
}
 

/** {@inheritDoc} */
@Override
public double getDistance(double[] v) {
    checkVectorDimensions(v.length);
    double res = 0;
    for (int i = 0; i < v.length; i++) {
        double delta = entries.get(i) - v[i];
        res += delta * delta;
    }
    return FastMath.sqrt(res);
}
 

@Test
public void testDistance() {
    EuclideanIntegerPoint e1 = new EuclideanIntegerPoint(new int[] { -3, -2, -1, 0, 1 });
    EuclideanIntegerPoint e2 = new EuclideanIntegerPoint(new int[] {  1,  0, -1, 1, 1 });
    assertEquals(FastMath.sqrt(21.0), e1.distanceFrom(e2), 1.0e-15);
    assertEquals(0.0, e1.distanceFrom(e1), 1.0e-15);
    assertEquals(0.0, e2.distanceFrom(e2), 1.0e-15);
}
 

/**
 * Computes the n-th roots of this complex number.
 * The nth roots are defined by the formula:
 * <pre>
 *  <code>
 *   z<sub>k</sub> = abs<sup>1/n</sup> (cos(phi + 2&pi;k/n) + i (sin(phi + 2&pi;k/n))
 *  </code>
 * </pre>
 * for <i>{@code k=0, 1, ..., n-1}</i>, where {@code abs} and {@code phi}
 * are respectively the {@link #abs() modulus} and
 * {@link #getArgument() argument} of this complex number.
 * <br/>
 * If one or both parts of this complex number is NaN, a list with just
 * one element, {@link #NaN} is returned.
 * if neither part is NaN, but at least one part is infinite, the result
 * is a one-element list containing {@link #INF}.
 *
 * @param n Degree of root.
 * @return a List<Complex> of all {@code n}-th roots of {@code this}.
 * @throws NotPositiveException if {@code n <= 0}.
 * @since 2.0
 */
public List<Complex> nthRoot(int n) {

    if (n <= 0) {
        throw new NotPositiveException(LocalizedFormats.CANNOT_COMPUTE_NTH_ROOT_FOR_NEGATIVE_N,
                                       n);
    }

    final List<Complex> result = new ArrayList<Complex>();

    if (isNaN) {
        result.add(NaN);
        return result;
    }
    if (isInfinite()) {
        result.add(INF);
        return result;
    }

    // nth root of abs -- faster / more accurate to use a solver here?
    final double nthRootOfAbs = FastMath.pow(abs(), 1.0 / n);

    // Compute nth roots of complex number with k = 0, 1, ... n-1
    final double nthPhi = getArgument() / n;
    final double slice = 2 * FastMath.PI / n;
    double innerPart = nthPhi;
    for (int k = 0; k < n ; k++) {
        // inner part
        final double realPart = nthRootOfAbs *  FastMath.cos(innerPart);
        final double imaginaryPart = nthRootOfAbs *  FastMath.sin(innerPart);
        result.add(createComplex(realPart, imaginaryPart));
        innerPart += slice;
    }

    return result;
}