Examples with DoubleDouble org.apache.sis.internal.util.DoubleDouble used on opensource projects

Example 26 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class MatricesTest method testCreateFromNumbers.

/**
 * Tests {@link Matrices#create(int, int, Number[])}.
 */
public void testCreateFromNumbers() {
    final double SENTINEL_VALUE = Double.MIN_VALUE;
    final int SIZE = Matrix3.SIZE;
    final Matrix3 expected = new Matrix3(1, 2, 3, 0.1, 0.2, 0.3, -1, -2, -3);
    final Number[] elements = { 1, 2, 3, 0.1, 0.2, 0.3, -1, -2, -3 };
    /*
         * Mix of Integer and Double objects but without DoubleDouble objects.
         * The result shall be a matrix using the standard double Java type.
         */
    assertEquals(expected, Matrices.create(SIZE, SIZE, elements));
    /*
         * Now put some DoubleDouble instances in the diagonal. We set the error term to
         * Double.MIN_VALUE in order to differentiate them from automatically calculated
         * error terms. The result shall use double-double arithmetic, and we should be
         * able to find back our error terms.
         */
    for (int i = 0; i < elements.length; i += SIZE + 1) {
        elements[i] = new DoubleDouble(elements[i].doubleValue(), SENTINEL_VALUE);
    }
    final MatrixSIS matrix = Matrices.create(SIZE, SIZE, elements);
    assertInstanceOf("Created with DoubleDouble elements", GeneralMatrix.class, matrix);
    // Because not the same type.
    assertFalse(expected.equals(matrix));
    assertTrue(Matrices.equals(expected, matrix, ComparisonMode.BY_CONTRACT));
    final double[] errors = ((GeneralMatrix) matrix).elements;
    for (int i = 0; i < SIZE * SIZE; i++) {
        // Only elements on the diagonal shall be our sentinel value.
        assertEquals((i % (SIZE + 1)) == 0, errors[i + SIZE * SIZE] == SENTINEL_VALUE);
    }
}

Example 27 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class ContextualParameters method denormalizeGeographicOutputs.

/**
 * Appends a denormalization step after the non-linear kernel, which will convert input ordinates
 * in the two first dimensions from radians to degrees. After this conversion, the denormalization
 * can optionally add the given λ₀ value (in degrees) to the longitude.
 *
 * <p>Invoking this method is equivalent to {@linkplain java.awt.geom.AffineTransform#concatenate concatenating}
 * the denormalization matrix with the following matrix. This will have the effect of applying the conversion
 * described above after the non-linear kernel operation:</p>
 *
 * <center>{@include formulas.html#DenormalizeGeographic}</center>
 *
 * @param  λ0  longitude of the central meridian, in degrees.
 * @return the denormalization affine transform as a matrix.
 *         Callers can change that matrix directly if they want to apply additional denormalization operations.
 * @throws IllegalStateException if this {@code ContextualParameter} has been made unmodifiable.
 */
public synchronized MatrixSIS denormalizeGeographicOutputs(final double λ0) {
    ensureModifiable();
    final DoubleDouble toDegrees = DoubleDouble.createRadiansToDegrees();
    // Must be the same instance, not a copy.
    final MatrixSIS denormalize = (MatrixSIS) this.denormalize;
    denormalize.convertAfter(0, toDegrees, (λ0 != 0) ? λ0 : null);
    denormalize.convertAfter(1, toDegrees, null);
    return denormalize;
}

Example 28 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class TransverseMercator method computeCoefficients.

/**
 * Automatically invoked after deserialization for restoring transient fields.
 */
@Override
final void computeCoefficients() {
    /*
         * Double-double precision is necessary for computing 'n', otherwise we have rounding errors
         * in the 3 last digits. Note that we still have sometime a 1 ULP difference compared to the
         * 'n' value at serialization time.
         */
    final DoubleDouble t = new DoubleDouble(1, 0);
    t.subtract(eccentricitySquared, 0);
    t.sqrt();
    t.ratio_1m_1p();
    computeCoefficients(t.doubleValue());
    if (ALLOW_TRIGONOMETRIC_IDENTITIES) {
        // Same scaling than in the constructor.
        cf4 *= 4;
        ci4 *= 4;
        cf6 *= 16;
        ci6 *= 16;
        cf8 *= 64;
        ci8 *= 64;
    }
}

Example 29 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class Solver method solve.

