From e06941005a462eff3578ec863630edbc519b560f Mon Sep 17 00:00:00 2001 From: gmungoc Date: Mon, 16 Jan 2017 09:02:34 +0000 Subject: [PATCH] JAL-2380 fix Matrix.postMultiply() --- src/jalview/math/Matrix.java | 180 +++++++++------------------- test/jalview/math/MatrixTest.java | 239 +++++++++++++++++++++++++++++++++++++ 2 files changed, 295 insertions(+), 124 deletions(-) create mode 100644 test/jalview/math/MatrixTest.java diff --git a/src/jalview/math/Matrix.java b/src/jalview/math/Matrix.java index bb7ed37..28b9d67 100755 --- a/src/jalview/math/Matrix.java +++ b/src/jalview/math/Matrix.java @@ -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! + * 
+   *   new Matrix(new double[][] {{2, 3}, {4, 5}, 2, 2)
+   * constructs
+   *   (2 3)
+   *   (4 5)
+   *
+ * + * 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. + *

+ * 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 index 0000000..795b2fa --- /dev/null +++ b/test/jalview/math/MatrixTest.java @@ -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(); + } +} -- 1.7.10.2