JAL-2379 made SparseMatrix.divide() so as to exactly match
authorgmungoc <g.m.carstairs@dundee.ac.uk>
Mon, 23 Jan 2017 20:25:42 +0000 (20:25 +0000)
committergmungoc <g.m.carstairs@dundee.ac.uk>
Mon, 23 Jan 2017 20:25:42 +0000 (20:25 +0000)
Matrix.divide()

src/jalview/ext/android/SparseDoubleArray.java
src/jalview/math/Matrix.java
src/jalview/math/SparseMatrix.java
test/jalview/ext/android/SparseDoubleArrayTest.java
test/jalview/math/SparseMatrixTest.java

index 5ae5afe..eaf059c 100644 (file)
@@ -425,13 +425,17 @@ public class SparseDoubleArray implements Cloneable
    * @oparam toAdd
    * @return the new value for the key
    */
-  public double multiply(int key, double d)
+  public double divide(int key, double divisor)
   {
     double newValue = 0d;
+    if (divisor == 0d)
+    {
+      return newValue;
+    }
     int i = ContainerHelpers.binarySearch(mKeys, mSize, key);
     if (i >= 0)
     {
-      mValues[i] *= d;
+      mValues[i] /= divisor;
       newValue = mValues[i];
     }
     return newValue;
index 0f93eb5..2b2314b 100755 (executable)
@@ -282,20 +282,16 @@ public class Matrix implements MatrixI
 
         if (scale == 0.0)
         {
-          // e[i - 1] = value[i - 1][l - 1];
           e[i - 1] = getValue(i - 1, l - 1);
         }
         else
         {
           for (k = 1; k <= l; k++)
           {
-            // value[i - 1][k - 1] /= scale;
-            // h += (value[i - 1][k - 1] * value[i - 1][k - 1]);
             double v = divideValue(i - 1, k - 1, scale);
             h += v * v;
           }
 
-          // f = value[i - 1][l - 1];
           f = getValue(i - 1, l - 1);
 
           if (f > 0)
@@ -309,34 +305,26 @@ public class Matrix implements MatrixI
 
           e[i - 1] = scale * g;
           h -= (f * g);
-          // value[i - 1][l - 1] = f - g;
           setValue(i - 1, l - 1, f - g);
-          // System.out.println(String.format("%d %d %f %f %f %f %f %.5e", i,
-          // l,
-          // scale, f, g, h, getValue(i - 1, l - 1), checksum()));
           f = 0.0;
 
           for (j = 1; j <= l; j++)
           {
-            // value[j - 1][i - 1] = value[i - 1][j - 1] / h;
             double val = getValue(i - 1, j - 1) / h;
             setValue(j - 1, i - 1, val);
             g = 0.0;
 
             for (k = 1; k <= j; k++)
             {
-              // g += (value[j - 1][k - 1] * value[i - 1][k - 1]);
               g += (getValue(j - 1, k - 1) * getValue(i - 1, k - 1));
             }
 
             for (k = j + 1; k <= l; k++)
             {
-              // g += (value[k - 1][j - 1] * value[i - 1][k - 1]);
               g += (getValue(k - 1, j - 1) * getValue(i - 1, k - 1));
             }
 
             e[j - 1] = g / h;
-            // f += (e[j - 1] * value[i - 1][j - 1]);
             f += (e[j - 1] * getValue(i - 1, j - 1));
           }
 
@@ -344,15 +332,12 @@ public class Matrix implements MatrixI
 
           for (j = 1; j <= l; j++)
           {
-            // f = value[i - 1][j - 1];
             f = getValue(i - 1, j - 1);
             g = e[j - 1] - (hh * f);
             e[j - 1] = g;
 
             for (k = 1; k <= j; k++)
             {
-              // value[j - 1][k - 1] -= ((f * e[k - 1]) + (g * value[i - 1][k -
-              // 1]));
               double val = (f * e[k - 1]) + (g * getValue(i - 1, k - 1));
               addValue(j - 1, k - 1, -val);
             }
@@ -361,7 +346,6 @@ public class Matrix implements MatrixI
       }
       else
       {
-        // e[i - 1] = value[i - 1][l - 1];
         e[i - 1] = getValue(i - 1, l - 1);
       }
 
@@ -383,27 +367,21 @@ public class Matrix implements MatrixI
 
           for (k = 1; k <= l; k++)
           {
-            // g += (value[i - 1][k - 1] * value[k - 1][j - 1]);
             g += (getValue(i - 1, k - 1) * getValue(k - 1, j - 1));
           }
 
           for (k = 1; k <= l; k++)
           {
-            // value[k - 1][j - 1] -= (g * value[k - 1][i - 1]);
-            addValue(k - 1, j - 1, -(g * getValue(k - 1, i - 1)));
+            double x = addValue(k - 1, j - 1, -(g * getValue(k - 1, i - 1)));
           }
         }
       }
 