/**
 * Implementation of {@code solve} and {@code inverse} methods.
 * This method contains the code ported from the JAMA package.
 * Use a "left-looking", dot-product, Crout/Doolittle algorithm.
 *
 * <p>This method does <strong>not</strong> checks the matrix size.
 * It is caller's responsibility to ensure that the following hold:</p>
 *
 * {@preformat java
 *   X.getNumRow() == size;
 *   X.getNumCol() == size;
 *   Y.getNumRow() == size;
 *   Y.getNumCol() == innerSize;
 * }
 *
 * @param  LU         elements of the {@code X} matrix to invert, including error terms.
 * @param  Y          the desired result of {@code X} × <var>U</var>.
 * @param  eltY       elements and error terms of the {@code Y} matrix, or {@code null} if not available.
 * @param  size       the value of {@code X.getNumRow()}, {@code X.getNumCol()} and {@code Y.getNumRow()}.
 * @param  innerSize  the value of {@code Y.getNumCol()}.
 * @throws NoninvertibleMatrixException if the {@code X} matrix is not square or singular.
 */
private static MatrixSIS solve(final double[] LU, final Matrix Y, final double[] eltY, final int size, final int innerSize) throws NoninvertibleMatrixException {
    final int errorLU = size * size;
    assert errorLU == GeneralMatrix.indexOfErrors(size, size, LU);
    final int[] pivot = new int[size];
    for (int j = 0; j < size; j++) {
        pivot[j] = j;
    }
    // [0 … size-1] : column values; [size … 2*size-1] : error terms.
    final double[] column = new double[size * 2];
    // Temporary variable for sum ("accumulator") and subtraction.
    final DoubleDouble acc = new DoubleDouble();
    // Temporary variable for products and ratios.
    final DoubleDouble rat = new DoubleDouble();
    for (int i = 0; i < size; i++) {
        /*
             * Make a copy of the i-th column.
             */
        for (int j = 0; j < size; j++) {
            final int k = j * size + i;
            // Value
            column[j] = LU[k];
            // Error
            column[j + size] = LU[k + errorLU];
        }
        /*
             * Apply previous transformations. This part is equivalent to the following code,
             * but using double-double arithmetic instead than the primitive 'double' type:
             *
             *     double sum = 0;
             *     for (int k=0; k<kmax; k++) {
             *         sum += LU[rowOffset + k] * column[k];
             *     }
             *     LU[rowOffset + i] = (column[j] -= sum);
             */
        for (int j = 0; j < size; j++) {
            final int rowOffset = j * size;
            final int kmax = Math.min(j, i);
            acc.clear();
            for (int k = 0; k < kmax; k++) {
                rat.setFrom(LU, rowOffset + k, errorLU);
                rat.multiply(column, k, size);
                acc.add(rat);
            }
            acc.subtract(column, j, size);
            acc.negate();
            acc.storeTo(column, j, size);
            acc.storeTo(LU, rowOffset + i, errorLU);
        }
        /*
             * Find pivot and exchange if necessary. There is no floating-point arithmetic here
             * (ignoring the comparison for magnitude order), only work on index values.
             */
        int p = i;
        for (int j = i; ++j < size; ) {
            if (Math.abs(column[j]) > Math.abs(column[p])) {
                p = j;
            }
        }
        if (p != i) {
            final int pRow = p * size;
            final int iRow = i * size;
            for (int k = 0; k < size; k++) {
                // Swap two full rows.
                DoubleDouble.swap(LU, pRow + k, iRow + k, errorLU);
            }
            ArraysExt.swap(pivot, p, i);
        }
        /*
             * Compute multipliers. This part is equivalent to the following code, but
             * using double-double arithmetic instead than the primitive 'double' type:
             *
             *     final double sum = LU[i*size + i];
             *     if (sum != 0.0) {
             *         for (int j=i; ++j < size;) {
             *             LU[j*size + i] /= sum;
             *         }
             *     }
             */
        acc.setFrom(LU, i * size + i, errorLU);
        if (!acc.isZero()) {
            for (int j = i; ++j < size; ) {
                final int t = j * size + i;
                rat.setFrom(acc);
                rat.inverseDivide(LU, t, errorLU);
                rat.storeTo(LU, t, errorLU);
            }
        }
    }
    /*
         * At this point, we are done computing LU.
         * Ensure that the matrix is not singular.
         */
    for (int j = 0; j < size; j++) {
        rat.setFrom(LU, j * size + j, errorLU);
        if (rat.isZero()) {
            throw new NoninvertibleMatrixException(Resources.format(Resources.Keys.SingularMatrix));
        }
    }
    /*
         * Copy right hand side with pivoting. Write the result directly in the elements array
         * of the result matrix. This block does not perform floating-point arithmetic operations.
         */
    final GeneralMatrix result = GeneralMatrix.createExtendedPrecision(size, innerSize, false);
    final double[] elements = result.elements;
    final int errorOffset = size * innerSize;
    for (int k = 0, j = 0; j < size; j++) {
        final int p = pivot[j];
        for (int i = 0; i < innerSize; i++) {
            if (eltY != null) {
                final int t = p * innerSize + i;
                elements[k] = eltY[t];
                elements[k + errorOffset] = eltY[t + errorOffset];
            } else {
                elements[k] = Y.getElement(p, i);
            }
            k++;
        }
    }
    /*
         * Solve L*Y = B(pivot, :). The inner block is equivalent to the following line,
         * but using double-double arithmetic instead of 'double' primitive type:
         *
         *     elements[loRowOffset + i] -= (elements[rowOffset + i] * LU[luRowOffset + k]);
         */
    for (int k = 0; k < size; k++) {
        // Offset of row computed by current iteration.
        final int rowOffset = k * innerSize;
        for (int j = k; ++j < size; ) {
            // Offset of some row after the current row.
            final int loRowOffset = j * innerSize;
            // Offset of the corresponding row in the LU matrix.
            final int luRowOffset = j * size;
            for (int i = 0; i < innerSize; i++) {
                acc.setFrom(elements, loRowOffset + i, errorOffset);
                rat.setFrom(elements, rowOffset + i, errorOffset);
                rat.multiply(LU, luRowOffset + k, errorLU);
                acc.subtract(rat);
                acc.storeTo(elements, loRowOffset + i, errorOffset);
            }
        }
    }
    /*
         * Solve U*X = Y. The content of the loop is equivalent to the following line,
         * but using double-double arithmetic instead of 'double' primitive type:
         *
         *     double sum = LU[k*size + k];
         *     for (int i=0; i<innerSize; i++) {
         *         elements[rowOffset + i] /= sum;
         *     }
         *     for (int j=0; j<k; j++) {
         *         sum = LU[j*size + k];
         *         for (int i=0; i<innerSize; i++) {
         *             elements[upRowOffset + i] -= (elements[rowOffset + i] * sum);
         *         }
         *     }
         */
    for (int k = size; --k >= 0; ) {
        // Offset of row computed by current iteration.
        final int rowOffset = k * innerSize;
        // A diagonal element on the current row.
        acc.setFrom(LU, k * size + k, errorLU);
        for (int i = 0; i < innerSize; i++) {
            // Apply to all columns in the current row.
            rat.setFrom(acc);
            rat.inverseDivide(elements, rowOffset + i, errorOffset);
            rat.storeTo(elements, rowOffset + i, errorOffset);
        }
        for (int j = 0; j < k; j++) {
            // Offset of a row before (locate upper) the current row.
            final int upRowOffset = j * innerSize;
            // Same column than the diagonal element, but in the upper row.
            acc.setFrom(LU, j * size + k, errorLU);
            for (int i = 0; i < innerSize; i++) {
                // Apply to all columns in the upper row.
                rat.setFrom(elements, rowOffset + i, errorOffset);
                rat.multiply(acc);
                rat.subtract(elements, upRowOffset + i, errorOffset);
                rat.negate();
                rat.storeTo(elements, upRowOffset + i, errorOffset);
            }
        }
    }
    return result;
}

