test/jalview/math/MatrixTest.java

   1 package jalview.math;
   2
   3 import static org.testng.Assert.assertEquals;
   4 import static org.testng.Assert.assertFalse;
   5 import static org.testng.Assert.assertNotSame;
   6 import static org.testng.Assert.assertNull;
   7 import static org.testng.Assert.assertTrue;
   8 import static org.testng.Assert.fail;
   9
  10 import java.util.Arrays;
  11 import java.util.Random;
  12
  13 import org.testng.annotations.Test;
  14 import org.testng.internal.junit.ArrayAsserts;
  15
  16 public class MatrixTest
  17 {
  18   final static double DELTA = 0.000001d;
  19
  20   @Test(groups = "Timing")
  21   public void testPreMultiply_timing()
  22   {
  23     int rows = 50; // increase to stress test timing
  24     int cols = 100;
  25     double[][] d1 = new double[rows][cols];
  26     double[][] d2 = new double[cols][rows];
  27     Matrix m1 = new Matrix(d1);
  28     Matrix m2 = new Matrix(d2);
  29     long start = System.currentTimeMillis();
  30     m1.preMultiply(m2);
  31     long elapsed = System.currentTimeMillis() - start;
  32     System.out.println(rows + "x" + cols
  33             + " multiplications of double took " + elapsed + "ms");
  34   }
  35
  36   @Test(groups = "Functional")
  37   public void testPreMultiply()
  38   {
  39     Matrix m1 = new Matrix(new double[][] { { 2, 3, 4 } }); // 1x3
  40     Matrix m2 = new Matrix(new double[][] { { 5 }, { 6 }, { 7 } }); // 3x1
  41
  42     /*
  43      * 1x3 times 3x1 is 1x1
  44      * 2x5 + 3x6 + 4*7 =  56
  45      */
  46     MatrixI m3 = m2.preMultiply(m1);
  47     assertEquals(m3.height(), 1);
  48     assertEquals(m3.width(), 1);
  49     assertEquals(m3.getValue(0, 0), 56d);
  50
  51     /*
  52      * 3x1 times 1x3 is 3x3
  53      */
  54     m3 = m1.preMultiply(m2);
  55     assertEquals(m3.height(), 3);
  56     assertEquals(m3.width(), 3);
  57     assertEquals(Arrays.toString(m3.getRow(0)), "[10.0, 15.0, 20.0]");
  58     assertEquals(Arrays.toString(m3.getRow(1)), "[12.0, 18.0, 24.0]");
  59     assertEquals(Arrays.toString(m3.getRow(2)), "[14.0, 21.0, 28.0]");
  60   }
  61
  62   @Test(
  63     groups = "Functional",
  64     expectedExceptions = { IllegalArgumentException.class })
  65   public void testPreMultiply_tooManyColumns()
  66   {
  67     Matrix m1 = new Matrix(new double[][] { { 2, 3, 4 }, { 3, 4, 5 } }); // 2x3
  68
  69     /*
  70      * 2x3 times 2x3 invalid operation -
  71      * multiplier has more columns than multiplicand has rows
  72      */
  73     m1.preMultiply(m1);
  74     fail("Expected exception");
  75   }
  76
  77   @Test(
  78     groups = "Functional",
  79     expectedExceptions = { IllegalArgumentException.class })
  80   public void testPreMultiply_tooFewColumns()
  81   {
  82     Matrix m1 = new Matrix(new double[][] { { 2, 3, 4 }, { 3, 4, 5 } }); // 2x3
  83
  84     /*
  85      * 3x2 times 3x2 invalid operation -
  86      * multiplier has more columns than multiplicand has row
  87      */
  88     m1.preMultiply(m1);
  89     fail("Expected exception");
  90   }
  91
  92
  93   private boolean matrixEquals(Matrix m1, Matrix m2) {
  94     if (m1.width() != m2.width() || m1.height() != m2.height())
  95     {
  96       return false;
  97     }
  98     for (int i = 0; i < m1.height(); i++)
  99     {
 100       if (!Arrays.equals(m1.getRow(i), m2.getRow(i)))
 101       {
 102         return false;
 103       }
 104     }
 105     return true;
 106   }
 107
 108   @Test(groups = "Functional")
 109   public void testPostMultiply()
 110   {
 111     /*
 112      * Square matrices
 113      * (2 3) . (10   100)
 114      * (4 5)   (1000 10000)
 115      * =
 116      * (3020 30200)
 117      * (5040 50400)
 118      */
 119     MatrixI m1 = new Matrix(new double[][] { { 2, 3 }, { 4, 5 } });
 120     MatrixI m2 = new Matrix(new double[][] { { 10, 100 }, { 1000, 10000 } });
 121     MatrixI m3 = m1.postMultiply(m2);
 122     assertEquals(Arrays.toString(m3.getRow(0)), "[3020.0, 30200.0]");
 123     assertEquals(Arrays.toString(m3.getRow(1)), "[5040.0, 50400.0]");
 124
 125     /*
 126      * also check m2.preMultiply(m1) - should be same as m1.postMultiply(m2)
 127      */
 128     m3 = m2.preMultiply(m1);
 129     assertEquals(Arrays.toString(m3.getRow(0)), "[3020.0, 30200.0]");
 130     assertEquals(Arrays.toString(m3.getRow(1)), "[5040.0, 50400.0]");
 131
 132     /*
 133      * m1 has more rows than columns
 134      * (2).(10 100 1000) = (20 200 2000)
 135      * (3)                 (30 300 3000)
 136      */
 137     m1 = new Matrix(new double[][] { { 2 }, { 3 } });
 138     m2 = new Matrix(new double[][] { { 10, 100, 1000 } });
 139     m3 = m1.postMultiply(m2);
 140     assertEquals(m3.height(), 2);
 141     assertEquals(m3.width(), 3);
 142     assertEquals(Arrays.toString(m3.getRow(0)), "[20.0, 200.0, 2000.0]");
 143     assertEquals(Arrays.toString(m3.getRow(1)), "[30.0, 300.0, 3000.0]");
 144     m3 = m2.preMultiply(m1);
 145     assertEquals(m3.height(), 2);
 146     assertEquals(m3.width(), 3);
 147     assertEquals(Arrays.toString(m3.getRow(0)), "[20.0, 200.0, 2000.0]");
 148     assertEquals(Arrays.toString(m3.getRow(1)), "[30.0, 300.0, 3000.0]");
 149
 150     /*
 151      * m1 has more columns than rows
 152      * (2 3 4) . (5 4) = (56 25)
 153      *           (6 3)
 154      *           (7 2)
 155      * [0, 0] = 2*5 + 3*6 + 4*7 = 56
 156      * [0, 1] = 2*4 + 3*3 + 4*2 = 25
 157      */
 158     m1 = new Matrix(new double[][] { { 2, 3, 4 } });
 159     m2 = new Matrix(new double[][] { { 5, 4 }, { 6, 3 }, { 7, 2 } });
 160     m3 = m1.postMultiply(m2);
 161     assertEquals(m3.height(), 1);
 162     assertEquals(m3.width(), 2);
 163     assertEquals(m3.getRow(0)[0], 56d);
 164     assertEquals(m3.getRow(0)[1], 25d);
 165
 166     /*
 167      * and check premultiply equivalent
 168      */
 169     m3 = m2.preMultiply(m1);
 170     assertEquals(m3.height(), 1);
 171     assertEquals(m3.width(), 2);
 172     assertEquals(m3.getRow(0)[0], 56d);
 173     assertEquals(m3.getRow(0)[1], 25d);
 174   }
 175
 176   @Test(groups = "Functional")
 177   public void testCopy()
 178   {
 179     Random r = new Random();
 180     int rows = 5;
 181     int cols = 11;
 182     double[][] in = new double[rows][cols];
 183
 184     for (int i = 0; i < rows; i++)
 185     {
 186       for (int j = 0; j < cols; j++)
 187       {
 188         in[i][j] = r.nextDouble();
 189       }
 190     }
 191     Matrix m1 = new Matrix(in);
 192
 193     Matrix m2 = (Matrix) m1.copy();
 194     assertNotSame(m1, m2);
 195     assertTrue(matrixEquals(m1, m2));
 196     assertNull(m2.d);
 197     assertNull(m2.e);
 198
 199     /*
 200      * now add d and e vectors and recopy
 201      */
 202     m1.d = Arrays.copyOf(in[2], in[2].length);
 203     m1.e = Arrays.copyOf(in[4], in[4].length);
 204     m2 = (Matrix) m1.copy();
 205     assertNotSame(m2.d, m1.d);
 206     assertNotSame(m2.e, m1.e);
 207     assertEquals(m2.d, m1.d);
 208     assertEquals(m2.e, m1.e);
 209   }
 210
 211   /**
 212    * main method extracted from Matrix
 213    *
 214    * @param args
 215    */
 216   public static void main(String[] args) throws Exception
 217   {
 218     int n = Integer.parseInt(args[0]);
 219     double[][] in = new double[n][n];
 220
 221     for (int i = 0; i < n; i++)
 222     {
 223       for (int j = 0; j < n; j++)
 224       {
 225         in[i][j] = Math.random();
 226       }
 227     }
 228
 229     Matrix origmat = new Matrix(in);
 230
 231     // System.out.println(" --- Original matrix ---- ");
 232     // / origmat.print(System.out);
 233     // System.out.println();
 234     // System.out.println(" --- transpose matrix ---- ");
 235     MatrixI trans = origmat.transpose();
 236
 237     // trans.print(System.out);
 238     // System.out.println();
 239     // System.out.println(" --- OrigT * Orig ---- ");
 240     MatrixI symm = trans.postMultiply(origmat);
 241
 242     // symm.print(System.out);
 243     // System.out.println();
 244     // Copy the symmetric matrix for later
 245     // Matrix origsymm = symm.copy();
 246
 247     // This produces the tridiagonal transformation matrix
 248     // long tstart = System.currentTimeMillis();
 249     symm.tred();
 250
 251     // long tend = System.currentTimeMillis();
 252
 253     // System.out.println("Time take for tred = " + (tend-tstart) + "ms");
 254     // System.out.println(" ---Tridiag transform matrix ---");
 255     // symm.print(System.out);
 256     // System.out.println();
 257     // System.out.println(" --- D vector ---");
 258     // symm.printD(System.out);
 259     // System.out.println();
 260     // System.out.println(" --- E vector ---");
 261     // symm.printE(System.out);
 262     // System.out.println();
 263     // Now produce the diagonalization matrix
 264     // tstart = System.currentTimeMillis();
 265     symm.tqli();
 266     // tend = System.currentTimeMillis();
 267
 268     // System.out.println("Time take for tqli = " + (tend-tstart) + " ms");
 269     // System.out.println(" --- New diagonalization matrix ---");
 270     // symm.print(System.out);
 271     // System.out.println();
 272     // System.out.println(" --- D vector ---");
 273     // symm.printD(System.out);
 274     // System.out.println();
 275     // System.out.println(" --- E vector ---");
 276     // symm.printE(System.out);
 277     // System.out.println();
 278     // System.out.println(" --- First eigenvector --- ");
 279     // double[] eigenv = symm.getColumn(0);
 280     // for (int i=0; i < eigenv.length;i++) {
 281     // Format.print(System.out,"%15.4f",eigenv[i]);
 282     // }
 283     // System.out.println();
 284     // double[] neigenv = origsymm.vectorPostMultiply(eigenv);
 285     // for (int i=0; i < neigenv.length;i++) {
 286     // Format.print(System.out,"%15.4f",neigenv[i]/symm.d[0]);
 287     // }
 288     // System.out.println();
 289   }
 290
 291   @Test(groups = "Timing")
 292   public void testSign()
 293   {
 294     assertEquals(Matrix.sign(-1, -2), -1d);
 295     assertEquals(Matrix.sign(-1, 2), 1d);
 296     assertEquals(Matrix.sign(-1, 0), 1d);
 297     assertEquals(Matrix.sign(1, -2), -1d);
 298     assertEquals(Matrix.sign(1, 2), 1d);
 299     assertEquals(Matrix.sign(1, 0), 1d);
 300   }
 301
 302   /**
 303    * Helper method to make values for a sparse, pseudo-random symmetric matrix
 304    *
 305    * @param rows
 306    * @param cols
 307    * @param occupancy
 308    *          one in 'occupancy' entries will be non-zero
 309    * @return
 310    */
 311   public double[][] getSparseValues(int rows, int cols, int occupancy)
 312   {
 313     Random r = new Random(1729);
 314
 315     /*
 316      * generate whole number values between -12 and +12
 317      * (to mimic score matrices used in Jalview)
 318      */
 319     double[][] d = new double[rows][cols];
 320     int m = 0;
 321     for (int i = 0; i < rows; i++)
 322     {
 323       if (++m % occupancy == 0)
 324       {
 325         d[i][i] = r.nextInt() % 13; // diagonal
 326       }
 327       for (int j = 0; j < i; j++)
 328       {
 329         if (++m % occupancy == 0)
 330         {
 331           d[i][j] = r.nextInt() % 13;
 332           d[j][i] = d[i][j];
 333         }
 334       }
 335     }
 336     return d;
 337
 338   }
 339
 340   /**
 341    * Verify that the results of method tred() are the same if the calculation is
 342    * redone
 343    */
 344   @Test(groups = "Functional")
 345   public void testTred_reproducible()
 346   {
 347     /*
 348      * make a pseudo-random symmetric matrix as required for tred/tqli
 349      */
 350     int rows = 10;
 351     int cols = rows;
 352     double[][] d = getSparseValues(rows, cols, 3);
 353
 354     /*
 355      * make a copy of the values so m1, m2 are not
 356      * sharing arrays!
 357      */
 358     double[][] d1 = new double[rows][cols];
 359     for (int row = 0; row < rows; row++)
 360     {
 361       for (int col = 0; col < cols; col++)
 362       {
 363         d1[row][col] = d[row][col];
 364       }
 365     }
 366     Matrix m1 = new Matrix(d);
 367     Matrix m2 = new Matrix(d1);
 368     assertMatricesMatch(m1, m2); // sanity check
 369     m1.tred();
 370     m2.tred();
 371     assertMatricesMatch(m1, m2);
 372   }
 373
 374   public static void assertMatricesMatch(MatrixI m1, MatrixI m2)
 375   {
 376     if (m1.height() != m2.height())
 377     {
 378       fail("height mismatch");
 379     }
 380     if (m1.width() != m2.width())
 381     {
 382       fail("width mismatch");
 383     }
 384     for (int row = 0; row < m1.height(); row++)
 385     {
 386       for (int col = 0; col < m1.width(); col++)
 387       {
 388         double v2 = m2.getValue(row, col);
 389         double v1 = m1.getValue(row, col);
 390         if (Math.abs(v1 - v2) > DELTA)
 391         {
 392           fail(String.format("At [%d, %d] %f != %f", row, col, v1, v2));
 393         }
 394       }
 395     }
 396     ArrayAsserts.assertArrayEquals("D vector", m1.getD(), m2.getD(),
 397             0.00001d);
 398     ArrayAsserts.assertArrayEquals("E vector", m1.getE(), m2.getE(),
 399             0.00001d);
 400   }
 401
 402   @Test(groups = "Functional")
 403   public void testFindMinMax()
 404   {
 405     /*
 406      * empty matrix case
 407      */
 408     Matrix m = new Matrix(new double[][] { {} });
 409     assertNull(m.findMinMax());
 410
 411     /*
 412      * normal case
 413      */
 414     double[][] vals = new double[2][];
 415     vals[0] = new double[] {7d, 1d, -2.3d};
 416     vals[1] = new double[] {-12d, 94.3d, -102.34d};
 417     m = new Matrix(vals);
 418     double[] minMax = m.findMinMax();
 419     assertEquals(minMax[0], -102.34d);
 420     assertEquals(minMax[1], 94.3d);
 421   }
 422
 423   @Test(groups = { "Functional", "Timing" })
 424   public void testFindMinMax_timing()
 425   {
 426     Random r = new Random();
 427     int size = 1000; // increase to stress test timing
 428     double[][] vals = new double[size][size];
 429     double max = -Double.MAX_VALUE;
 430     double min = Double.MAX_VALUE;
 431     for (int i = 0; i < size; i++)
 432     {
 433       vals[i] = new double[size];
 434       for (int j = 0; j < size; j++)
 435       {
 436         // use nextLong rather than nextDouble to include negative values
 437         double d = r.nextLong();
 438         if (d > max)
 439         {
 440           max = d;
 441         }
 442         if (d < min)
 443         {
 444           min = d;
 445         }
 446         vals[i][j] = d;
 447       }
 448     }
 449     Matrix m = new Matrix(vals);
 450     long now = System.currentTimeMillis();
 451     double[] minMax = m.findMinMax();
 452     System.out.println(String.format("findMinMax for %d x %d took %dms",
 453             size, size, (System.currentTimeMillis() - now)));
 454     assertEquals(minMax[0], min);
 455     assertEquals(minMax[1], max);
 456   }
 457
 458   /**
 459    * Test range reversal with maximum value becoming zero
 460    */
 461   @Test(groups = "Functional")
 462   public void testReverseRange_maxToZero()
 463   {
 464     Matrix m1 = new Matrix(
 465             new double[][] { { 2, 3.5, 4 }, { -3.4, 4, 15 } });
 466
 467     /*
 468      * subtract all from max: range -3.4 to 15 becomes 18.4 to 0
 469      */
 470     m1.reverseRange(true);
 471     assertEquals(m1.getValue(0, 0), 13d, DELTA);
 472     assertEquals(m1.getValue(0, 1), 11.5d, DELTA);
 473     assertEquals(m1.getValue(0, 2), 11d, DELTA);
 474     assertEquals(m1.getValue(1, 0), 18.4d, DELTA);
 475     assertEquals(m1.getValue(1, 1), 11d, DELTA);
 476     assertEquals(m1.getValue(1, 2), 0d, DELTA);
 477
 478     /*
 479      * repeat operation - range is now 0 to 18.4
 480      */
 481     m1.reverseRange(true);
 482     assertEquals(m1.getValue(0, 0), 5.4d, DELTA);
 483     assertEquals(m1.getValue(0, 1), 6.9d, DELTA);
 484     assertEquals(m1.getValue(0, 2), 7.4d, DELTA);
 485     assertEquals(m1.getValue(1, 0), 0d, DELTA);
 486     assertEquals(m1.getValue(1, 1), 7.4d, DELTA);
 487     assertEquals(m1.getValue(1, 2), 18.4d, DELTA);
 488   }
 489
 490   /**
 491    * Test range reversal with minimum and maximum values swapped
 492    */
 493   @Test(groups = "Functional")
 494   public void testReverseRange_swapMinMax()
 495   {
 496     Matrix m1 = new Matrix(
 497             new double[][] { { 2, 3.5, 4 }, { -3.4, 4, 15 } });
 498
 499     /*
 500      * swap all values in min-max range
 501      * = subtract from (min + max = 11.6)
 502      * range -3.4 to 15 becomes 18.4 to -3.4
 503      */
 504     m1.reverseRange(false);
 505     assertEquals(m1.getValue(0, 0), 9.6d, DELTA);
 506     assertEquals(m1.getValue(0, 1), 8.1d, DELTA);
 507     assertEquals(m1.getValue(0, 2), 7.6d, DELTA);
 508     assertEquals(m1.getValue(1, 0), 15d, DELTA);
 509     assertEquals(m1.getValue(1, 1), 7.6d, DELTA);
 510     assertEquals(m1.getValue(1, 2), -3.4d, DELTA);
 511
 512     /*
 513      * repeat operation - original values restored
 514      */
 515     m1.reverseRange(false);
 516     assertEquals(m1.getValue(0, 0), 2d, DELTA);
 517     assertEquals(m1.getValue(0, 1), 3.5d, DELTA);
 518     assertEquals(m1.getValue(0, 2), 4d, DELTA);
 519     assertEquals(m1.getValue(1, 0), -3.4d, DELTA);
 520     assertEquals(m1.getValue(1, 1), 4d, DELTA);
 521     assertEquals(m1.getValue(1, 2), 15d, DELTA);
 522   }
 523
 524   @Test(groups = "Functional")
 525   public void testMultiply()
 526   {
 527     Matrix m = new Matrix(new double[][] { { 2, 3.5, 4 }, { -3.4, 4, 15 } });
 528     m.multiply(2d);
 529     assertEquals(m.getValue(0, 0), 4d, DELTA);
 530     assertEquals(m.getValue(0, 1), 7d, DELTA);
 531     assertEquals(m.getValue(0, 2), 8d, DELTA);
 532     assertEquals(m.getValue(1, 0), -6.8d, DELTA);
 533     assertEquals(m.getValue(1, 1), 8d, DELTA);
 534     assertEquals(m.getValue(1, 2), 30d, DELTA);
 535   }
 536
 537   @Test(groups = "Functional")
 538   public void testConstructor()
 539   {
 540     double[][] values = new double[][] { { 1, 2, 3 }, { 4, 5, 6 } };
 541     Matrix m = new Matrix(values);
 542     assertEquals(m.getValue(0, 0), 1d, DELTA);
 543
 544     /*
 545      * verify the matrix has a copy of the original array
 546      */
 547     assertNotSame(values[0], m.getRow(0));
 548     values[0][0] = -1d;
 549     assertEquals(m.getValue(0, 0), 1d, DELTA); // unchanged
 550   }
 551
 552   @Test(groups = "Functional")
 553   public void testEquals()
 554   {
 555     double[][] values = new double[][] { { 1, 2, 3 }, { 4, 5, 6 } };
 556     Matrix m1 = new Matrix(values);
 557     double[][] values2 = new double[][] { { 1, 2, 3 }, { 4, 5, 6 } };
 558     Matrix m2 = new Matrix(values2);
 559
 560     double delta = 0.0001d;
 561     assertTrue(m1.equals(m1, delta));
 562     assertTrue(m1.equals(m2, delta));
 563     assertTrue(m2.equals(m1, delta));
 564
 565     double[][] values3 = new double[][] { { 1, 2, 3 }, { 4, 5, 7 } };
 566     m2 = new Matrix(values3);
 567     assertFalse(m1.equals(m2, delta));
 568     assertFalse(m2.equals(m1, delta));
 569
 570     // must be same shape
 571     values2 = new double[][] { { 1, 2, 3 } };
 572     m2 = new Matrix(values2);
 573     assertFalse(m2.equals(m1, delta));
 574
 575     assertFalse(m1.equals(null, delta));
 576   }
 577 }