JAL-2380 fix Matrix.postMultiply()
authorgmungoc <g.m.carstairs@dundee.ac.uk>
Mon, 16 Jan 2017 09:02:34 +0000 (09:02 +0000)
committergmungoc <g.m.carstairs@dundee.ac.uk>
Mon, 16 Jan 2017 09:02:34 +0000 (09:02 +0000)
src/jalview/math/Matrix.java
test/jalview/math/MatrixTest.java [new file with mode: 0644]

index bb7ed37..28b9d67 100755 (executable)
@@ -57,24 +57,33 @@ public class Matrix
   int maxIter = 45; // fudge - add 15 iterations, just in case
 
   /**
-   * Creates a new Matrix object.
+   * Creates a new Matrix object. For example
    * 
-   * @param value
-   *          DOCUMENT ME!
+   * <pre>
+   *   new Matrix(new double[][] {{2, 3}, {4, 5}, 2, 2)
+   * constructs
+   *   (2 3)
+   *   (4 5)
+   * </pre>
+   * 
+   * Note that ragged arrays (with not all rows, or columns, of the same
+   * length), are not supported by this class. They can be constructed, but
+   * results of operations on them are undefined and may throw exceptions.
+   * 
+   * @param values
+   *          the matrix values in row-major order
    * @param rows
-   *          DOCUMENT ME!
    * @param cols
-   *          DOCUMENT ME!
    */
-  public Matrix(double[][] value, int rows, int cols)
+  public Matrix(double[][] values, int rows, int cols)
   {
     this.rows = rows;
     this.cols = cols;
-    this.value = value;
+    this.value = values;
   }
 
   /**
-   * DOCUMENT ME!
+   * Returns a new matrix which is the transposes of this one
    * 
    * @return DOCUMENT ME!
    */
@@ -113,15 +122,24 @@ public class Matrix
   }
 
   /**
-   * DOCUMENT ME!
+   * Returns a new matrix which is the result of premultiplying this matrix by
+   * the supplied argument. If this of size AxB (A rows and B columns), and the
+   * argument is CxA (C rows and A columns), the result is of size CxB.
    * 
    * @param in
-   *          DOCUMENT ME!
    * 
-   * @return DOCUMENT ME!
+   * @return
+   * @throws IllegalArgumentException
+   *           if the number of columns in the pre-multiplier is not equal to
+   *           the number of rows in the multiplicand (this)
    */
   public Matrix preMultiply(Matrix in)
   {
+    if (in.cols != this.rows)
+    {
+      throw new IllegalArgumentException("Can't pre-multiply " + this.rows
+              + " rows by " + in.cols + " columns");
+    }
     double[][] tmp = new double[in.rows][this.cols];
 
     for (int i = 0; i < in.rows; i++)
@@ -141,12 +159,10 @@ public class Matrix
   }
 
   /**
-   * DOCUMENT ME!
    * 
    * @param in
-   *          DOCUMENT ME!
    * 
-   * @return DOCUMENT ME!
+   * @return
    */
   public double[] vectorPostMultiply(double[] in)
   {
@@ -166,37 +182,34 @@ public class Matrix
   }
 
   /**
-   * DOCUMENT ME!
+   * Returns a new matrix which is the result of postmultiplying this matrix by
+   * the supplied argument. If this of size AxB (A rows and B columns), and the
+   * argument is BxC (B rows and C columns), the result is of size AxC.
+   * <p>
+   * This method simply returns the result of in.preMultiply(this)
    * 
    * @param in
-   *          DOCUMENT ME!
    * 
-   * @return DOCUMENT ME!
+   * @return
+   * @throws IllegalArgumentException
+   *           if the number of rows in the post-multiplier is not equal to the
+   *           number of columns in the multiplicand (this)
+   * @see #preMultiply(Matrix)
    */
   public Matrix postMultiply(Matrix in)
   {
-    double[][] out = new double[this.rows][in.cols];
-
-    for (int i = 0; i < this.rows; i++)
+    if (in.rows != this.cols)
     {
-      for (int j = 0; j < in.cols; j++)
-      {
-        out[i][j] = 0.0;
-
-        for (int k = 0; k < rows; k++)
-        {
-          out[i][j] = out[i][j] + (value[i][k] * in.value[k][j]);
-        }
-      }
+      throw new IllegalArgumentException("Can't post-multiply " + this.cols
+              + " columns by " + in.rows + " rows");
     }
-
-    return new Matrix(out, this.cols, in.rows);
+    return in.preMultiply(this);
   }
 
   /**
-   * DOCUMENT ME!
+   * Answers a new matrix with a copy of the values in this one
    * 
-   * @return DOCUMENT ME!
+   * @return
    */
   public Matrix copy()
   {
@@ -204,10 +217,11 @@ public class Matrix
 
     for (int i = 0; i < rows; i++)
     {
-      for (int j = 0; j < cols; j++)
-      {
-        newmat[i][j] = value[i][j];
-      }
+      System.arraycopy(value[i], 0, newmat[i], 0, value[i].length);
+      // for (int j = 0; j < cols; j++)
+      // {
+      // newmat[i][j] = value[i][j];
+      // }
     }
 
     return new Matrix(newmat, rows, cols);
@@ -738,20 +752,19 @@ public class Matrix
   }
 
   /**
-   * DOCUMENT ME!
+   * Returns an array containing the values in the specified column
    * 
-   * @param n
-   *          DOCUMENT ME!
+   * @param col
    * 
-   * @return DOCUMENT ME!
+   * @return
    */
