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

private void doTestTransformFunction(final int n, final double tol,
    final TransformType type) {
    final RealTransformer transformer = createRealTransformer();
    final UnivariateFunction f = getValidFunction();
    final double a = getValidLowerBound();
    final double b = getValidUpperBound();
    final double[] x = createRealData(n);
    for (int i = 0; i < n; i++) {
        final double t = a + i * (b - a) / n;
        x[i] = f.value(t);
    }
    final double[] expected = transform(x, type);
    final double[] actual = transformer.transform(f, a, b, n, type);
    for (int i = 0; i < n; i++) {
        final String msg = String.format("%d, %d", n, i);
        final double delta = tol * FastMath.abs(expected[i]);
        Assert.assertEquals(msg, expected[i], actual[i], delta);
    }
}
 

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 / 450.0));
    Assert.assertTrue(maxStep > (1.0 / 4.2));
  }
}
 

@Test
public void testNoError() {
    Random randomizer = new Random(64925784252l);
    for (int degree = 1; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
        for (int i = 0; i <= degree; ++i) {
            fitter.addObservedPoint(1.0, i, p.value(i));
        }

        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)));
            Assert.assertEquals(0.0, error, 1.0e-6);
        }
    }
}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve()
    throws TooManyEvaluationsException {
    final double startValue = getStartValue();
    final double absoluteAccuracy = getAbsoluteAccuracy();

    double x0 = startValue;
    double x1;
    while (true) {
        final DerivativeStructure y0 = computeObjectiveValueAndDerivative(x0);
        x1 = x0 - (y0.getValue() / y0.getPartialDerivative(1));
        if (FastMath.abs(x1 - x0) <= absoluteAccuracy) {
            return x1;
        }

        x0 = x1;
    }
}
 

@Test
public void testNoError() {
    Random randomizer = new Random(64925784252l);
    for (int degree = 1; degree < 10; ++degree) {
        PolynomialFunction p = buildRandomPolynomial(degree, randomizer);

        PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer());
        for (int i = 0; i <= degree; ++i) {
            fitter.addObservedPoint(1.0, i, p.value(i));
        }

        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)));
            Assert.assertEquals(0.0, error, 1.0e-6);
        }
    }
}
 

/**
 * A part of the deviance portion of the saddle point approximation.
 * <p>
 * References:
 * <ol>
 * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
 * Probabilities.". <a target="_blank"
 * href="http://www.herine.net/stat/papers/dbinom.pdf">
 * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
 * </ol>
 * </p>
 *
 * @param x the x value.
 * @param mu the average.
 * @return a part of the deviance.
 */
static double getDeviancePart(double x, double mu) {
    double ret;
    if (FastMath.abs(x - mu) < 0.1 * (x + mu)) {
        double d = x - mu;
        double v = d / (x + mu);
        double s1 = v * d;
        double s = Double.NaN;
        double ej = 2.0 * x * v;
        v *= v;
        int j = 1;
        while (s1 != s) {
            s = s1;
            ej *= v;
            s1 = s + ej / ((j * 2) + 1);
            ++j;
        }
        ret = s1;
    } else {
        ret = x * FastMath.log(x / mu) + mu - x;
    }
    return ret;
}
 

/**
 * Returns true iff there is a column that is a scalar multiple of column
 * in searchMatrix (modulo tolerance)
 */
private boolean isIncludedColumn(double[] column, RealMatrix searchMatrix,
        double tolerance) {
    boolean found = false;
    int i = 0;
    while (!found && i < searchMatrix.getColumnDimension()) {
        double multiplier = 1.0;
        boolean matching = true;
        int j = 0;
        while (matching && j < searchMatrix.getRowDimension()) {
            double colEntry = searchMatrix.getEntry(j, i);
            // Use the first entry where both are non-zero as scalar
            if (FastMath.abs(multiplier - 1.0) <= FastMath.ulp(1.0) && FastMath.abs(colEntry) > 1E-14
                    && FastMath.abs(column[j]) > 1e-14) {
                multiplier = colEntry / column[j];
            }
            if (FastMath.abs(column[j] * multiplier - colEntry) > tolerance) {
                matching = false;
            }
            j++;
        }
        found = matching;
        i++;
    }
    return found;
}
 

/** Get the intersection of a line with the instance.
 * @param line line intersecting the instance
 * @return intersection point between between the line and the
 * instance (null if the line is parallel to the instance)
 */
