test/jalview/math/MatrixTest.java

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