sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject svn commit: r1530255 - in /sis/branches/JDK7/core: sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/ sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/ sis-utility/src/main/java/org/apache/sis/util/resou...
Date Tue, 08 Oct 2013 12:52:33 GMT
Author: desruisseaux
Date: Tue Oct  8 12:52:33 2013
New Revision: 1530255

URL: http://svn.apache.org/r1530255
Log:
Special case for inversion of matrix having less columns than rows.

Modified:
    sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrix.java
    sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Solver.java
    sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrixTest.java
    sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java
    sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties
    sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties

Modified: sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrix.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrix.java?rev=1530255&r1=1530254&r2=1530255&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrix.java
[UTF-8] (original)
+++ sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrix.java
[UTF-8] Tue Oct  8 12:52:33 2013
@@ -17,6 +17,7 @@
 package org.apache.sis.referencing.operation.matrix;
 
 import org.opengis.referencing.operation.Matrix;
+import org.apache.sis.util.resources.Errors;
 
 
 /**
@@ -146,105 +147,182 @@ final class NonSquareMatrix extends Gene
      * There is an issue about whether the full row shall contains NaN, or only the last
element (the translation
      * term) as in the above example.  The current implementation inserts a NaN value in
the translation term and
      * sets all other values to 0 on the assumption that if (x₂,y₂) do not depend on
(z,t), then conversely (z,t)
-     * do not depend on (x₂,y₂) neither.
+     * do not depend on (x₂,y₂) neither. Setting the scale factor to zero expresses that
fact, while setting them
+     * to NaN would mean "don't know".
+     *
+     * <p>Note that the above strategy assumes that the matrix is used for an affine
transform, which is not always
+     * true (it could be the matrix of a map projection derivative for instance). If the
matrix is not for an affine
+     * transform, then the last column has no special meaning and the above strategy is somewhat
asymmetric.
+     * However it will still produce NaN for the full row in matrix multiplications.</p>
+     *
+     * <p>Conversely, if the matrix has more rows than columns (in a system of linear
equations, the system would
+     * be <cite>overdetermined</cite>), then we omit the rows containing only
zero or NaN values. After the matrix
+     * inversion, we insert columns having only zero values for the dimensions associated
to those rows.
+     * Semantically, the inverse matrix is a (x₁,y₁,z,t) → (x₂,y₂) transform that
just discards the ordinate values
+     * at the dimensions corresponding to those rows.</p>
      */
     @Override
     public MatrixSIS inverse() throws NoninvertibleMatrixException {
+        if (numRow < numCol) {
+            return inverseDimensionReduction();
+        } else {
+            return inverseDimensionIncrease();
+        }
+    }
+
+    /**
+     * Inverses a matrix for a transform where target points have fewer ordinates than source
points.
+     * If a column contains only zero values, then this means that the ordinate at the corresponding
+     * column is simply deleted. We can omit that column. We check the last columns before
the first
+     * columns on the assumption that last dimensions are more likely to be independent dimensions
+     * like time.
+     */
+    private MatrixSIS inverseDimensionReduction() throws NoninvertibleMatrixException {
         final int numRow = this.numRow; // Protection against accidental changes.
         final int numCol = this.numCol;
-        if (numRow < numCol) {
-            final int length = numRow * numCol;
-            /*
-             * Target points have fewer ordinates than source points. If a column contains
only zero values,
-             * then this means that the ordinate at the corresponding column is simply deleted.
We can omit
-             * that column. We check the last columns before the first columns on the assumption
that last
-             * dimensions are more likely to be independant dimensions like time.
-             */
-            int oi = numCol - numRow;
-            final int[] omitted = new int[oi];
-skipColumn: for (int i=numCol; --i>=0;) {
-                for (int j=length + i; (j -= numCol) >= 0;) {
-                    if (elements[j] != 0) {
-                        continue skipColumn;
-                    }
+        final int length = numRow * numCol;
+        int i  = numCol;
+        int oi = numCol - numRow; // Initialized to the maximal amount of columns that we
may omit.
+        final int[] omitted = new int[oi];
+next:   do {
+            if (--i < 0) {
+                throw nonInvertible(); // Not enough columns that we can omit.
+            }
+            for (int j=length + i; (j -= numCol) >= 0;) {
+                if (elements[j] != 0) {
+                    continue next;
                 }
-                omitted[--oi] = i; // Found a column which contains only 0 elements.
-                if (oi == 0) {
-                    /*
-                     * Found enough columns containing only zero elements. Create a square
matrix omitting those
-                     * columns, and invert that matrix. Note that we also need to either
copy the error terms,
-                     * or to infer them.
-                     */
-                    GeneralMatrix squareMatrix = new GeneralMatrix(numRow, numRow, false,
2);
-                    int j=0;
-                    for (i=0; i<numCol; i++) {
-                        if (oi != omitted.length && i == omitted[oi]) oi++;
-                        else copyColumnTo(i, squareMatrix, j++); // Copy only if not skipped.
-                    }
-                    // If the source matrix does not use double-double arithmetic, infer
the error terms.
-                    if (indexOfErrors(numRow, numCol, elements) == 0) {
-                        inferErrors(squareMatrix.elements);
-                    }
-                    squareMatrix = (GeneralMatrix) Solver.inverse(squareMatrix, false);
-                    /*
-                     * Create a new matrix with new rows added for the omitted ordinates.
-                     * From this point, the meaning of 'numCol' and 'numRow' are interchanged.
-                     */
-                    final NonSquareMatrix inverse = new NonSquareMatrix(numCol, numRow, false,
2);
-                    for (oi=0, j=0, i=0; i<numCol; i++) {
-                        if (oi != omitted.length && i == omitted[oi]) {
-                            inverse.setElement(i, numRow-1, Double.NaN);
-                            oi++;
-                        } else {
-                            inverse.copyRowFrom(squareMatrix, j++, i);
-                        }
-                    }
-                    return inverse;
+            }
+            omitted[--oi] = i; // Found a column which contains only 0 elements.
+        } while (oi != 0);
+        /*
+         * Found enough columns containing only zero elements. Create a square matrix omitting
those columns,
+         * and invert that matrix. Note that we also need to either copy the error terms,
or to infer them.
+         */
+        GeneralMatrix squareMatrix = new GeneralMatrix(numRow, numRow, false, 2);
+        int j = 0;
+        for (i=0; i<numCol; i++) {
+            if (oi != omitted.length && i == omitted[oi]) oi++;
+            else copyColumn(this, i, squareMatrix, j++); // Copy only if not skipped.
+        }
+        // If the source matrix does not use double-double arithmetic, infer the error terms.
+        if (indexOfErrors(numRow, numCol, elements) == 0) {
+            inferErrors(squareMatrix.elements);
+        }
+        squareMatrix = (GeneralMatrix) Solver.inverse(squareMatrix, false);
+        /*
+         * Create a new matrix with new rows added for the omitted ordinates.
+         * From this point, the meaning of 'numCol' and 'numRow' are interchanged.
+         */
+        final NonSquareMatrix inverse = new NonSquareMatrix(numCol, numRow, false, 2);
+        for (oi=0, j=0, i=0; i<numCol; i++) {
+            if (oi != omitted.length && i == omitted[oi]) {
+                inverse.setElement(i, numRow-1, Double.NaN); // Translation term to NaN,
remaining to 0.
+                oi++;
+            } else {
+                copyRow(squareMatrix, j++, inverse, i);
+            }
+        }
+        return inverse;
+    }
+
+    /**
+     * Inverse a matrix for a transform where target points has more ordinates than source
points.
+     * In other words, the target matrices will be a transform that discard some ordinate
values.
+     * We will discard the ones for which the row contains only 0 or NaN elements.
+     */
+    private MatrixSIS inverseDimensionIncrease() throws NoninvertibleMatrixException {
+        final int numRow = this.numRow; // Protection against accidental changes.
+        final int numCol = this.numCol;
+        int j  = numRow;
+        int oi = numRow - numCol; // Initialized to the maximal amount of rows that we may
discard.
+        final int[] omitted = new int[oi];
+next:   do {
+            if (--j < 0) {
+                throw nonInvertible(); // Not enough rows that we can omit.
+            }
+            final int offset = j * numCol;
+            for (int i=offset + numCol; --i >= offset;) {
+                final double element = elements[i];
+                if (element != 0 && !Double.isNaN(element)) {
+                    continue next;
                 }
             }
+            omitted[--oi] = j; // Found a row which contains only 0 or NaN elements.
+        } while (oi != 0);
+        /*
+         * Found enough rows containing only zero elements. Create a square matrix omitting
those rows,
+         * and invert that matrix. Note that we also need to either copy the error terms,
or to infer them.
+         */
+        GeneralMatrix squareMatrix = new GeneralMatrix(numCol, numCol, false, 2);
+        int i = 0;
+        for (j=0; j<numRow; j++) {
+            if (oi != omitted.length && j == omitted[oi]) oi++;
+            else copyRow(this, j, squareMatrix, i++); // Copy only if not skipped.
         }
+        if (indexOfErrors(numRow, numCol, elements) == 0) {
+            inferErrors(squareMatrix.elements);
+        }
+        squareMatrix = (GeneralMatrix) Solver.inverse(squareMatrix, false);
         /*
-         * If we reach this point, we have not been able to replace the non-square matrix
by a square one.
-         * Delegate to the super-class method as a matter of principle, but that method is
expected to fail.
+         * Create a new matrix with new columns added for the omitted rows.
+         * From this point, the meaning of 'numCol' and 'numRow' are interchanged.
          */
-        return super.inverse();
+        final NonSquareMatrix inverse = new NonSquareMatrix(numCol, numRow, false, 2);
+        for (oi=0, i=0, j=0; j<numRow; j++) {
+            if (oi != omitted.length && j == omitted[oi]) oi++;
+            else copyColumn(squareMatrix, i++, inverse, j);
+        }
+        return inverse;
     }
 
     /**
-     * Copies a column from this matrix to the given matrix, including the double-double
arithmetic error terms
-     * if any. The given matrix must have the same number of rows than this matrix, and must
have enough room for
-     * error terms (this is not verified).
+     * Copies a column between two matrices, including the double-double arithmetic error
terms if any.
+     * The target matrix must have the same number of rows than the source matrix, and must
have enough
+     * room for error terms (this is not verified).
      *
-     * @param srcIndex  Index of the column to copy from this matrix.
+     * @param source    The matrix from which to copy a column.
+     * @param srcIndex  Index of the column to copy from the source matrix.
      * @param target    The matrix where to copy the column.
      * @param dstIndex  Index of the column where to copy in the target matrix.
      */
-    private void copyColumnTo(int srcIndex, final GeneralMatrix target, int dstIndex) {
-        assert target.numRow == numRow;
-        while (srcIndex < elements.length) {
-            target.elements[dstIndex] = elements[srcIndex];
-            srcIndex += this.numCol;
+    private static void copyColumn(final GeneralMatrix source, int srcIndex, final GeneralMatrix
target, int dstIndex) {
+        assert target.numRow == source.numRow;
+        while (srcIndex < source.elements.length) {
+            target.elements[dstIndex] = source.elements[srcIndex];
+            srcIndex += source.numCol;
             dstIndex += target.numCol;
         }
     }
 
     /**
-     * Copies a row from the given matrix to this matrix, including the double-double arithmetic
error terms.
-     * The two matrix must have the same number of columns, and both of them must have room
for the error terms
-     * (this is not verified).
-     *
-     * @param source   The matrix from which to copy a row.
-     * @param srcIndex Index of the row to copy from the source matrix.
-     * @param dstIndex Index of the row where to copy in this matrix.
+     * Copies a row between two matrices, including the double-double arithmetic error terms
if any.
+     * The target matrix must have the same number of columns than the source matrix, and
must have
+     * enough room for error terms (this is not verified).
+     *
+     * @param source    The matrix from which to copy a row.
+     * @param srcIndex  Index of the row to copy from the source matrix.
+     * @param target    The matrix where to copy the row.
+     * @param dstIndex  Index of the row where to copy in the target matrix.
      */
-    private void copyRowFrom(final GeneralMatrix source, int srcIndex, int dstIndex) {
-        final int numCol = this.numCol;
+    private static void copyRow(final GeneralMatrix source, int srcIndex, final GeneralMatrix
target, int dstIndex) {
+        final int numCol = target.numCol;
         assert numCol == source.numCol;
         srcIndex *= numCol;
         dstIndex *= numCol;
-        System.arraycopy(source.elements, srcIndex, elements, dstIndex, numCol);  // Copy
main values
-        System.arraycopy(source.elements, srcIndex + numCol * source.numRow,      // Copy
error terms
-                                elements, dstIndex + numCol * numRow, numCol);
+        System.arraycopy(source.elements, srcIndex, target.elements, dstIndex, numCol);
+        srcIndex += numCol * source.numRow;
+        if (srcIndex < source.elements.length) {
+            dstIndex += numCol * target.numRow;
+            System.arraycopy(source.elements, srcIndex, target.elements, dstIndex, numCol);
+        }
+    }
+
+    /**
+     * Returns the exception for a non-invertible transform.
+     */
+    private NoninvertibleMatrixException nonInvertible() {
+        return new NoninvertibleMatrixException(Errors.format(Errors.Keys.NonInvertibleMatrix_2,
numRow, numCol));
     }
 
     /**

Modified: sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Solver.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Solver.java?rev=1530255&r1=1530254&r2=1530255&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Solver.java
[UTF-8] (original)
+++ sis/branches/JDK7/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Solver.java
[UTF-8] Tue Oct  8 12:52:33 2013
@@ -403,8 +403,7 @@ searchNaN:  for (int flatIndex = (size -
         for (int j=0; j<size; j++) {
             rat.setFrom(LU, j*size + j, errorLU);
             if (rat.isZero()) {
-                final Integer n = size;
-                throw new NoninvertibleMatrixException(Errors.format(Errors.Keys.NonInvertibleMatrix_2,
n, n));
+                throw new NoninvertibleMatrixException(Errors.format(Errors.Keys.SingularMatrix));
             }
         }
         /*

Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrixTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrixTest.java?rev=1530255&r1=1530254&r2=1530255&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrixTest.java
[UTF-8] (original)
+++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/NonSquareMatrixTest.java
[UTF-8] Tue Oct  8 12:52:33 2013
@@ -79,6 +79,7 @@ public final strictfp class NonSquareMat
     @Override
     public void testInverse() throws NoninvertibleMatrixException {
         testDimensionReduction(null, 1, 0);
+        testDimensionIncrease (null, 1);
     }
 
     /**
@@ -86,11 +87,16 @@ public final strictfp class NonSquareMat
      */
     @Override
     public void testSolve() throws NoninvertibleMatrixException {
-        final Matrix3 Y = new Matrix3(
+        testDimensionReduction(new Matrix3(
                 2, 0, 0,
                 0, 2, 0,
-                0, 0, 1);
-        testDimensionReduction(Y, 2, NaN);
+                0, 0, 1), 2, NaN);
+        testDimensionIncrease(new GeneralMatrix(5, 5, new double[] {
+                2, 0, 0, 0, 0,
+                0, 2, 0, 0, 0,
+                0, 0, 2, 0, 0,
+                0, 0, 0, 2, 0,
+                0, 0, 0, 0, 1}), 2);
     }
 
     /**
@@ -111,14 +117,41 @@ public final strictfp class NonSquareMat
             0, 0, 0, 0, 1
         });
         final double[] expected = {
-            0.5*sf, 0,          -4,
-            uks,    uks,       NaN,
-            0,      0.25*sf, -1.25,
-            uks,    uks,       NaN,
-            0,      0,           1
+            0.5*sf,  0,           -4,
+            uks,     uks,        NaN,
+            0,       0.25*sf,  -1.25,
+            uks,     uks,        NaN,
+            0,       0,            1
         };
         final MatrixSIS inverse = (Y != null) ? matrix.solve(Y) : matrix.inverse();
-        assertMatrixEquals(expected, 5, 3, inverse, SolverTest.TOLERANCE);
+        assertMatrixEquals(expected, 5, 3, inverse, TOLERANCE);
+    }
+
+    /**
+     * Tests {@link NonSquareMatrix#inverse()} or {@link NonSquareMatrix#solve(Matrix)} with
a conversion
+     * matrix having more target dimensions (rows) than source dimensions (columns).
+     *
+     * @param  Y    The matrix to give to {@code solve(Y)}, {@code null} for testing {@code
inverse()}.
+     * @param  sf   The scale factor by which to multiply all expected scale elements.
+     * @throws NoninvertibleMatrixException Should never happen.
+     */
+    private static void testDimensionIncrease(final MatrixSIS Y, final double sf)
+            throws NoninvertibleMatrixException
+    {
+        final MatrixSIS matrix = Matrices.create(5, 3, new double[] {
+              2,   0,   8,
+            NaN, NaN, NaN,
+              0,   4,   5,
+              0,   0,   0,
+              0,   0,   1
+        });
+        final double[] expected = {
+            0.5*sf,  0,  0,        0,  -4,
+            0,       0,  0.25*sf,  0,  -1.25,
+            0,       0,  0,        0,   1
+        };
+        final MatrixSIS inverse = (Y != null) ? matrix.solve(Y) : matrix.inverse();
+        assertMatrixEquals(expected, 3, 5, inverse, TOLERANCE);
     }
 
     /**

Modified: sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java?rev=1530255&r1=1530254&r2=1530255&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java
[UTF-8] (original)
+++ sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.java
[UTF-8] Tue Oct  8 12:52:33 2013
@@ -536,6 +536,11 @@ public final class Errors extends Indexe
         public static final int RequireDecimalSeparator = 33;
 
         /**
+         * Matrix is singular.
+         */
+        public static final int SingularMatrix = 125;
+
+        /**
          * Thread “{0}” seems stalled.
          */
         public static final int StalledThread_1 = 63;

Modified: sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties?rev=1530255&r1=1530254&r2=1530255&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties
[ISO-8859-1] (original)
+++ sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors.properties
[ISO-8859-1] Tue Oct  8 12:52:33 2013
@@ -117,6 +117,7 @@ NullArgument_1                  = Argume
 NullMapKey                      = Null key is not allowed in this dictionary.
 NullMapValue                    = Null values are not allowed in this dictionary.
 OddArrayLength_1                = Array length is {0}, while we expected an even length.
+SingularMatrix                  = Matrix is singular.
 RecursiveCreateCallForKey_1     = Recursive call while creating an object for the \u201c{0}\u201d
key.
 RequireDecimalSeparator         = A decimal separator is required.
 StalledThread_1                 = Thread \u201c{0}\u201d seems stalled.

Modified: sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties?rev=1530255&r1=1530254&r2=1530255&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties
[ISO-8859-1] (original)
+++ sis/branches/JDK7/core/sis-utility/src/main/java/org/apache/sis/util/resources/Errors_fr.properties
[ISO-8859-1] Tue Oct  8 12:52:33 2013
@@ -108,6 +108,7 @@ NullMapValue                    = Les va
 OddArrayLength_1                = La longueur du tableau est {0}, alors qu\u2019on attendait
une longueur paire.
 RecursiveCreateCallForKey_1     = Appel r\u00e9cursif lors de la cr\u00e9ation d\u2019un
objet pour la cl\u00e9 \u201c{0}\u201d.
 RequireDecimalSeparator         = Un s\u00e9parateur d\u00e9cimal est requis.
+SingularMatrix                  = La matrice est singuli\u00e8re.
 StalledThread_1                 = La t\u00e2che \u201c{0}\u201d semble bloqu\u00e9e.
 StreamIsForwardOnly_1           = Ne peut pas reculer dans le flux de donn\u00e9es \u201c{0}\u201d.
 TooFewArguments_2               = Au moins {0} argument{0,choice,1# \u00e9tait attendu|2#s
\u00e9taient attendus}, mais seulement {1} {1,choice,1#a \u00e9t\u00e9 sp\u00e9cifi\u00e9|2#ont
\u00e9t\u00e9 sp\u00e9cifi\u00e9s}.



Mime
View raw message