sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject svn commit: r1529310 - in /sis/branches/JDK7/core/sis-referencing/src: main/java/org/apache/sis/referencing/operation/matrix/Solver.java test/java/org/apache/sis/referencing/operation/matrix/SolverTest.java
Date Fri, 04 Oct 2013 21:32:43 GMT
Author: desruisseaux
Date: Fri Oct  4 21:32:42 2013
New Revision: 1529310

URL: http://svn.apache.org/r1529310
Log:
Ported from Geotk the handling of NaN values.

Modified:
    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/SolverTest.java

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=1529310&r1=1529309&r2=1529310&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] Fri Oct  4 21:32:42 2013
@@ -31,12 +31,20 @@ import org.apache.sis.util.ArraysExt;
  * <p>This class implements the {@link Matrix} interface as an implementation convenience.
  * This implementation details can be ignored.</p>
  *
+ * @author  JAMA team
+ * @author  Martin Desruisseaux (Geomatys)
  * @since   0.4
  * @version 0.4
  * @module
  */
 final class Solver implements Matrix {
     /**
+     * The size of the (i, j, s) tuples used internally by {@link #solve(MatrixSIS, Matrix,
double[], int, int)}
+     * for storing information about the NaN values.
+     */
+    private static final int TUPLE_SIZE = 3;
+
+    /**
      * A immutable identity matrix without defined size.
      * This is used only for computing the inverse.
      */
@@ -138,8 +146,12 @@ final class Solver implements Matrix {
     }
 
     /**
-     * Implementation of {@code solve} and {@code inverse} methods.
-     * Use a "left-looking", dot-product, Crout/Doolittle algorithm.
+     * Implementation of {@code solve} and {@code inverse} methods, with filtering of NaN
values.
+     * This method searches for NaN values before to attempt solving or inverting the matrix.
+     * If some NaN values are found but the matrix is written in such a way that each NaN
value
+     * is used for exactly one ordinate value (i.e. a matrix row is used for a one-dimensional
+     * conversion which is independent of all other dimensions), then we will edit the matrix
in
+     * such a way that this NaN value does not prevent the inverse matrix to be computed.
      *
      * <p>This method does <strong>not</strong> checks the matrix size.
      * Check for matrix size shall be performed by the caller like below:</p>
@@ -166,9 +178,139 @@ final class Solver implements Matrix {
     {
         assert (X.getNumRow() == size && X.getNumCol() == size) : size;
         assert (Y.getNumRow() == size && Y.getNumCol() == innerSize) || (Y instanceof
Solver);
+        final double[] LU = GeneralMatrix.getExtendedElements(X, size, size, true);
+        final int lastRowOrColumn = size - 1;
+        /*
+         * indexOfNaN array will be created only if at least one NaN value is found, and
those NaN meet
+         * the conditions documented in the code below. In such case, the array will contain
a sequence
+         * of (i,j,s) where (i,j) are the indices where the NaN value has been found and
s is the column
+         * of the scale factor.
+         */
+        int[] indexOfNaN = null;
+        int   indexCount = 0;
+        if (X.isAffine()) {
+            /*
+             * Conservatively search for NaN values only if the matrix looks like an affine
transform.
+             * If the matrix is affine, then we will assume that we can interpret the last
column as
+             * translation terms and other columns as scale factors.
+             */
+searchNaN:  for (int j=lastRowOrColumn; --j>=0;) {  // Skip the last row, since it is
(0, 0, ..., 1).
+                for (int i=size; --i>=0;) {         // Scan all columns, including the
translation terms.
+                    if (Double.isNaN(LU[j*size + i])) {
+                        /*
+                         * Found a NaN value. First, if we are not in the translation column,
ensure
+                         * that the column contains only zero values except on the current
line.
+                         */
+                        int columnOfScale = -1;
+                        if (i != lastRowOrColumn) {                // Enter only if this
column is not for translations.
+                            columnOfScale = i;                     // The non-translation
element is the scale factor.
+                            for (int k=lastRowOrColumn; --k>=0;) { // Scan all other rows
in the current column.
+                                if (k != j && LU[k*size + i] != 0) {
+                                    // Found a non-zero element in the current column.
+                                    // We can not proceed - cancel everything.
+                                    indexOfNaN = null;
+                                    indexCount = 0;
+                                    break searchNaN;
+                                }
+                            }
+                        }
+                        /*
+                         * Next, ensure that the row contains only zero elements except for
+                         * the scale factor and the offset (the element in the translation
+                         * column, which is not checked by the loop below).
+                         */
+                        for (int k=lastRowOrColumn; --k>=0;) {
+                            if (k != i && LU[j*size + k] != 0) {
+                                if (columnOfScale >= 0) {
+                                    // If there is more than 1 non-zero element,
+                                    // abandon the attempt to handle NaN values.
+                                    indexOfNaN = null;
+                                    indexCount = 0;
+                                    break searchNaN;
+                                }
+                                columnOfScale = k;
+                            }
+                        }
+                        /*
+                         * At this point, the NaN element has been determined as replaceable.
+                         * Remember its index; the replacement will be performed later.
+                         */
+                        if (indexOfNaN == null) {
+                            indexOfNaN = new int[lastRowOrColumn * (2*TUPLE_SIZE)]; // At
most one scale and one offset per row.
+                        }
+                        indexOfNaN[indexCount++] = i;
+                        indexOfNaN[indexCount++] = j;
+                        indexOfNaN[indexCount++] = columnOfScale; // May be -1 (while uncommon)
+                        assert (indexCount % TUPLE_SIZE) == 0;
+                    }
+                }
+            }
+            /*
+             * If there is any NaN value to edit, replace them by 0 if they appear in the
translation column
+             * or by 1 otherwise (scale or shear). We perform this replacement after the
loop searching for
+             * NaN, not inside the loop, in order to not change anything if the search has
been canceled.
+             */
+            for (int k=0; k<indexCount; k += TUPLE_SIZE) {
+                final int i = indexOfNaN[k  ];
+                final int j = indexOfNaN[k+1];
+                final int flatIndex = j*size + i;
+                LU[flatIndex] = (i == lastRowOrColumn) ? 0 : 1;
+                LU[flatIndex + size*size] = 0; // Error term (see 'errorLU') in next method.
+            }
+        }
+        /*
+         * Now apply the inversion.
+         */
+        final MatrixSIS matrix = solve(LU, Y, eltY, size, innerSize);
+        /*
+         * At this point, the matrix has been inverted. If they were any NaN value in the
original
+         * matrix, set the corresponding scale factor and offset to NaN in the resulting
matrix.
+         */
+        for (int k=0; k<indexCount;) {
+            assert (k % TUPLE_SIZE) == 0;
+            final int i = indexOfNaN[k++];
+            final int j = indexOfNaN[k++];
+            final int s = indexOfNaN[k++];
+            if (i != lastRowOrColumn) {
+                // Found a scale factor to set to NaN.
+                matrix.setElement(i, j, Double.NaN); // Note that i,j indices are interchanged.
+                if (matrix.getElement(i, lastRowOrColumn) != 0) {
+                    matrix.setElement(i, lastRowOrColumn, Double.NaN); // = -offset/scale,
so 0 stay 0.
+                }
+            } else if (s >= 0) {
+                // Found a translation factory to set to NaN.
+                matrix.setElement(s, lastRowOrColumn, Double.NaN);
+            }
+        }
+        return matrix;
+    }
 