-      // d[i - 1] = value[i - 1][i - 1];
-      // value[i - 1][i - 1] = 1.0;
       d[i - 1] = getValue(i - 1, i - 1);
       setValue(i - 1, i - 1, 1.0);
 
       for (j = 1; j <= l; j++)
       {
-        // value[j - 1][i - 1] = 0.0;
-        // value[i - 1][j - 1] = 0.0;
         setValue(j - 1, i - 1, 0.0);
         setValue(i - 1, j - 1, 0.0);
       }
index 4a9d427..72f0963 100644 (file)
@@ -189,7 +189,7 @@ public class SparseMatrix extends Matrix
     {
       return getValue(i, j);
     }
-    double v = sparseColumns[j].multiply(i, 1 / divisor);
+    double v = sparseColumns[j].divide(i, divisor);
     return v;
   }
 
index 58c36e6..7d64a28 100644 (file)
@@ -36,4 +36,18 @@ public class SparseDoubleArrayTest
     s.add(12, 9.8d);
     assertEquals(s.get(12), 9.8d);
   }
+
+  @Test
+  public void testDivide()
+  {
+    double delta = 1.0e-10;
+    double[] d = new double[] { 0d, 2.4d, 1.2d, 0d, -4.8d, -3.6d };
+    SparseDoubleArray s = new SparseDoubleArray(d);
+    assertEquals(s.divide(0, 1d), 0d); // no such entry
+    assertEquals(s.divide(2, 0d), 0d); // zero divisor
+    assertEquals(s.divide(1, 2d), 1.2d, delta); // + / +
+    assertEquals(s.divide(2, -2d), -0.6d, delta); // + / -
+    assertEquals(s.divide(4, 3d), -1.6d, delta); // - / +
+    assertEquals(s.divide(5, -3d), 1.2d, delta); // - / -
+  }
 }
index 5e32e18..3c2ccaa 100644 (file)
@@ -134,7 +134,7 @@ public class SparseMatrixTest
      * also check m2.preMultiply(m1) - should be same as m1.postMultiply(m2) 
      */
     MatrixI m4 = m2.preMultiply(m1);
-    assertMatricesMatch(m3, m4);
+    assertMatricesMatch(m3, m4, 0.00001d);
 
     /*
      * m1 has more rows than columns
@@ -154,7 +154,7 @@ public class SparseMatrixTest
     assertEquals(m3.getValue(1, 2), 3000d);
 
     m4 = m2.preMultiply(m1);
-    assertMatricesMatch(m3, m4);
+    assertMatricesMatch(m3, m4, 0.00001d);
 
     /*
      * m1 has more columns than rows
@@ -176,7 +176,7 @@ public class SparseMatrixTest
      * and check premultiply equivalent
      */
     m4 = m2.preMultiply(m1);
-    assertMatricesMatch(m3, m4);
+    assertMatricesMatch(m3, m4, 0.00001d);
   }
 
   @Test(groups = "Timing")
@@ -201,8 +201,7 @@ public class SparseMatrixTest
     /*
      * make a pseudo-random symmetric matrix as required for tred/tqli
      */
-    // TODO why does this fail for rows > 9??
-    int rows = 9;
+    int rows = 10;
     int cols = rows;
     double[][] d = getSparseValues(rows, cols, 3);
 
@@ -220,13 +219,13 @@ public class SparseMatrixTest
     }
     Matrix m1 = new Matrix(d);
     Matrix m2 = new SparseMatrix(d1);
-    assertMatricesMatch(m1, m2); // sanity check
+    assertMatricesMatch(m1, m2, 0.00001d); // sanity check
     m1.tred();
     m2.tred();
-    assertMatricesMatch(m1, m2);
+    assertMatricesMatch(m1, m2, 0.00001d);
   }
 
-  private void assertMatricesMatch(MatrixI m1, MatrixI m2)
+  private void assertMatricesMatch(MatrixI m1, MatrixI m2, double delta)
   {
     if (m1.height() != m2.height())
     {
@@ -248,7 +247,7 @@ public class SparseMatrixTest
         }
       }
     }
-    ArrayAsserts.assertArrayEquals(m1.getD(), m2.getD(), 0.00001d);
+    ArrayAsserts.assertArrayEquals(m1.getD(), m2.getD(), delta);
     ArrayAsserts.assertArrayEquals(m1.getE(), m2.getE(), 0.00001d);
   }
 
@@ -303,11 +302,11 @@ public class SparseMatrixTest
     // have to do tred() before doing tqli()
     m1.tred();
     m2.tred();
-    assertMatricesMatch(m1, m2);
+    assertMatricesMatch(m1, m2, 0.00001d);
 
     m1.tqli();
     m2.tqli();
-    assertMatricesMatch(m1, m2);
+    assertMatricesMatch(m1, m2, 0.00001d);
   }
 
   /**
@@ -409,9 +408,9 @@ public class SparseMatrixTest
     }
     Matrix m1 = new SparseMatrix(d);
     Matrix m2 = new SparseMatrix(d1);
-    assertMatricesMatch(m1, m2); // sanity check
+    assertMatricesMatch(m1, m2, 1.0e16); // sanity check
     m1.tred();
     m2.tred();
-    assertMatricesMatch(m1, m2);
+    assertMatricesMatch(m1, m2, 0.00001d);
   }
 }
\ No newline at end of file