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

@Test
public void testNormalVariance() {
    final double twoOverSqrtPi = 2 / FastMath.sqrt(Math.PI);

    final double mu = 12345.6789;
    final double sigma = 987.654321;
    final double sigma2 = sigma * sigma;
    final int numPoints = 5;

    // Change of variable:
    //   y = (x - mu) / (sqrt(2) *  sigma)
    // such that the integrand
    //   (x - mu)^2 * N(x, mu, sigma)
    // is transformed to
    //   f(y) * exp(-y^2)
    final UnivariateFunction f = new UnivariateFunction() {
            public double value(double y) {
                return twoOverSqrtPi * sigma2 * y * y;
            }
        };

    final GaussIntegrator integrator = factory.hermite(numPoints);
    final double result = integrator.integrate(f);
    final double expected = sigma2;
    Assert.assertEquals(expected, result, 10 * Math.ulp(expected));
}
 

@Test
public void testStair() {
    Vector2D[][] vertices = new Vector2D[][] {
        new Vector2D[] {
            new Vector2D( 0.0, 0.0),
            new Vector2D( 0.0, 2.0),
            new Vector2D(-0.1, 2.0),
            new Vector2D(-0.1, 1.0),
            new Vector2D(-0.3, 1.0),
            new Vector2D(-0.3, 1.5),
            new Vector2D(-1.3, 1.5),
            new Vector2D(-1.3, 2.0),
            new Vector2D(-1.8, 2.0),
            new Vector2D(-1.8 - 1.0 / FastMath.sqrt(2.0),
                        2.0 - 1.0 / FastMath.sqrt(2.0))
        }
    };

    PolygonsSet set = buildSet(vertices);
    checkVertices(set.getVertices(), vertices);

    Assert.assertEquals(1.1 + 0.95 * FastMath.sqrt(2.0), set.getSize(), 1.0e-10);

}
 

/**
 * 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);

}
 

/**
 * Update the residuals array and cost function value.
 * @throws DimensionMismatchException if the dimension does not match the
 * problem dimension.
 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
 * if the maximal number of evaluations is exceeded.
 */
protected void updateResidualsAndCost() {
    objective = computeObjectiveValue(point);
    if (objective.length != rows) {
        throw new DimensionMismatchException(objective.length, rows);
    }

    final double[] targetValues = getTargetRef();
    final double[] residualsWeights = getWeightRef();

    cost = 0;
    for (int i = 0; i < rows; i++) {
        final double residual = targetValues[i] - objective[i];
        weightedResiduals[i]= residual*FastMath.sqrt(residualsWeights[i]);
        cost += residualsWeights[i] * residual * residual;
    }
    cost = FastMath.sqrt(cost);
}
 

/** {@inheritDoc} */
public double nextGaussian() {

    final double random;
    if (Double.isNaN(nextGaussian)) {
        // generate a new pair of gaussian numbers
        final double x = nextDouble();
        final double y = nextDouble();
        final double alpha = 2 * FastMath.PI * x;
        final double r      = FastMath.sqrt(-2 * FastMath.log(y));
        random       = r * FastMath.cos(alpha);
        nextGaussian = r * FastMath.sin(alpha);
    } else {
        // use the second element of the pair already generated
        random = nextGaussian;
        nextGaussian = Double.NaN;
    }

    return random;

}
 

/** Reset the instance as if built from two points.
 * @param p1 first point belonging to the line (this can be any point)
 * @param p2 second point belonging to the line (this can be any point, different from p1)
 * @exception MathIllegalArgumentException if the points are equal
 */
public void reset(final Vector3D p1, final Vector3D p2) throws MathIllegalArgumentException {
    final Vector3D delta = p2.subtract(p1);
    final double norm2 = delta.getNormSq();
    if (norm2 == 0.0) {
        throw new MathIllegalArgumentException(LocalizedFormats.ZERO_NORM);
    }
    this.direction = new Vector3D(1.0 / FastMath.sqrt(norm2), delta);
    zero = new Vector3D(1.0, p1, -p1.dotProduct(delta) / norm2, delta);
}
 