+    /**
+     * 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;
-        final double[] LU = GeneralMatrix.getExtendedElements(X, size, size, true);
         assert errorLU == GeneralMatrix.indexOfErrors(size, size, LU);
         final int[] pivot = new int[size];
         for (int j=0; j<size; j++) {

Modified: sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/SolverTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/SolverTest.java?rev=1529310&r1=1529309&r2=1529310&view=diff
==============================================================================
--- sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/SolverTest.java
[UTF-8] (original)
+++ sis/branches/JDK7/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/SolverTest.java
[UTF-8] Fri Oct  4 21:32:42 2013
@@ -19,10 +19,14 @@ package org.apache.sis.referencing.opera
 import java.util.Random;
 import Jama.Matrix;
 import org.apache.sis.test.DependsOn;
+import org.apache.sis.test.DependsOnMethod;
 import org.apache.sis.test.TestUtilities;
 import org.apache.sis.test.TestCase;
 import org.junit.Test;
 
+import static java.lang.Double.NaN;
+import static org.junit.Assert.*;
+
 
 /**
  * Tests the {@link Solver} class using <a href="http://math.nist.gov/javanumerics/jama">JAMA</a>
@@ -58,6 +62,19 @@ public final strictfp class SolverTest e
     protected static final double TOLERANCE = 100 * MatrixTestCase.TOLERANCE;
 
     /**
+     * Asserts that the given matrix is equals to the given expected values, up to the {@link
#TOLERANCE}
+     * threshold. This method compares the elements values in two slightly redundant ways.
+     */
+    private static void assertMatrixEquals(final double[] expected,
+            final int numRow, final int numCol, final MatrixSIS actual)
+    {
+        assertEquals("numRow", numRow, actual.getNumRow());
+        assertEquals("numCol", numCol, actual.getNumCol());
+        assertArrayEquals(expected, actual.getElements(), TOLERANCE); // First because more
informative in case of failure.
+        assertTrue(Matrices.create(numRow, numCol, expected).equals(actual, TOLERANCE));
+    }
+
+    /**
      * The matrix to test.
      */
     private MatrixSIS matrix;