Example 30 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class Initializer method rν2.

/**
 * Computes the square of the reciprocal of the radius of curvature of the ellipsoid
 * perpendicular to the meridian at latitude φ. That radius of curvature is:
 *
 * <blockquote>ν = 1 / √(1 - ℯ²⋅sin²φ)</blockquote>
 *
 * This method returns 1/ν², which is the (1 - ℯ²⋅sin²φ) part of above equation.
 *
 * <div class="section">Relationship with Snyder</div>
 * This is related to functions (14-15) from Snyder (used for computation of scale factors
 * at the true scale latitude) as below:
 *
 * <blockquote>m = cosφ / sqrt(rν²)</blockquote>
 *
 * Special cases:
 * <ul>
 *   <li>If φ is 0°, then <var>m</var> is 1.</li>
 *   <li>If φ is ±90°, then <var>m</var> is 0 provided that we are not in the spherical case
 *       (otherwise we get {@link Double#NaN}).</li>
 * </ul>
 *
 * @param  sinφ  the sine of the φ latitude.
 * @return reciprocal squared of the radius of curvature of the ellipsoid
 *         perpendicular to the meridian at latitude φ.
 */
private DoubleDouble rν2(final double sinφ) {
    if (DoubleDouble.DISABLED) {
        return verbatim(1 - eccentricitySquared.value * (sinφ * sinφ));
    }
    final DoubleDouble t = verbatim(sinφ);
    t.square();
    t.multiply(eccentricitySquared);
    /*
         * Compute 1 - ℯ²⋅sin²φ.  Since  ℯ²⋅sin²φ  may be small,
         * this is where double-double arithmetic has more value.
         */
    t.negate();
    t.add(1, 0);
    return t;
}