@Override
public void doComputeDerivatives(double t, double[] y, double[] yDot) {

  // current radius
  double r2 = y[0] * y[0] + y[1] * y[1];
  double invR3 = 1 / (r2 * FastMath.sqrt(r2));

  // compute the derivatives
  yDot[0] = y[2];
  yDot[1] = y[3];
  yDot[2] = -invR3  * y[0];
  yDot[3] = -invR3  * y[1];

}
 

@Test
public void testErf1960() {
    double x = 1.960 / FastMath.sqrt(2.0);
    double actual = Erf.erf(x);
    double expected = 0.95;
    Assert.assertEquals(expected, actual, 1.0e-5);
    Assert.assertEquals(1 - actual, Erf.erfc(x), 1.0e-15);

    actual = Erf.erf(-x);
    expected = -expected;
    Assert.assertEquals(expected, actual, 1.0e-5);
    Assert.assertEquals(1 - actual, Erf.erfc(-x), 1.0e-15);
}
 

/** {@inheritDoc} */
public Line apply(final Hyperplane<Euclidean2D> hyperplane) {
    final Line   line    = (Line) hyperplane;
    final double rOffset = c1X * line.cos + c1Y * line.sin + c11 * line.originOffset;
    final double rCos    = cXX * line.cos + cXY * line.sin;
    final double rSin    = cYX * line.cos + cYY * line.sin;
    final double inv     = 1.0 / FastMath.sqrt(rSin * rSin + rCos * rCos);
    return new Line(FastMath.PI + FastMath.atan2(-rSin, -rCos),
                    inv * rCos, inv * rSin,
                    inv * rOffset);
}
 

/** {@inheritDoc} */
public double getNumericalMean() {
    return Gamma.gamma(mu + 0.5) / Gamma.gamma(mu) * FastMath.sqrt(omega / mu);
}
 

public LinearRank1ZeroColsAndRowsFunction(int m, int n, double x0) {
    super(m, buildArray(n, x0),
          FastMath.sqrt((m * (m + 3) - 6) / (2.0 * (2 * m - 3))),
          null);
}
 

/**
 * Non-linear test case: fitting of decay curve (from Chapter 8 of
 * Bevington's textbook, "Data reduction and analysis for the physical sciences").
 * XXX The expected ("reference") values may not be accurate and the tolerance too
 * relaxed for this test to be currently really useful (the issue is under
 * investigation).
 */