@@ -113,4 +130,94 @@ public final strictfp class SolverTest e
             MatrixTestCase.assertMatrixEquals(jama, U, TOLERANCE);
         }
     }
+
+    /**
+     * Tests {@link Solver#inverse(MatrixSIS)} with a square matrix that contains a {@link
Double#NaN} value.
+     *
+     * @throws NoninvertibleMatrixException Should not happen.
+     */
+    @Test
+    @DependsOnMethod("testSolve")
+    public void testInverseWithNaN() throws NoninvertibleMatrixException {
+        /*
+         * Just for making sure that our matrix is correct.
+         */
+        matrix = Matrices.create(5, 5, new double[] {
+            20,  0,   0,   0, -3000,
+            0, -20,   0,   0,  4000,
+            0,   0,   0,   2,    20,
+            0,   0, 400,   0,  2000,
+            0,   0,   0,   0,     1
+        });
+        double[] expected = {
+            0.05,  0,  0,      0,  150,
+            0, -0.05,  0,      0,  200,
+            0,     0,  0, 0.0025,   -5,
+            0,     0,  0.5,    0,  -10,
+            0,     0,  0,      0,    1
+        };
+        MatrixSIS inverse = Solver.inverse(matrix);
+        assertMatrixEquals(expected, 5, 5, inverse);
+        /*
+         * Set a scale factor to NaN. The translation term for the corresponding
+         * dimension become unknown, so it most become NaN in the inverse matrix.
+         */
+        matrix = Matrices.create(5, 5, new double[] {
+            20,  0,   0,   0, -3000,
+            0, -20,   0,   0,  4000,
+            0,   0,   0, NaN,    20,  // Translation is 20: can not be converted.
+            0,   0, 400,   0,  2000,
+            0,   0,   0,   0,     1
+        });
+        expected = new double[] {
+            0.05,  0,  0,      0,  150,
+            0, -0.05,  0,      0,  200,
+            0,     0,  0, 0.0025,   -5,
+            0,     0,  NaN,    0,  NaN,
+            0,     0,  0,      0,    1
+        };
+        inverse = Solver.inverse(matrix);
+        assertMatrixEquals(expected, 5, 5, inverse);
+        /*
+         * Set a scale factor to NaN with translation equals to 0.
+         * The zero value should be preserved, since 0 × any == 0
+         * (ignoring infinities).
+         */
+        matrix = Matrices.create(5, 5, new double[] {
+            20,  0,   0,   0, -3000,
+            0, -20,   0,   0,  4000,
+            0,   0,   0, NaN,     0,  // Translation is 0: should be preserved.
+            0,   0, 400,   0,  2000,
+            0,   0,   0,   0,     1
+        });
+        expected = new double[] {
+            0.05,  0,  0,      0,  150,
+            0, -0.05,  0,      0,  200,
+            0,     0,  0, 0.0025,   -5,
+            0,     0,  NaN,    0,    0,
+            0,     0,  0,      0,    1
+        };
+        inverse = Solver.inverse(matrix);
+        assertMatrixEquals(expected, 5, 5, inverse);
+        /*
+         * Set a translation term to NaN. The translation should be NaN in
+         * the inverse matrix too, but the scale factor can still be compute.
+         */
+        matrix = Matrices.create(5, 5, new double[] {
+            20,  0,   0,   0, -3000,
+            0, -20,   0,   0,  4000,
+            0,   0,   0,   2,   NaN,
+            0,   0, 400,   0,  2000,
+            0,   0,   0,   0,     1
+        });
+        expected = new double[] {
+            0.05,  0,  0,      0,  150,
+            0, -0.05,  0,      0,  200,
+            0,     0,  0, 0.0025,   -5,
+            0,     0,  0.5,    0,  NaN,
+            0,     0,  0,      0,    1
+        };
+        inverse = Solver.inverse(matrix);
+        assertMatrixEquals(expected, 5, 5, inverse);
+    }
 }



Mime
View raw message