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

Example 11 with DoubleDouble

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

the class MatrixSIS method setMatrix.

/**
 * Sets this matrix to the values of another matrix.
 * The given matrix must have the same size.
 *
 * @param  matrix  the matrix to copy.
 * @throws MismatchedMatrixSizeException if the given matrix has a different size than this matrix.
 *
 * @since 0.7
 */
public void setMatrix(final Matrix matrix) throws MismatchedMatrixSizeException {
    ArgumentChecks.ensureNonNull("matrix", matrix);
    final int numRow = getNumRow();
    final int numCol = getNumCol();
    ensureSizeMatch(numRow, numCol, matrix);
    final int count = numRow * numCol;
    final double[] elements;
    /*
         * If both matrices use extended precision, the elements array will have twice the expected length
         * with the matrix values in the first half and the error terms in the second half.  If we want to
         * preserve the extended precision, we have to transfer the values between the two matrices with a
         * DoubleDouble object.
         */
    if (isExtendedPrecision() && matrix instanceof ExtendedPrecisionMatrix) {
        elements = ((ExtendedPrecisionMatrix) matrix).getExtendedElements();
        if (elements.length > count) {
            final DoubleDouble t = new DoubleDouble();
            for (int i = 0; i < count; i++) {
                t.value = elements[i];
                t.error = elements[i + count];
                set(i / numCol, i % numCol, t);
            }
            return;
        }
    } else {
        // Fallback for matrices that do not use extended precision.
        elements = new double[count];
        getElements(matrix, numRow, numCol, elements);
    }
    setElements(elements);
}

Example 12 with DoubleDouble

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

the class MatrixSIS method removeRows.

/**
 * Returns a new matrix with the same elements than this matrix except for the specified rows.
 * This method is useful for removing a range of <em>target</em> dimensions in an affine transform.
 *
 * @param  lower  index of the first row to remove (inclusive).
 * @param  upper  index after the last row to remove (exclusive).
 * @return a copy of this matrix with the specified rows removed.
 *
 * @since 0.7
 */
public MatrixSIS removeRows(final int lower, final int upper) {
    final int numRow = getNumRow();
    final int numCol = getNumCol();
    ArgumentChecks.ensureValidIndexRange(numRow, lower, upper);
    final DoubleDouble dd = isExtendedPrecision() ? new DoubleDouble() : null;
    final MatrixSIS reduced = Matrices.createZero(numRow - (upper - lower), numCol, dd != null);
    int dest = 0;
    for (int j = 0; j < numRow; j++) {
        if (j == lower) {
            j = upper;
            if (j == numRow)
                break;
        }
        for (int i = 0; i < numCol; i++) {
            if (dd != null) {
                get(j, i, dd);
                reduced.set(dest, i, dd);
            } else {
                reduced.setElement(dest, i, getElement(j, i));
            }
        }
        dest++;
    }
    return reduced;
}

Example 13 with DoubleDouble

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

the class Initializer method axisLengthRatio.

/**
 * Returns {@code b/a} where {@code a} is the semi-major axis length and {@code b} the semi-minor axis length.
 * We retrieve this value from the eccentricity with {@code b/a = sqrt(1-ℯ²)}.
 */
final DoubleDouble axisLengthRatio() {
    final DoubleDouble b = new DoubleDouble(1, 0);
    b.subtract(eccentricitySquared);
    b.sqrt();
    return b;
}

Example 14 with DoubleDouble

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

the class CoordinateSystems method swapAndScaleAxes.