-  public double[] getColumn(int n)
+  public double[] getColumn(int col)
   {
     double[] out = new double[rows];
 
     for (int i = 0; i < rows; i++)
     {
-      out[i] = value[i][n];
+      out[i] = value[i][col];
     }
 
     return out;
@@ -784,85 +797,4 @@ public class Matrix
       Format.print(ps, "%15.4e", e[j]);
     }
   }
-
-  /**
-   * DOCUMENT ME!
-   * 
-   * @param args
-   *          DOCUMENT ME!
-   */
-  public static void main(String[] args) throws Exception
-  {
-    int n = Integer.parseInt(args[0]);
-    double[][] in = new double[n][n];
-
-    for (int i = 0; i < n; i++)
-    {
-      for (int j = 0; j < n; j++)
-      {
-        in[i][j] = (double) Math.random();
-      }
-    }
-
-    Matrix origmat = new Matrix(in, n, n);
-
-    // System.out.println(" --- Original matrix ---- ");
-    // / origmat.print(System.out);
-    // System.out.println();
-    // System.out.println(" --- transpose matrix ---- ");
-    Matrix trans = origmat.transpose();
-
-    // trans.print(System.out);
-    // System.out.println();
-    // System.out.println(" --- OrigT * Orig ---- ");
-    Matrix symm = trans.postMultiply(origmat);
-
-    // symm.print(System.out);
-    // System.out.println();
-    // Copy the symmetric matrix for later
-    // Matrix origsymm = symm.copy();
-
-    // This produces the tridiagonal transformation matrix
-    // long tstart = System.currentTimeMillis();
-    symm.tred();
-
-    // long tend = System.currentTimeMillis();
-
-    // System.out.println("Time take for tred = " + (tend-tstart) + "ms");
-    // System.out.println(" ---Tridiag transform matrix ---");
-    // symm.print(System.out);
-    // System.out.println();
-    // System.out.println(" --- D vector ---");
-    // symm.printD(System.out);
-    // System.out.println();
-    // System.out.println(" --- E vector ---");
-    // symm.printE(System.out);
-    // System.out.println();
-    // Now produce the diagonalization matrix
-    // tstart = System.currentTimeMillis();
-    symm.tqli();
-    // tend = System.currentTimeMillis();
-
-    // System.out.println("Time take for tqli = " + (tend-tstart) + " ms");
-    // System.out.println(" --- New diagonalization matrix ---");
-    // symm.print(System.out);
-    // System.out.println();
-    // System.out.println(" --- D vector ---");
-    // symm.printD(System.out);
-    // System.out.println();
-    // System.out.println(" --- E vector ---");
-    // symm.printE(System.out);
-    // System.out.println();
-    // System.out.println(" --- First eigenvector --- ");
-    // double[] eigenv = symm.getColumn(0);
-    // for (int i=0; i < eigenv.length;i++) {
-    // Format.print(System.out,"%15.4f",eigenv[i]);
-    // }
-    // System.out.println();
-    // double[] neigenv = origsymm.vectorPostMultiply(eigenv);
-    // for (int i=0; i < neigenv.length;i++) {
-    // Format.print(System.out,"%15.4f",neigenv[i]/symm.d[0]);
-    // }
-    // System.out.println();
-  }
 }