public Vector3D intersection(final Line line) {
    final Vector3D direction = line.getDirection();
    final double   dot       = w.dotProduct(direction);
    if (FastMath.abs(dot) < 1.0e-10) {
        return null;
    }
    final Vector3D point = line.toSpace(Vector1D.ZERO);
    final double   k     = -(originOffset + w.dotProduct(point)) / dot;
    return new Vector3D(1.0, point, k, direction);
}
 

/**
 * Test of interpolator for the sine function.
 * <p>
 * |sin^(n)(zeta)| <= 1.0, zeta in [0, 2*PI]
 */
@Test
public void testSinFunction() {
    UnivariateFunction f = new SinFunction();
    UnivariateInterpolator interpolator = new DividedDifferenceInterpolator();
    double x[], y[], z, expected, result, tolerance;

    // 6 interpolating points on interval [0, 2*PI]
    int n = 6;
    double min = 0.0, max = 2 * FastMath.PI;
    x = new double[n];
    y = new double[n];
    for (int i = 0; i < n; i++) {
        x[i] = min + i * (max - min) / n;
        y[i] = f.value(x[i]);
    }
    double derivativebound = 1.0;
    UnivariateFunction p = interpolator.interpolate(x, y);

    z = FastMath.PI / 4; expected = f.value(z); result = p.value(z);
    tolerance = FastMath.abs(derivativebound * partialerror(x, z));
    Assert.assertEquals(expected, result, tolerance);

    z = FastMath.PI * 1.5; expected = f.value(z); result = p.value(z);
    tolerance = FastMath.abs(derivativebound * partialerror(x, z));
    Assert.assertEquals(expected, result, tolerance);
}
 

/** Get the intersection point of the instance and another line.
 * @param other other line
 * @return intersection point of the instance and the other line
 * or null if there are no intersection points
 */
public Vector2D intersection(final Line other) {
    final double d = sin * other.cos - other.sin * cos;
    if (FastMath.abs(d) < 1.0e-10) {
        return null;
    }
    return new Vector2D((cos * other.originOffset - other.cos * originOffset) / d,
                        (sin * other.originOffset - other.sin * originOffset) / d);
}
 

/**
 * {@inheritDoc}
 *
 * 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.
 */
public double cumulativeProbability(double x)  {
    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 * SQRT2)));
}
 

/**
 * {@inheritDoc}
 *
 * 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.
 */
public double cumulativeProbability(double x)  {
    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 * SQRT2)));
}
 

private float getCoefficient(List<Float2FloatKeyValue> scores) {
    float coefficient = 0F;
    for (Float2FloatKeyValue term : scores) {
        float distance = term.getKey() - term.getValue();
        coefficient += FastMath.abs(distance);
    }
    return coefficient;
}
 

/** {@inheritDoc} */
public double distance1(Vector<Euclidean2D> p) {
    Vector2D p3 = (Vector2D) p;
    final double dx = FastMath.abs(p3.x - x);
    final double dy = FastMath.abs(p3.y - y);
    return dx + dy;
}
 

/**
 * Perform a double QR step involving rows l:idx and columns m:n
 *
 * @param il the index of the small sub-diagonal element
 * @param im the start index for the QR step
 * @param iu the current eigenvalue index
 * @param shift shift information holder
 * @param hVec the initial houseHolder vector
 */
