test/jalview/math/MatrixTest.java

   1 package jalview.math;
   2
   3 import static org.testng.Assert.assertEquals;
   4 import static org.testng.Assert.fail;
   5
   6 import java.util.Arrays;
   7 import java.util.Random;
   8
   9 import org.testng.annotations.Test;
  10 import org.testng.internal.junit.ArrayAsserts;
  11
  12 public class MatrixTest
  13 {
  14   final static double DELTA = 0.0001d;
  15
  16   @Test(groups = "Timing")
  17   public void testPreMultiply_timing()
  18   {
  19     int rows = 500;
  20     int cols = 1000;
  21     double[][] d1 = new double[rows][cols];
  22     double[][] d2 = new double[cols][rows];
  23     Matrix m1 = new Matrix(d1);
  24     Matrix m2 = new Matrix(d2);
  25     long start = System.currentTimeMillis();
  26     m1.preMultiply(m2);
  27     long elapsed = System.currentTimeMillis() - start;
  28     System.out.println(rows + "x" + cols
  29             + " multiplications of double took " + elapsed + "ms");
  30   }
  31
  32   @Test(groups = "Functional")
  33   public void testPreMultiply()
  34   {
  35     MatrixI m1 = new Matrix(new double[][] { { 2, 3, 4 } }); // 1x3
  36     MatrixI m2 = new Matrix(new double[][] { { 5 }, { 6 }, { 7 } }); // 3x1
  37
  38     /*
  39      * 1x3 times 3x1 is 1x1
  40      * 2x5 + 3x6 + 4*7 =  56
  41      */
  42     MatrixI m3 = m2.preMultiply(m1);
  43     assertEquals(m3.height(), 1);
  44     assertEquals(m3.width(), 1);
  45     assertEquals(m3.getValue(0, 0), 56d);
  46
  47     /*
  48      * 3x1 times 1x3 is 3x3
  49      */
  50     m3 = m1.preMultiply(m2);
  51     assertEquals(m3.height(), 3);
  52     assertEquals(m3.width(), 3);
  53     assertEquals(Arrays.toString(m3.getRow(0)), "[10.0, 15.0, 20.0]");
  54     assertEquals(Arrays.toString(m3.getRow(1)), "[12.0, 18.0, 24.0]");
  55     assertEquals(Arrays.toString(m3.getRow(2)), "[14.0, 21.0, 28.0]");
  56   }
  57
  58   @Test(
  59     groups = "Functional",
  60     expectedExceptions = { IllegalArgumentException.class })
  61   public void testPreMultiply_tooManyColumns()
  62   {
  63     Matrix m1 = new Matrix(new double[][] { { 2, 3, 4 }, { 3, 4, 5 } }); // 2x3
  64
  65     /*
  66      * 2x3 times 2x3 invalid operation -
  67      * multiplier has more columns than multiplicand has rows
  68      */
  69     m1.preMultiply(m1);
  70     fail("Expected exception");
  71   }
  72
  73   @Test(
  74     groups = "Functional",
  75     expectedExceptions = { IllegalArgumentException.class })
  76   public void testPreMultiply_tooFewColumns()
  77   {
  78     Matrix m1 = new Matrix(new double[][] { { 2, 3, 4 }, { 3, 4, 5 } }); // 2x3
  79
  80     /*
  81      * 3x2 times 3x2 invalid operation -
  82      * multiplier has more columns than multiplicand has row
  83      */
  84     m1.preMultiply(m1);
  85     fail("Expected exception");
  86   }
  87
  88
  89   private boolean matrixEquals(Matrix m1, Matrix m2) {
  90     if (m1.width() != m2.width() || m1.height() != m2.height())
  91     {
  92       return false;
  93     }
  94     for (int i = 0; i < m1.height(); i++)
  95     {
  96       if (!Arrays.equals(m1.getRow(i), m2.getRow(i)))
  97       {
  98         return false;
  99       }
 100     }
 101     return true;
 102   }
 103
 104   @Test(groups = "Functional")
 105   public void testPostMultiply()
 106   {
 107     /*
 108      * Square matrices
 109      * (2 3) . (10   100)
 110      * (4 5)   (1000 10000)
 111      * =
 112      * (3020 30200)
 113      * (5040 50400)
 114      */
 115     MatrixI m1 = new Matrix(new double[][] { { 2, 3 }, { 4, 5 } });
 116     MatrixI m2 = new Matrix(new double[][] { { 10, 100 }, { 1000, 10000 } });
 117     MatrixI m3 = m1.postMultiply(m2);
 118     assertEquals(Arrays.toString(m3.getRow(0)), "[3020.0, 30200.0]");
 119     assertEquals(Arrays.toString(m3.getRow(1)), "[5040.0, 50400.0]");
 120
 121     /*
 122      * also check m2.preMultiply(m1) - should be same as m1.postMultiply(m2)
 123      */
 124     m3 = m2.preMultiply(m1);
 125     assertEquals(Arrays.toString(m3.getRow(0)), "[3020.0, 30200.0]");
 126     assertEquals(Arrays.toString(m3.getRow(1)), "[5040.0, 50400.0]");
 127
 128     /*
 129      * m1 has more rows than columns
 130      * (2).(10 100 1000) = (20 200 2000)
 131      * (3)                 (30 300 3000)
 132      */
 133     m1 = new Matrix(new double[][] { { 2 }, { 3 } });
 134     m2 = new Matrix(new double[][] { { 10, 100, 1000 } });
 135     m3 = m1.postMultiply(m2);
 136     assertEquals(m3.height(), 2);
 137     assertEquals(m3.width(), 3);
 138     assertEquals(Arrays.toString(m3.getRow(0)), "[20.0, 200.0, 2000.0]");
 139     assertEquals(Arrays.toString(m3.getRow(1)), "[30.0, 300.0, 3000.0]");
 140     m3 = m2.preMultiply(m1);
 141     assertEquals(m3.height(), 2);
 142     assertEquals(m3.width(), 3);
 143     assertEquals(Arrays.toString(m3.getRow(0)), "[20.0, 200.0, 2000.0]");
 144     assertEquals(Arrays.toString(m3.getRow(1)), "[30.0, 300.0, 3000.0]");
 145
 146     /*
 147      * m1 has more columns than rows
 148      * (2 3 4) . (5 4) = (56 25)
 149      *           (6 3)
 150      *           (7 2)
 151      * [0, 0] = 2*5 + 3*6 + 4*7 = 56
 152      * [0, 1] = 2*4 + 3*3 + 4*2 = 25
 153      */
 154     m1 = new Matrix(new double[][] { { 2, 3, 4 } });
 155     m2 = new Matrix(new double[][] { { 5, 4 }, { 6, 3 }, { 7, 2 } });
 156     m3 = m1.postMultiply(m2);
 157     assertEquals(m3.height(), 1);
 158     assertEquals(m3.width(), 2);
 159     assertEquals(m3.getRow(0)[0], 56d);
 160     assertEquals(m3.getRow(0)[1], 25d);
 161
 162     /*
 163      * and check premultiply equivalent
 164      */
 165     m3 = m2.preMultiply(m1);
 166     assertEquals(m3.height(), 1);
 167     assertEquals(m3.width(), 2);
 168     assertEquals(m3.getRow(0)[0], 56d);
 169     assertEquals(m3.getRow(0)[1], 25d);
 170   }
 171
 172   /**
 173    * main method extracted from Matrix
 174    *
 175    * @param args
 176    */
 177   public static void main(String[] args) throws Exception
 178   {
 179     int n = Integer.parseInt(args[0]);
 180     double[][] in = new double[n][n];
 181
 182     for (int i = 0; i < n; i++)
 183     {
 184       for (int j = 0; j < n; j++)
 185       {
 186         in[i][j] = Math.random();
 187       }
 188     }
 189
 190     Matrix origmat = new Matrix(in);
 191
 192     // System.out.println(" --- Original matrix ---- ");
 193     // / origmat.print(System.out);
 194     // System.out.println();
 195     // System.out.println(" --- transpose matrix ---- ");
 196     MatrixI trans = origmat.transpose();
 197
 198     // trans.print(System.out);
 199     // System.out.println();
 200     // System.out.println(" --- OrigT * Orig ---- ");
 201     MatrixI symm = trans.postMultiply(origmat);
 202
 203     // symm.print(System.out);
 204     // System.out.println();
 205     // Copy the symmetric matrix for later
 206     // Matrix origsymm = symm.copy();
 207
 208     // This produces the tridiagonal transformation matrix
 209     // long tstart = System.currentTimeMillis();
 210     symm.tred();
 211
 212     // long tend = System.currentTimeMillis();
 213
 214     // System.out.println("Time take for tred = " + (tend-tstart) + "ms");
 215     // System.out.println(" ---Tridiag transform matrix ---");
 216     // symm.print(System.out);
 217     // System.out.println();
 218     // System.out.println(" --- D vector ---");
 219     // symm.printD(System.out);
 220     // System.out.println();
 221     // System.out.println(" --- E vector ---");
 222     // symm.printE(System.out);
 223     // System.out.println();
 224     // Now produce the diagonalization matrix
 225     // tstart = System.currentTimeMillis();
 226     symm.tqli();
 227     // tend = System.currentTimeMillis();
 228
 229     // System.out.println("Time take for tqli = " + (tend-tstart) + " ms");
 230     // System.out.println(" --- New diagonalization matrix ---");
 231     // symm.print(System.out);
 232     // System.out.println();
 233     // System.out.println(" --- D vector ---");
 234     // symm.printD(System.out);
 235     // System.out.println();
 236     // System.out.println(" --- E vector ---");
 237     // symm.printE(System.out);
 238     // System.out.println();
 239     // System.out.println(" --- First eigenvector --- ");
 240     // double[] eigenv = symm.getColumn(0);
 241     // for (int i=0; i < eigenv.length;i++) {
 242     // Format.print(System.out,"%15.4f",eigenv[i]);
 243     // }
 244     // System.out.println();
 245     // double[] neigenv = origsymm.vectorPostMultiply(eigenv);
 246     // for (int i=0; i < neigenv.length;i++) {
 247     // Format.print(System.out,"%15.4f",neigenv[i]/symm.d[0]);
 248     // }
 249     // System.out.println();
 250   }
 251
 252   @Test(groups = "Timing")
 253   public void testSign()
 254   {
 255     assertEquals(Matrix.sign(-1, -2), -1d);
 256     assertEquals(Matrix.sign(-1, 2), 1d);
 257     assertEquals(Matrix.sign(-1, 0), 1d);
 258     assertEquals(Matrix.sign(1, -2), -1d);
 259     assertEquals(Matrix.sign(1, 2), 1d);
 260     assertEquals(Matrix.sign(1, 0), 1d);
 261   }
 262
 263   /**
 264    * Helper method to make values for a sparse, pseudo-random symmetric matrix
 265    *
 266    * @param rows
 267    * @param cols
 268    * @param occupancy
 269    *          one in 'occupancy' entries will be non-zero
 270    * @return
 271    */
 272   public double[][] getSparseValues(int rows, int cols, int occupancy)
 273   {
 274     Random r = new Random(1729);
 275
 276     /*
 277      * generate whole number values between -12 and +12
 278      * (to mimic score matrices used in Jalview)
 279      */
 280     double[][] d = new double[rows][cols];
 281     int m = 0;
 282     for (int i = 0; i < rows; i++)
 283     {
 284       if (++m % occupancy == 0)
 285       {
 286         d[i][i] = r.nextInt() % 13; // diagonal
 287       }
 288       for (int j = 0; j < i; j++)
 289       {
 290         if (++m % occupancy == 0)
 291         {
 292           d[i][j] = r.nextInt() % 13;
 293           d[j][i] = d[i][j];
 294         }
 295       }
 296     }
 297     return d;
 298
 299   }
 300
 301   /**
 302    * Verify that the results of method tred() are the same if the calculation is
 303    * redone
 304    */
 305   @Test(groups = "Functional")
 306   public void testTred_reproducible()
 307   {
 308     /*
 309      * make a pseudo-random symmetric matrix as required for tred/tqli
 310      */
 311     int rows = 10;
 312     int cols = rows;
 313     double[][] d = getSparseValues(rows, cols, 3);
 314
 315     /*
 316      * make a copy of the values so m1, m2 are not
 317      * sharing arrays!
 318      */
 319     double[][] d1 = new double[rows][cols];
 320     for (int row = 0; row < rows; row++)
 321     {
 322       for (int col = 0; col < cols; col++)
 323       {
 324         d1[row][col] = d[row][col];
 325       }
 326     }
 327     Matrix m1 = new Matrix(d);
 328     Matrix m2 = new Matrix(d1);
 329     assertMatricesMatch(m1, m2); // sanity check
 330     m1.tred();
 331     m2.tred();
 332     assertMatricesMatch(m1, m2);
 333   }
 334
 335   private void assertMatricesMatch(MatrixI m1, MatrixI m2)
 336   {
 337     if (m1.height() != m2.height())
 338     {
 339       fail("height mismatch");
 340     }
 341     if (m1.width() != m2.width())
 342     {
 343       fail("width mismatch");
 344     }
 345     for (int row = 0; row < m1.height(); row++)
 346     {
 347       for (int col = 0; col < m1.width(); col++)
 348       {
 349         double v2 = m2.getValue(row, col);
 350         double v1 = m1.getValue(row, col);
 351         if (Math.abs(v1 - v2) > DELTA)
 352         {
 353           fail(String.format("At [%d, %d] %f != %f", row, col, v1, v2));
 354         }
 355       }
 356     }
 357     ArrayAsserts.assertArrayEquals(m1.getD(), m2.getD(), 0.00001d);
 358     ArrayAsserts.assertArrayEquals(m1.getE(), m2.getE(), 0.00001d);
 359   }
 360 }