/**
 * Returns an affine transform between two coordinate systems.
 * Only units and axes order (e.g. transforming from
 * ({@linkplain AxisDirection#NORTH North}, {@linkplain AxisDirection#WEST West}) to
 * ({@linkplain AxisDirection#EAST East}, {@linkplain AxisDirection#NORTH North})
 * are taken in account by this method.
 *
 * <div class="section">Conditions</div>
 * The two coordinate systems must implement the same GeoAPI coordinate system interface.
 * For example if {@code sourceCS} is a {@link org.opengis.referencing.cs.CartesianCS},
 * then {@code targetCS} must be a {@code CartesianCS} too.
 *
 * <div class="note"><b>Example:</b>
 * If coordinates in {@code sourceCS} are (<var>x</var>,<var>y</var>) tuples in metres
 * and coordinates in {@code targetCS} are (<var>-y</var>,<var>x</var>) tuples in centimetres,
 * then the transformation can be performed as below:
 *
 * {@preformat math
 *     ┌      ┐   ┌                ┐ ┌     ┐
 *     │-y(cm)│   │   0  -100    0 │ │ x(m)│
 *     │ x(cm)│ = │ 100     0    0 │ │ y(m)│
 *     │ 1    │   │   0     0    1 │ │ 1   │
 *     └      ┘   └                ┘ └     ┘
 * }
 * </div>
 *
 * @param  sourceCS  the source coordinate system.
 * @param  targetCS  the target coordinate system.
 * @return the conversion from {@code sourceCS} to {@code targetCS} as an affine transform.
 *         Only axis direction and units are taken in account.
 * @throws IllegalArgumentException if the CS are not of the same type, or axes do not match.
 * @throws IncommensurableException if the units are not compatible, or the conversion is non-linear.
 *
 * @see Matrices#createTransform(AxisDirection[], AxisDirection[])
 */
@SuppressWarnings("fallthrough")
public static Matrix swapAndScaleAxes(final CoordinateSystem sourceCS, final CoordinateSystem targetCS) throws IllegalArgumentException, IncommensurableException {
    ensureNonNull("sourceCS", sourceCS);
    ensureNonNull("targetCS", targetCS);
    if (!Classes.implementSameInterfaces(sourceCS.getClass(), targetCS.getClass(), CoordinateSystem.class)) {
        throw new IllegalArgumentException(Resources.format(Resources.Keys.IncompatibleCoordinateSystemTypes));
    }
    final AxisDirection[] srcAxes = getAxisDirections(sourceCS);
    final AxisDirection[] dstAxes = getAxisDirections(targetCS);
    final MatrixSIS matrix = Matrices.createTransform(srcAxes, dstAxes);
    assert Arrays.equals(srcAxes, dstAxes) == matrix.isIdentity() : matrix;
    /*
         * The previous code computed a matrix for swapping axes. Usually, this
         * matrix contains only 0 and 1 values with only one "1" value by row.
         * For example, the matrix operation for swapping x and y axes is:
         *          ┌ ┐   ┌         ┐ ┌ ┐
         *          │y│   │ 0  1  0 │ │x│
         *          │x│ = │ 1  0  0 │ │y│
         *          │1│   │ 0  0  1 │ │1│
         *          └ ┘   └         ┘ └ ┘
         * Now, take in account units conversions. Each matrix's element (j,i)
         * is multiplied by the conversion factor from sourceCS.getUnit(i) to
         * targetCS.getUnit(j). This is an element-by-element multiplication,
         * not a matrix multiplication. The last column is processed in a special
         * way, since it contains the offset values.
         */
    // == sourceCS.getDimension()
    final int sourceDim = matrix.getNumCol() - 1;
    // == targetCS.getDimension()
    final int targetDim = matrix.getNumRow() - 1;
    for (int j = 0; j < targetDim; j++) {
        final Unit<?> targetUnit = targetCS.getAxis(j).getUnit();
        for (int i = 0; i < sourceDim; i++) {
            if (matrix.getElement(j, i) == 0) {
                // (i.e. axes are orthogonal).
                continue;
            }
            final Unit<?> sourceUnit = sourceCS.getAxis(i).getUnit();
            if (Objects.equals(sourceUnit, targetUnit)) {
                // between source[i] and target[j].
                continue;
            }
            Number scale = 1;
            Number offset = 0;
            final Number[] coefficients = Units.coefficients(sourceUnit.getConverterToAny(targetUnit));
            switch(coefficients != null ? coefficients.length : -1) {
                // Fall through
                case 2:
                    scale = coefficients[1];
                // Fall through
                case 1:
                    offset = coefficients[0];
                case 0:
                    break;
                default:
                    throw new IncommensurableException(Resources.format(Resources.Keys.NonLinearUnitConversion_2, sourceUnit, targetUnit));
            }
            final DoubleDouble element = DoubleDouble.castOrCopy(matrix.getNumber(j, i));
            final DoubleDouble r = new DoubleDouble(element);
            r.multiply(scale);
            matrix.setNumber(j, i, r);
            r.setFrom(element);
            r.multiply(offset);
            r.add(matrix.getNumber(j, sourceDim));
            matrix.setNumber(j, sourceDim, r);
        }
    }
    return matrix;
}