diff --git a/test/jalview/math/MatrixTest.java b/test/jalview/math/MatrixTest.java
new file mode 100644 (file)
index 0000000..795b2fa
--- /dev/null
@@ -0,0 +1,239 @@
+package jalview.math;
+
+import static org.testng.Assert.assertEquals;
+import static org.testng.Assert.fail;
+
+import java.util.Arrays;
+
+import org.testng.annotations.Test;
+
+public class MatrixTest
+{
+  @Test(groups = "Timing")
+  public void testPreMultiply_timing()
+  {
+    int rows = 500;
+    int cols = 1000;
+    double[][] d1 = new double[rows][cols];
+    double[][] d2 = new double[cols][rows];
+    Matrix m1 = new Matrix(d1, rows, cols);
+    Matrix m2 = new Matrix(d2, cols, rows);
+    long start = System.currentTimeMillis();
+    m1.preMultiply(m2);
+    long elapsed = System.currentTimeMillis() - start;
+    System.out.println(rows + "x" + cols
+            + " multiplications of double took " + elapsed + "ms");
+  }
+
+  @Test(groups = "Functional")
+  public void testPreMultiply()
+  {
+    Matrix m1 = new Matrix(new double[][] { { 2, 3, 4 } }, 1, 3); // 1x3
+    Matrix m2 = new Matrix(new double[][] { { 5 }, { 6 }, { 7 } }, 3, 1); // 3x1
+
+    /*
+     * 1x3 times 3x1 is 1x1
+     * 2x5 + 3x6 + 4*7 =  56
+     */
+    Matrix m3 = m2.preMultiply(m1);
+    assertEquals(m3.rows, 1);
+    assertEquals(m3.cols, 1);
+    assertEquals(m3.value[0][0], 56d);
+
+    /*
+     * 3x1 times 1x3 is 3x3
+     */
+    m3 = m1.preMultiply(m2);
+    assertEquals(m3.rows, 3);
+    assertEquals(m3.cols, 3);
+    assertEquals(Arrays.toString(m3.value[0]), "[10.0, 15.0, 20.0]");
+    assertEquals(Arrays.toString(m3.value[1]), "[12.0, 18.0, 24.0]");
+    assertEquals(Arrays.toString(m3.value[2]), "[14.0, 21.0, 28.0]");
+  }
+
+  @Test(
+    groups = "Functional",
+    expectedExceptions = { IllegalArgumentException.class })
+  public void testPreMultiply_tooManyColumns()
+  {
+    Matrix m1 = new Matrix(new double[][] { { 2, 3, 4 }, { 3, 4, 5 } }, 2,
+            3); // 2x3
+
+    /*
+     * 2x3 times 2x3 invalid operation - 
+     * multiplier has more columns than multiplicand has rows
+     */
+    m1.preMultiply(m1);
+    fail("Expected exception");
+  }
+
+  @Test(
+    groups = "Functional",
+    expectedExceptions = { IllegalArgumentException.class })
+  public void testPreMultiply_tooFewColumns()
+  {
+    Matrix m1 = new Matrix(new double[][] { { 2, 3, 4 }, { 3, 4, 5 } }, 2,
+            3); // 2x3
+
+    /*
+     * 3x2 times 3x2 invalid operation - 
+     * multiplier has more columns than multiplicand has row
+     */
+    m1.preMultiply(m1);
+    fail("Expected exception");
+  }
+  
+  
+  private boolean matrixEquals(Matrix m1, Matrix m2) {
+    return Arrays.deepEquals(m1.value, m2.value);
+  }
+
+  @Test(groups = "Functional")
+  public void testPostMultiply()
+  {
+    /*
+     * Square matrices
+     * (2 3) . (10   100)
+     * (4 5)   (1000 10000)
+     * =
+     * (3020 30200)
+     * (5040 50400)
+     */
+    Matrix m1 = new Matrix(new double[][] { { 2, 3 }, { 4, 5 } }, 2, 2);
+    Matrix m2 = new Matrix(new double[][] { { 10, 100 }, { 1000, 10000 } },
+            2, 2);
+    Matrix m3 = m1.postMultiply(m2);
+    assertEquals(Arrays.toString(m3.value[0]), "[3020.0, 30200.0]");
+    assertEquals(Arrays.toString(m3.value[1]), "[5040.0, 50400.0]");
+
+    /*
+     * also check m2.preMultiply(m1) - should be same as m1.postMultiply(m2) 
+     */
+    m3 = m2.preMultiply(m1);
+    assertEquals(Arrays.toString(m3.value[0]), "[3020.0, 30200.0]");
+    assertEquals(Arrays.toString(m3.value[1]), "[5040.0, 50400.0]");
+
+    /*
+     * m1 has more rows than columns
+     * (2).(10 100 1000) = (20 200 2000)
+     * (3)                 (30 300 3000)
+     */
+    m1 = new Matrix(new double[][] { { 2 }, { 3 } }, 2, 1);
+    m2 = new Matrix(new double[][] { { 10, 100, 1000 } }, 1, 3);
+    m3 = m1.postMultiply(m2);
+    assertEquals(m3.rows, 2);
+    assertEquals(m3.cols, 3);
+    assertEquals(Arrays.toString(m3.value[0]), "[20.0, 200.0, 2000.0]");
+    assertEquals(Arrays.toString(m3.value[1]), "[30.0, 300.0, 3000.0]");
+    m3 = m2.preMultiply(m1);
+    assertEquals(m3.rows, 2);
+    assertEquals(m3.cols, 3);
+    assertEquals(Arrays.toString(m3.value[0]), "[20.0, 200.0, 2000.0]");
+    assertEquals(Arrays.toString(m3.value[1]), "[30.0, 300.0, 3000.0]");
+
+    /*
+     * m1 has more columns than rows
+     * (2 3 4) . (5 4) = (56 25)
+     *           (6 3) 
+     *           (7 2)
+     * [0, 0] = 2*5 + 3*6 + 4*7 = 56
+     * [0, 1] = 2*4 + 3*3 + 4*2 = 25  
+     */
+    m1 = new Matrix(new double[][] { { 2, 3, 4 } }, 1, 3);
+    m2 = new Matrix(new double[][] { { 5, 4 }, { 6, 3 }, { 7, 2 } }, 3, 2);
+    m3 = m1.postMultiply(m2);
+    assertEquals(m3.rows, 1);
+    assertEquals(m3.cols, 2);
+    assertEquals(m3.value[0][0], 56d);
+    assertEquals(m3.value[0][1], 25d);
+
+    /*
+     * and check premultiply equivalent
+     */
+    m3 = m2.preMultiply(m1);
+    assertEquals(m3.rows, 1);
+    assertEquals(m3.cols, 2);
+    assertEquals(m3.value[0][0], 56d);
+    assertEquals(m3.value[0][1], 25d);
+  }
+
+  /**
+   * main method extracted from Matrix
+   * 
+   * @param args
+   */
+  public static void main(String[] args) throws Exception
+  {
+    int n = Integer.parseInt(args[0]);
+    double[][] in = new double[n][n];
+  
+    for (int i = 0; i < n; i++)
+    {
+      for (int j = 0; j < n; j++)
+      {
+        in[i][j] = Math.random();
+      }
+    }
+  
+    Matrix origmat = new Matrix(in, n, n);
+  
+    // System.out.println(" --- Original matrix ---- ");
+    // / origmat.print(System.out);
+    // System.out.println();
+    // System.out.println(" --- transpose matrix ---- ");
+    Matrix trans = origmat.transpose();
+  
+    // trans.print(System.out);
+    // System.out.println();
+    // System.out.println(" --- OrigT * Orig ---- ");
+    Matrix symm = trans.postMultiply(origmat);
+  
+    // symm.print(System.out);
+    // System.out.println();
+    // Copy the symmetric matrix for later
+    // Matrix origsymm = symm.copy();
+  
+    // This produces the tridiagonal transformation matrix
+    // long tstart = System.currentTimeMillis();
+    symm.tred();
+  
+    // long tend = System.currentTimeMillis();
+  
+    // System.out.println("Time take for tred = " + (tend-tstart) + "ms");
+    // System.out.println(" ---Tridiag transform matrix ---");
+    // symm.print(System.out);
+    // System.out.println();
+    // System.out.println(" --- D vector ---");
+    // symm.printD(System.out);
+    // System.out.println();
+    // System.out.println(" --- E vector ---");
+    // symm.printE(System.out);
+    // System.out.println();
+    // Now produce the diagonalization matrix
+    // tstart = System.currentTimeMillis();
+    symm.tqli();
+    // tend = System.currentTimeMillis();
+  
+    // System.out.println("Time take for tqli = " + (tend-tstart) + " ms");
+    // System.out.println(" --- New diagonalization matrix ---");
+    // symm.print(System.out);
+    // System.out.println();
+    // System.out.println(" --- D vector ---");
+    // symm.printD(System.out);
+    // System.out.println();
+    // System.out.println(" --- E vector ---");
+    // symm.printE(System.out);
+    // System.out.println();
+    // System.out.println(" --- First eigenvector --- ");
+    // double[] eigenv = symm.getColumn(0);
+    // for (int i=0; i < eigenv.length;i++) {
+    // Format.print(System.out,"%15.4f",eigenv[i]);
+    // }
+    // System.out.println();
+    // double[] neigenv = origsymm.vectorPostMultiply(eigenv);
+    // for (int i=0; i < neigenv.length;i++) {
+    // Format.print(System.out,"%15.4f",neigenv[i]/symm.d[0]);
+    // }
+    // System.out.println();
+  }
+}