private void performDoubleQRStep(final int il, final int im, final int iu,
                                 final ShiftInfo shift, final double[] hVec) {

    final int n = matrixT.length;
    double p = hVec[0];
    double q = hVec[1];
    double r = hVec[2];

    for (int k = im; k <= iu - 1; k++) {
        boolean notlast = k != (iu - 1);
        if (k != im) {
            p = matrixT[k][k - 1];
            q = matrixT[k + 1][k - 1];
            r = notlast ? matrixT[k + 2][k - 1] : 0.0;
            shift.x = FastMath.abs(p) + FastMath.abs(q) + FastMath.abs(r);
            if (Precision.equals(shift.x, 0.0, epsilon)) {
                continue;
            }
            p /= shift.x;
            q /= shift.x;
            r /= shift.x;
        }
        double s = FastMath.sqrt(p * p + q * q + r * r);
        if (p < 0.0) {
            s = -s;
        }
        if (s != 0.0) {
            if (k != im) {
                matrixT[k][k - 1] = -s * shift.x;
            } else if (il != im) {
                matrixT[k][k - 1] = -matrixT[k][k - 1];
            }
            p += s;
            shift.x = p / s;
            shift.y = q / s;
            double z = r / s;
            q /= p;
            r /= p;

            // Row modification
            for (int j = k; j < n; j++) {
                p = matrixT[k][j] + q * matrixT[k + 1][j];
                if (notlast) {
                    p += r * matrixT[k + 2][j];
                    matrixT[k + 2][j] -= p * z;
                }
                matrixT[k][j] -= p * shift.x;
                matrixT[k + 1][j] -= p * shift.y;
            }

            // Column modification
            for (int i = 0; i <= FastMath.min(iu, k + 3); i++) {
                p = shift.x * matrixT[i][k] + shift.y * matrixT[i][k + 1];
                if (notlast) {
                    p += z * matrixT[i][k + 2];
                    matrixT[i][k + 2] -= p * r;
                }
                matrixT[i][k] -= p;
                matrixT[i][k + 1] -= p * q;
            }

            // Accumulate transformations
            final int high = matrixT.length - 1;
            for (int i = 0; i <= high; i++) {
                p = shift.x * matrixP[i][k] + shift.y * matrixP[i][k + 1];
                if (notlast) {
                    p += z * matrixP[i][k + 2];
                    matrixP[i][k + 2] -= p * r;
                }
                matrixP[i][k] -= p;
                matrixP[i][k + 1] -= p * q;
            }
        }  // (s != 0)
    }  // k loop

    // clean up pollution due to round-off errors
    for (int i = im + 2; i <= iu; i++) {
        matrixT[i][i-2] = 0.0;
        if (i > im + 2) {
            matrixT[i][i-3] = 0.0;
        }
    }
}
 

/** Append another model at the end of the instance.
 * @param model model to add at the end of the instance
 * @exception MathIllegalArgumentException if the model to append is not
 * compatible with the instance (dimension of the state vector,
 * propagation direction, hole between the dates)
 */
public void append(final ContinuousOutputModel model)
  throws MathIllegalArgumentException {

  if (model.steps.size() == 0) {
    return;
  }

  if (steps.size() == 0) {
    initialTime = model.initialTime;
    forward     = model.forward;
  } else {

    if (getInterpolatedState().length != model.getInterpolatedState().length) {
        throw new DimensionMismatchException(model.getInterpolatedState().length,
                                             getInterpolatedState().length);
    }

    if (forward ^ model.forward) {
        throw new MathIllegalArgumentException(LocalizedFormats.PROPAGATION_DIRECTION_MISMATCH);
    }

    final StepInterpolator lastInterpolator = steps.get(index);
    final double current  = lastInterpolator.getCurrentTime();
    final double previous = lastInterpolator.getPreviousTime();
    final double step = current - previous;
    final double gap = model.getInitialTime() - current;
    if (FastMath.abs(gap) > 1.0e-3 * FastMath.abs(step)) {
      throw new MathIllegalArgumentException(LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES,
                                             FastMath.abs(gap));
    }

  }

  for (StepInterpolator interpolator : model.steps) {
    steps.add(interpolator.copy());
  }

  index = steps.size() - 1;
  finalTime = (steps.get(index)).getCurrentTime();

}
 

/**
 * {@inheritDoc}
 */
@Override
protected double doSolve()
    throws TooManyEvaluationsException,
           NumberIsTooLargeException,
           NoBracketingException {
    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;
    }
}
 

/**
 * 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.math3.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;
    }
    final double ret = Gamma.regularizedGammaP(0.5, x * x, 1.0e-15, 10000);
    return x < 0 ? -ret : ret;
}
 

/**
 * Checks whether the instance is a pure quaternion within a given
 * tolerance.
 *
 * @param eps Tolerance (absolute error).
 * @return {@code true} if the scalar part of the quaternion is zero.
 */
public boolean isPureQuaternion(double eps) {
    return FastMath.abs(getQ0()) <= eps;
}
 

/**
 * Returns {@code true} if there is no double value strictly between the
 * arguments or the difference between them is within the range of allowed
 * error (inclusive).
 *
 * @param x First value.
 * @param y Second value.
 * @param eps Amount of allowed absolute error.
 * @return {@code true} if the values are two adjacent floating point
 * numbers or they are within range of each other.
 */
public static boolean equals(double x, double y, double eps) {
    return equals(x, y, 1) || FastMath.abs(y - x) <= eps;
}