Example 15 with DoubleDouble

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

the class Line method fit.

/**
 * Given a sequence of points, fits them to a straight line
 * <var>y</var> = <var>slope</var>⋅<var>x</var> + <var>y₀</var> in a least-squares senses.
 * Points shall be two dimensional with ordinate values in the (<var>x</var>,<var>y</var>) order.
 * This method assumes that the <var>x</var> values are precise and all uncertainty is in <var>y</var>.
 * {@link Double#NaN} ordinate values are ignored.
 *
 * @param  points  the two-dimensional points.
 * @return estimation of the correlation coefficient. The closer this coefficient is to +1 or -1, the better the fit.
 * @throws MismatchedDimensionException if a point is not two-dimensional.
 */
public double fit(final Iterable<? extends DirectPosition> points) {
    ArgumentChecks.ensureNonNull("points", points);
    int i = 0, n = 0;
    final DoubleDouble mean_x = new DoubleDouble();
    final DoubleDouble mean_y = new DoubleDouble();
    for (final DirectPosition p : points) {
        final int dimension = p.getDimension();
        if (dimension != DIMENSION) {
            throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimension_3, "points[" + i + ']', DIMENSION, dimension));
        }
        i++;
        final double x, y;
        if (// Test first the dimension which is most likely to contain NaN.
        !isNaN(y = p.getOrdinate(1)) && !isNaN(x = p.getOrdinate(0))) {
            mean_x.add(x);
            mean_y.add(y);
            n++;
        }
    }
    mean_x.divide(n, 0);
    mean_y.divide(n, 0);
    /*
         * We have to solve two equations with two unknowns:
         *
         *   1)    mean(y)   = m⋅mean(x) + y₀
         *   2)    mean(x⋅y) = m⋅mean(x²) + y₀⋅mean(x)
         *
         * Those formulas lead to a quadratic equation. However, the formulas become very simples
         * if we set 'mean(x) = 0'. We can achieve this result by computing instead of (2):
         *
         *   2b)   mean(Δx⋅y) = m⋅mean(Δx²)
         *
         * where dx = x-mean(x). In this case mean(Δx) == 0.
         */
    final DoubleDouble a = new DoubleDouble();
    final DoubleDouble b = new DoubleDouble();
    final DoubleDouble mean_x2 = new DoubleDouble();
    final DoubleDouble mean_y2 = new DoubleDouble();
    final DoubleDouble mean_xy = new DoubleDouble();
    for (final DirectPosition p : points) {
        final double y;
        if (// Test first the dimension which is most likely to contain NaN.
        !isNaN(y = p.getOrdinate(1)) && !isNaN(a.value = p.getOrdinate(0))) {
            a.error = 0;
            // Δx = x - mean_x
            a.subtract(mean_x);
            b.setFrom(a);
            // (Δx)² = Δx * Δx
            b.multiply(a);
            // mean_x² += (Δx)²
            mean_x2.add(b);
            // xy = Δx * y
            a.multiply(y);
            // mean_xy += xy
            mean_xy.add(a);
            // y² = y * y
            a.setToProduct(y, y);
            // mean_y² += y²
            mean_y2.add(a);
        }
    }
    mean_x2.divide(n, 0);
    mean_y2.divide(n, 0);
    mean_xy.divide(n, 0);
    /*
         * Assuming that 'mean(x) == 0', then the correlation
         * coefficient can be approximate by:
         *
         * R = mean(xy) / sqrt( mean(x²) * (mean(y²) - mean(y)²) )
         */
    a.setFrom(mean_xy);
    // slope = mean_xy / mean_x²
    a.divide(mean_x2);
    b.setFrom(mean_x);
    b.multiply(a);
    b.negate();
    // y₀ = mean_y - mean_x * slope
    b.add(mean_y);
    setEquation(a, b);
    /*
         * Compute the correlation coefficient:
         * mean_xy / sqrt(mean_x2 * (mean_y2 - mean_y * mean_y))
         */
    a.setFrom(mean_y);
    a.multiply(mean_y);
    a.negate();
    a.add(mean_y2);
    a.multiply(mean_x2);
    a.sqrt();
    a.inverseDivide(mean_xy);
    return a.value;
}