@Test
public void testBevington() {
    final double[][] dataPoints = {
        // column 1 = times
        { 15, 30, 45, 60, 75, 90, 105, 120, 135, 150,
          165, 180, 195, 210, 225, 240, 255, 270, 285, 300,
          315, 330, 345, 360, 375, 390, 405, 420, 435, 450,
          465, 480, 495, 510, 525, 540, 555, 570, 585, 600,
          615, 630, 645, 660, 675, 690, 705, 720, 735, 750,
          765, 780, 795, 810, 825, 840, 855, 870, 885, },
        // column 2 = measured counts
        { 775, 479, 380, 302, 185, 157, 137, 119, 110, 89,
          74, 61, 66, 68, 48, 54, 51, 46, 55, 29,
          28, 37, 49, 26, 35, 29, 31, 24, 25, 35,
          24, 30, 26, 28, 21, 18, 20, 27, 17, 17,
          14, 17, 24, 11, 22, 17, 12, 10, 13, 16,
          9, 9, 14, 21, 17, 13, 12, 18, 10, },
    };
    final double[] start = {10, 900, 80, 27, 225};

    final BevingtonProblem problem = new BevingtonProblem();

    final int len = dataPoints[0].length;
    final double[] weights = new double[len];
    for (int i = 0; i < len; i++) {
        problem.addPoint(dataPoints[0][i],
                         dataPoints[1][i]);

        weights[i] = 1 / dataPoints[1][i];
    }

    final Optimum optimum = optimizer.optimize(
            builder(problem)
                    .target(dataPoints[1])
                    .weight(new DiagonalMatrix(weights))
                    .start(start)
                    .maxIterations(20)
                    .build()
    );

    final RealVector solution = optimum.getPoint();
    final double[] expectedSolution = { 10.4, 958.3, 131.4, 33.9, 205.0 };

    final RealMatrix covarMatrix = optimum.getCovariances(1e-14);
    final double[][] expectedCovarMatrix = {
        { 3.38, -3.69, 27.98, -2.34, -49.24 },
        { -3.69, 2492.26, 81.89, -69.21, -8.9 },
        { 27.98, 81.89, 468.99, -44.22, -615.44 },
        { -2.34, -69.21, -44.22, 6.39, 53.80 },
        { -49.24, -8.9, -615.44, 53.8, 929.45 }
    };

    final int numParams = expectedSolution.length;

    // Check that the computed solution is within the reference error range.
    for (int i = 0; i < numParams; i++) {
        final double error = FastMath.sqrt(expectedCovarMatrix[i][i]);
        Assert.assertEquals("Parameter " + i, expectedSolution[i], solution.getEntry(i), error);
    }

    // Check that each entry of the computed covariance matrix is within 10%
    // of the reference matrix entry.
    for (int i = 0; i < numParams; i++) {
        for (int j = 0; j < numParams; j++) {
            Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]",
                                expectedCovarMatrix[i][j],
                                covarMatrix.getEntry(i, j),
                                FastMath.abs(0.1 * expectedCovarMatrix[i][j]));
        }
    }
}
 

/** {@inheritDoc} */
public SparseGradient asinh() {
    return new SparseGradient(FastMath.asinh(value), 1.0 / FastMath.sqrt(value * value + 1.0), derivatives);
}
 

/** Compute arc sine of a derivative structure.
 * @param operand array holding the operand
 * @param operandOffset offset of the operand in its array
 * @param result array where result must be stored (for
 * arc sine the result array <em>cannot</em> be the input
 * array)
 * @param resultOffset offset of the result in its array
 */
public void asin(final double[] operand, final int operandOffset,
                final double[] result, final int resultOffset) {

    // create the function value and derivatives
    double[] function = new double[1 + order];
    final double x = operand[operandOffset];
    function[0] = FastMath.asin(x);
    if (order > 0) {
        // the nth order derivative of asin has the form:
        // dn(asin(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
        // where P_n(x) is a degree n-1 polynomial with same parity as n-1
        // P_1(x) = 1, P_2(x) = x, P_3(x) = 2x^2 + 1 ...
        // the general recurrence relation for P_n is:
        // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
        // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
        final double[] p = new double[order];
        p[0] = 1;
        final double x2    = x * x;
        final double f     = 1.0 / (1 - x2);
        double coeff = FastMath.sqrt(f);
        function[1] = coeff * p[0];
        for (int n = 2; n <= order; ++n) {

            // update and evaluate polynomial P_n(x)
            double v = 0;
            p[n - 1] = (n - 1) * p[n - 2];
            for (int k = n - 1; k >= 0; k -= 2) {
                v = v * x2 + p[k];
                if (k > 2) {
                    p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
                } else if (k == 2) {
                    p[0] = p[1];
                }
            }
            if ((n & 0x1) == 0) {
                v *= x;
            }

            coeff *= f;
            function[n] = coeff * v;

        }
    }

    // apply function composition
    compose(operand, operandOffset, function, result, resultOffset);

}
 

/**
 * <p>
 * Returns an estimate of the standard deviation of each parameter. The
 * returned values are the so-called (asymptotic) standard errors on the
 * parameters, defined as {@code sd(a[i]) = sqrt(S / (n - m) * C[i][i])},
 * where {@code a[i]} is the optimized value of the {@code i}-th parameter,
 * {@code S} is the minimized value of the sum of squares objective function
 * (as returned by {@link #getChiSquare()}), {@code n} is the number of
 * observations, {@code m} is the number of parameters and {@code C} is the
 * covariance matrix.
 * </p>
 * <p>
 * See also
 * <a href="http://en.wikipedia.org/wiki/Least_squares">Wikipedia</a>,
 * or
 * <a href="http://mathworld.wolfram.com/LeastSquaresFitting.html">MathWorld</a>,
 * equations (34) and (35) for a particular case.
 * </p>
 *
 * @return an estimate of the standard deviation of the optimized parameters
 * @throws org.apache.commons.math3.linear.SingularMatrixException
 * if the covariance matrix cannot be computed.
 * @throws NumberIsTooSmallException if the number of degrees of freedom is not
 * positive, i.e. the number of measurements is less or equal to the number of
 * parameters.
 * @deprecated as of version 3.1, {@link #getSigma()} should be used
 * instead. It should be emphasized that {@link #guessParametersErrors()} and
 * {@link #getSigma()} are <em>not</em> strictly equivalent.
 */
@Deprecated
public double[] guessParametersErrors() {
    if (rows <= cols) {
        throw new NumberIsTooSmallException(LocalizedFormats.NO_DEGREES_OF_FREEDOM,
                                            rows, cols, false);
    }
    double[] errors = new double[cols];
    final double c = FastMath.sqrt(getChiSquare() / (rows - cols));
    double[][] covar = getCovariances();
    for (int i = 0; i < errors.length; ++i) {
        errors[i] = FastMath.sqrt(covar[i][i]) * c;
    }
    return errors;
}
 

/**
 * Computes t test statistic for 2-sample t-test.
 * <p>
 * Does not assume that subpopulation variances are equal.</p>
 *
 * @param m1 first sample mean
 * @param m2 second sample mean
 * @param v1 first sample variance
 * @param v2 second sample variance
 * @param n1 first sample n
 * @param n2 second sample n
 * @return t test statistic
 */
protected double t(final double m1, final double m2,
                   final double v1, final double v2,
                   final double n1, final double n2)  {
    return (m1 - m2) / FastMath.sqrt((v1 / n1) + (v2 / n2));
}
 

/**
 * Compute the <tt>L<sub>2</sub></tt> (Euclidean) norm, or the sqrt
 * of the sum of squared terms in the vector
 * @param a
 * @return the norm
 */
public static double l2Norm(final double[] a) {
	return FastMath.sqrt(innerProduct(a, a));
}
 

/**
 * Computes t test statistic for 1-sample t-test.
 *
 * @param m sample mean
 * @param mu constant to test against
 * @param v sample variance
 * @param n sample n
 * @return t test statistic
 */
protected double t(final double m, final double mu,
                   final double v, final double n) {
    return (m - mu) / FastMath.sqrt(v / n);
}
 

/**
 * Computes t test statistic for 1-sample t-test.
 *
 * @param m sample mean
 * @param mu constant to test against
 * @param v sample variance
 * @param n sample n
 * @return t test statistic
 */
protected double t(final double m, final double mu,
                   final double v, final double n) {
    return (m - mu) / FastMath.sqrt(v / n);
}
 

/**
 * Returns the <a href="http://www.xycoon.com/standerrorb(1).htm">standard
 * error of the slope estimate</a>,
 * usually denoted s(b1).
 * <p>
 * If there are fewer that <strong>three</strong> data pairs in the model,
 * or if there is no variation in x, this returns <code>Double.NaN</code>.
 * </p>
 *
 * @return standard error associated with slope estimate
 */
public double getSlopeStdErr() {
    return FastMath.sqrt(getMeanSquareError() / sumXX);
}