src/jalview/math/Matrix.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.math;
  22
  23 import jalview.util.Format;
  24 import jalview.util.MessageManager;
  25
  26 import java.io.PrintStream;
  27
  28 /**
  29  * DOCUMENT ME!
  30  *
  31  * @author $author$
  32  * @version $Revision$
  33  */
  34 public class Matrix implements MatrixI
  35 {
  36   /*
  37    * the [row][column] values in the matrix
  38    */
  39   private double[][] value;
  40
  41   /*
  42    * the number of rows
  43    */
  44   protected int rows;
  45
  46   /*
  47    * the number of columns
  48    */
  49   protected int cols;
  50
  51   /** DOCUMENT ME!! */
  52   protected double[] d; // Diagonal
  53
  54   /** DOCUMENT ME!! */
  55   protected double[] e; // off diagonal
  56
  57   /**
  58    * maximum number of iterations for tqli
  59    *
  60    */
  61   private static final int maxIter = 45; // fudge - add 15 iterations, just in
  62                                          // case
  63
  64   /**
  65    * Default constructor
  66    */
  67   public Matrix()
  68   {
  69
  70   }
  71
  72   /**
  73    * Creates a new Matrix object. For example
  74    *
  75    * <pre>
  76    *   new Matrix(new double[][] {{2, 3, 4}, {5, 6, 7})
  77    * constructs
  78    *   (2 3 4)
  79    *   (5 6 7)
  80    * </pre>
  81    *
  82    * Note that ragged arrays (with not all rows, or columns, of the same
  83    * length), are not supported by this class. They can be constructed, but
  84    * results of operations on them are undefined and may throw exceptions.
  85    *
  86    * @param values
  87    *          the matrix values in row-major order
  88    */
  89   public Matrix(double[][] values)
  90   {
  91     this.rows = values.length;
  92     if (rows > 0)
  93     {
  94       this.cols = values[0].length;
  95     }
  96     this.value = values;
  97   }
  98
  99   /**
 100    * Returns a new matrix which is the transpose of this one
 101    *
 102    * @return DOCUMENT ME!
 103    */
 104   @Override
 105   public MatrixI transpose()
 106   {
 107     double[][] out = new double[cols][rows];
 108
 109     for (int i = 0; i < cols; i++)
 110     {
 111       for (int j = 0; j < rows; j++)
 112       {
 113         out[i][j] = value[j][i];
 114       }
 115     }
 116
 117     return new Matrix(out);
 118   }
 119
 120   /**
 121    * DOCUMENT ME!
 122    *
 123    * @param ps
 124    *          DOCUMENT ME!
 125    * @param format
 126    */
 127   @Override
 128   public void print(PrintStream ps, String format)
 129   {
 130     for (int i = 0; i < rows; i++)
 131     {
 132       for (int j = 0; j < cols; j++)
 133       {
 134         Format.print(ps, format, getValue(i, j));
 135       }
 136
 137       ps.println();
 138     }
 139   }
 140
 141   /**
 142    * Returns a new matrix which is the result of premultiplying this matrix by
 143    * the supplied argument. If this of size AxB (A rows and B columns), and the
 144    * argument is CxA (C rows and A columns), the result is of size CxB.
 145    *
 146    * @param in
 147    *
 148    * @return
 149    * @throws IllegalArgumentException
 150    *           if the number of columns in the pre-multiplier is not equal to
 151    *           the number of rows in the multiplicand (this)
 152    */
 153   @Override
 154   public MatrixI preMultiply(MatrixI in)
 155   {
 156     if (in.width() != rows)
 157     {
 158       throw new IllegalArgumentException("Can't pre-multiply " + this.rows
 159               + " rows by " + in.width() + " columns");
 160     }
 161     double[][] tmp = new double[in.height()][this.cols];
 162
 163     for (int i = 0; i < in.height(); i++)
 164     {
 165       for (int j = 0; j < this.cols; j++)
 166       {
 167         /*
 168          * result[i][j] is the vector product of
 169          * in.row[i] and this.column[j]
 170          */
 171         for (int k = 0; k < in.width(); k++)
 172         {
 173           tmp[i][j] += (in.getValue(i, k) * this.value[k][j]);
 174         }
 175       }
 176     }
 177
 178     return new Matrix(tmp);
 179   }
 180
 181   /**
 182    *
 183    * @param in
 184    *
 185    * @return
 186    */
 187   public double[] vectorPostMultiply(double[] in)
 188   {
 189     double[] out = new double[in.length];
 190
 191     for (int i = 0; i < in.length; i++)
 192     {
 193       out[i] = 0.0;
 194
 195       for (int k = 0; k < in.length; k++)
 196       {
 197         out[i] += (value[i][k] * in[k]);
 198       }
 199     }
 200
 201     return out;
 202   }
 203
 204   /**
 205    * Returns a new matrix which is the result of postmultiplying this matrix by
 206    * the supplied argument. If this of size AxB (A rows and B columns), and the
 207    * argument is BxC (B rows and C columns), the result is of size AxC.
 208    * <p>
 209    * This method simply returns the result of in.preMultiply(this)
 210    *
 211    * @param in
 212    *
 213    * @return
 214    * @throws IllegalArgumentException
 215    *           if the number of rows in the post-multiplier is not equal to the
 216    *           number of columns in the multiplicand (this)
 217    * @see #preMultiply(Matrix)
 218    */
 219   @Override
 220   public MatrixI postMultiply(MatrixI in)
 221   {
 222     if (in.height() != this.cols)
 223     {
 224       throw new IllegalArgumentException("Can't post-multiply " + this.cols
 225               + " columns by " + in.height() + " rows");
 226     }
 227     return in.preMultiply(this);
 228   }
 229
 230   /**
 231    * Answers a new matrix with a copy of the values in this one
 232    *
 233    * @return
 234    */
 235   @Override
 236   public MatrixI copy()
 237   {
 238     double[][] newmat = new double[rows][cols];
 239
 240     for (int i = 0; i < rows; i++)
 241     {
 242       System.arraycopy(value[i], 0, newmat[i], 0, value[i].length);
 243     }
 244
 245     return new Matrix(newmat);
 246   }
 247
 248   /**
 249    * DOCUMENT ME!
 250    */
 251   @Override
 252   public void tred()
 253   {
 254     int n = rows;
 255     int k;
 256     int j;
 257     int i;
 258
 259     double scale;
 260     double hh;
 261     double h;
 262     double g;
 263     double f;
 264
 265     this.d = new double[rows];
 266     this.e = new double[rows];
 267
 268     for (i = n; i >= 2; i--)
 269     {
 270       final int l = i - 1;
 271       h = 0.0;
 272       scale = 0.0;
 273
 274       if (l > 1)
 275       {
 276         for (k = 1; k <= l; k++)
 277         {
 278           // double v = Math.abs(value[i - 1][k - 1]);
 279           double v = Math.abs(getValue(i - 1, k - 1));
 280           scale += v;
 281         }
 282
 283         if (scale == 0.0)
 284         {
 285           // e[i - 1] = value[i - 1][l - 1];
 286           e[i - 1] = getValue(i - 1, l - 1);
 287         }
 288         else
 289         {
 290           for (k = 1; k <= l; k++)
 291           {
 292             // value[i - 1][k - 1] /= scale;
 293             // h += (value[i - 1][k - 1] * value[i - 1][k - 1]);
 294             double v = divideValue(i - 1, k - 1, scale);
 295             h += v * v;
 296           }
 297
 298           // f = value[i - 1][l - 1];
 299           f = getValue(i - 1, l - 1);
 300
 301           if (f > 0)
 302           {
 303             g = -1.0 * Math.sqrt(h);
 304           }
 305           else
 306           {
 307             g = Math.sqrt(h);
 308           }
 309
 310           e[i - 1] = scale * g;
 311           h -= (f * g);
 312           // value[i - 1][l - 1] = f - g;
 313           setValue(i - 1, l - 1, f - g);
 314           // System.out.println(String.format("%d %d %f %f %f %f %f %.5e", i,
 315           // l,
 316           // scale, f, g, h, getValue(i - 1, l - 1), checksum()));
 317           f = 0.0;
 318
 319           for (j = 1; j <= l; j++)
 320           {
 321             // value[j - 1][i - 1] = value[i - 1][j - 1] / h;
 322             double val = getValue(i - 1, j - 1) / h;
 323             setValue(j - 1, i - 1, val);
 324             g = 0.0;
 325
 326             for (k = 1; k <= j; k++)
 327             {
 328               // g += (value[j - 1][k - 1] * value[i - 1][k - 1]);
 329               g += (getValue(j - 1, k - 1) * getValue(i - 1, k - 1));
 330             }
 331
 332             for (k = j + 1; k <= l; k++)
 333             {
 334               // g += (value[k - 1][j - 1] * value[i - 1][k - 1]);
 335               g += (getValue(k - 1, j - 1) * getValue(i - 1, k - 1));
 336             }
 337
 338             e[j - 1] = g / h;
 339             // f += (e[j - 1] * value[i - 1][j - 1]);
 340             f += (e[j - 1] * getValue(i - 1, j - 1));
 341           }
 342
 343           hh = f / (h + h);
 344
 345           for (j = 1; j <= l; j++)
 346           {
 347             // f = value[i - 1][j - 1];
 348             f = getValue(i - 1, j - 1);
 349             g = e[j - 1] - (hh * f);
 350             e[j - 1] = g;
 351
 352             for (k = 1; k <= j; k++)
 353             {
 354               // value[j - 1][k - 1] -= ((f * e[k - 1]) + (g * value[i - 1][k -
 355               // 1]));
 356               double val = (f * e[k - 1]) + (g * getValue(i - 1, k - 1));
 357               addValue(j - 1, k - 1, -val);
 358             }
 359           }
 360         }
 361       }
 362       else
 363       {
 364         // e[i - 1] = value[i - 1][l - 1];
 365         e[i - 1] = getValue(i - 1, l - 1);
 366       }
 367
 368       d[i - 1] = h;
 369     }
 370
 371     d[0] = 0.0;
 372     e[0] = 0.0;
 373
 374     for (i = 1; i <= n; i++)
 375     {
 376       final int l = i - 1;
 377
 378       if (d[i - 1] != 0.0)
 379       {
 380         for (j = 1; j <= l; j++)
 381         {
 382           g = 0.0;
 383
 384           for (k = 1; k <= l; k++)
 385           {
 386             // g += (value[i - 1][k - 1] * value[k - 1][j - 1]);
 387             g += (getValue(i - 1, k - 1) * getValue(k - 1, j - 1));
 388           }
 389
 390           for (k = 1; k <= l; k++)
 391           {
 392             // value[k - 1][j - 1] -= (g * value[k - 1][i - 1]);
 393             addValue(k - 1, j - 1, -(g * getValue(k - 1, i - 1)));
 394           }
 395         }
 396       }
 397
 398       // d[i - 1] = value[i - 1][i - 1];
 399       // value[i - 1][i - 1] = 1.0;
 400       d[i - 1] = getValue(i - 1, i - 1);
 401       setValue(i - 1, i - 1, 1.0);
 402
 403       for (j = 1; j <= l; j++)
 404       {
 405         // value[j - 1][i - 1] = 0.0;
 406         // value[i - 1][j - 1] = 0.0;
 407         setValue(j - 1, i - 1, 0.0);
 408         setValue(i - 1, j - 1, 0.0);
 409       }
 410     }
 411   }
 412
 413   /**
 414    * Adds f to the value at [i, j] and returns the new value
 415    *
 416    * @param i
 417    * @param j
 418    * @param f
 419    */
 420   protected double addValue(int i, int j, double f)
 421   {
 422     double v = value[i][j] + f;
 423     value[i][j] = v;
 424     return v;
 425   }
 426
 427   /**
 428    * Divides the value at [i, j] by divisor and returns the new value. If d is
 429    * zero, returns the unchanged value.
 430    *
 431    * @param i
 432    * @param j
 433    * @param divisor
 434    * @return
 435    */
 436   protected double divideValue(int i, int j, double divisor)
 437   {
 438     if (divisor == 0d)
 439     {
 440       return getValue(i, j);
 441     }
 442     double v = value[i][j];
 443     v = v / divisor;
 444     value[i][j] = v;
 445     return v;
 446   }
 447
 448   /**
 449    * DOCUMENT ME!
 450    */
 451   @Override
 452   public void tqli() throws Exception
 453   {
 454     int n = rows;
 455
 456     int m;
 457     int l;
 458     int iter;
 459     int i;
 460     int k;
 461     double s;
 462     double r;
 463     double p;
 464
 465     double g;
 466     double f;
 467     double dd;
 468     double c;
 469     double b;
 470
 471     for (i = 2; i <= n; i++)
 472     {
 473       e[i - 2] = e[i - 1];
 474     }
 475
 476     e[n - 1] = 0.0;
 477
 478     for (l = 1; l <= n; l++)
 479     {
 480       iter = 0;
 481
 482       do
 483       {
 484         for (m = l; m <= (n - 1); m++)
 485         {
 486           dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
 487
 488           if ((Math.abs(e[m - 1]) + dd) == dd)
 489           {
 490             break;
 491           }
 492         }
 493
 494         if (m != l)
 495         {
 496           iter++;
 497
 498           if (iter == maxIter)
 499           {
 500             throw new Exception(MessageManager.formatMessage(
 501                     "exception.matrix_too_many_iteration", new String[] {
 502                         "tqli", Integer.valueOf(maxIter).toString() }));
 503           }
 504           else
 505           {
 506             // System.out.println("Iteration " + iter);
 507           }
 508
 509           g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
 510           r = Math.sqrt((g * g) + 1.0);
 511           g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
 512           c = 1.0;
 513           s = c;
 514           p = 0.0;
 515
 516           for (i = m - 1; i >= l; i--)
 517           {
 518             f = s * e[i - 1];
 519             b = c * e[i - 1];
 520
 521             if (Math.abs(f) >= Math.abs(g))
 522             {
 523               c = g / f;
 524               r = Math.sqrt((c * c) + 1.0);
 525               e[i] = f * r;
 526               s = 1.0 / r;
 527               c *= s;
 528             }
 529             else
 530             {
 531               s = f / g;
 532               r = Math.sqrt((s * s) + 1.0);
 533               e[i] = g * r;
 534               c = 1.0 / r;
 535               s *= c;
 536             }
 537
 538             g = d[i] - p;
 539             r = ((d[i - 1] - g) * s) + (2.0 * c * b);
 540             p = s * r;
 541             d[i] = g + p;
 542             g = (c * r) - b;
 543
 544             for (k = 1; k <= n; k++)
 545             {
 546               // f = value[k - 1][i];
 547               // value[k - 1][i] = (s * value[k - 1][i - 1]) + (c * f);
 548               // value[k - 1][i - 1] = (c * value[k - 1][i - 1]) - (s * f);
 549               f = getValue(k - 1, i);
 550               setValue(k - 1, i, (s * getValue(k - 1, i - 1)) + (c * f));
 551               setValue(k - 1, i - 1, (c * getValue(k - 1, i - 1)) - (s * f));
 552             }
 553           }
 554
 555           d[l - 1] = d[l - 1] - p;
 556           e[l - 1] = g;
 557           e[m - 1] = 0.0;
 558         }
 559       } while (m != l);
 560     }
 561   }
 562
 563   @Override
 564   public double getValue(int i, int j)
 565   {
 566     return value[i][j];
 567   }
 568
 569   public void setValue(int i, int j, double val)
 570   {
 571     value[i][j] = val;
 572   }
 573
 574   /**
 575    * DOCUMENT ME!
 576    */
 577   public void tred2()
 578   {
 579     int n = rows;
 580     int l;
 581     int k;
 582     int j;
 583     int i;
 584
 585     double scale;
 586     double hh;
 587     double h;
 588     double g;
 589     double f;
 590
 591     this.d = new double[rows];
 592     this.e = new double[rows];
 593
 594     for (i = n - 1; i >= 1; i--)
 595     {
 596       l = i - 1;
 597       h = 0.0;
 598       scale = 0.0;
 599
 600       if (l > 0)
 601       {
 602         for (k = 0; k < l; k++)
 603         {
 604           scale += Math.abs(value[i][k]);
 605         }
 606
 607         if (scale == 0.0)
 608         {
 609           e[i] = value[i][l];
 610         }
 611         else
 612         {
 613           for (k = 0; k < l; k++)
 614           {
 615             value[i][k] /= scale;
 616             h += (value[i][k] * value[i][k]);
 617           }
 618
 619           f = value[i][l];
 620
 621           if (f > 0)
 622           {
 623             g = -1.0 * Math.sqrt(h);
 624           }
 625           else
 626           {
 627             g = Math.sqrt(h);
 628           }
 629
 630           e[i] = scale * g;
 631           h -= (f * g);
 632           value[i][l] = f - g;
 633           f = 0.0;
 634
 635           for (j = 0; j < l; j++)
 636           {
 637             value[j][i] = value[i][j] / h;
 638             g = 0.0;
 639
 640             for (k = 0; k < j; k++)
 641             {
 642               g += (value[j][k] * value[i][k]);
 643             }
 644
 645             for (k = j; k < l; k++)
 646             {
 647               g += (value[k][j] * value[i][k]);
 648             }
 649
 650             e[j] = g / h;
 651             f += (e[j] * value[i][j]);
 652           }
 653
 654           hh = f / (h + h);
 655
 656           for (j = 0; j < l; j++)
 657           {
 658             f = value[i][j];
 659             g = e[j] - (hh * f);
 660             e[j] = g;
 661
 662             for (k = 0; k < j; k++)
 663             {
 664               value[j][k] -= ((f * e[k]) + (g * value[i][k]));
 665             }
 666           }
 667         }
 668       }
 669       else
 670       {
 671         e[i] = value[i][l];
 672       }
 673
 674       d[i] = h;
 675     }
 676
 677     d[0] = 0.0;
 678     e[0] = 0.0;
 679
 680     for (i = 0; i < n; i++)
 681     {
 682       l = i - 1;
 683
 684       if (d[i] != 0.0)
 685       {
 686         for (j = 0; j < l; j++)
 687         {
 688           g = 0.0;
 689
 690           for (k = 0; k < l; k++)
 691           {
 692             g += (value[i][k] * value[k][j]);
 693           }
 694
 695           for (k = 0; k < l; k++)
 696           {
 697             value[k][j] -= (g * value[k][i]);
 698           }
 699         }
 700       }
 701
 702       d[i] = value[i][i];
 703       value[i][i] = 1.0;
 704
 705       for (j = 0; j < l; j++)
 706       {
 707         value[j][i] = 0.0;
 708         value[i][j] = 0.0;
 709       }
 710     }
 711   }
 712
 713   /**
 714    * DOCUMENT ME!
 715    */
 716   public void tqli2() throws Exception
 717   {
 718     int n = rows;
 719
 720     int m;
 721     int l;
 722     int iter;
 723     int i;
 724     int k;
 725     double s;
 726     double r;
 727     double p;
 728     ;
 729
 730     double g;
 731     double f;
 732     double dd;
 733     double c;
 734     double b;
 735
 736     for (i = 2; i <= n; i++)
 737     {
 738       e[i - 2] = e[i - 1];
 739     }
 740
 741     e[n - 1] = 0.0;
 742
 743     for (l = 1; l <= n; l++)
 744     {
 745       iter = 0;
 746
 747       do
 748       {
 749         for (m = l; m <= (n - 1); m++)
 750         {
 751           dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
 752
 753           if ((Math.abs(e[m - 1]) + dd) == dd)
 754           {
 755             break;
 756           }
 757         }
 758
 759         if (m != l)
 760         {
 761           iter++;
 762
 763           if (iter == maxIter)
 764           {
 765             throw new Exception(MessageManager.formatMessage(
 766                     "exception.matrix_too_many_iteration", new String[] {
 767                         "tqli2", Integer.valueOf(maxIter).toString() }));
 768           }
 769           else
 770           {
 771             // System.out.println("Iteration " + iter);
 772           }
 773
 774           g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
 775           r = Math.sqrt((g * g) + 1.0);
 776           g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
 777           c = 1.0;
 778           s = c;
 779           p = 0.0;
 780
 781           for (i = m - 1; i >= l; i--)
 782           {
 783             f = s * e[i - 1];
 784             b = c * e[i - 1];
 785
 786             if (Math.abs(f) >= Math.abs(g))
 787             {
 788               c = g / f;
 789               r = Math.sqrt((c * c) + 1.0);
 790               e[i] = f * r;
 791               s = 1.0 / r;
 792               c *= s;
 793             }
 794             else
 795             {
 796               s = f / g;
 797               r = Math.sqrt((s * s) + 1.0);
 798               e[i] = g * r;
 799               c = 1.0 / r;
 800               s *= c;
 801             }
 802
 803             g = d[i] - p;
 804             r = ((d[i - 1] - g) * s) + (2.0 * c * b);
 805             p = s * r;
 806             d[i] = g + p;
 807             g = (c * r) - b;
 808
 809             for (k = 1; k <= n; k++)
 810             {
 811               f = value[k - 1][i];
 812               value[k - 1][i] = (s * value[k - 1][i - 1]) + (c * f);
 813               value[k - 1][i - 1] = (c * value[k - 1][i - 1]) - (s * f);
 814             }
 815           }
 816
 817           d[l - 1] = d[l - 1] - p;
 818           e[l - 1] = g;
 819           e[m - 1] = 0.0;
 820         }
 821       } while (m != l);
 822     }
 823   }
 824
 825   /**
 826    * Answers the first argument with the sign of the second argument
 827    *
 828    * @param a
 829    * @param b
 830    *
 831    * @return
 832    */
 833   static double sign(double a, double b)
 834   {
 835     if (b < 0)
 836     {
 837       return -Math.abs(a);
 838     }
 839     else
 840     {
 841       return Math.abs(a);
 842     }
 843   }
 844
 845   /**
 846    * Returns an array containing the values in the specified column
 847    *
 848    * @param col
 849    *
 850    * @return
 851    */
 852   public double[] getColumn(int col)
 853   {
 854     double[] out = new double[rows];
 855
 856     for (int i = 0; i < rows; i++)
 857     {
 858       out[i] = value[i][col];
 859     }
 860
 861     return out;
 862   }
 863
 864   /**
 865    * DOCUMENT ME!
 866    *
 867    * @param ps
 868    *          DOCUMENT ME!
 869    * @param format
 870    */
 871   @Override
 872   public void printD(PrintStream ps, String format)
 873   {
 874     for (int j = 0; j < rows; j++)
 875     {
 876       Format.print(ps, format, d[j]);
 877     }
 878   }
 879
 880   /**
 881    * DOCUMENT ME!
 882    *
 883    * @param ps
 884    *          DOCUMENT ME!
 885    * @param format TODO
 886    */
 887   @Override
 888   public void printE(PrintStream ps, String format)
 889   {
 890     for (int j = 0; j < rows; j++)
 891     {
 892       Format.print(ps, format, e[j]);
 893     }
 894   }
 895
 896   @Override
 897   public double[] getD()
 898   {
 899     return d;
 900   }
 901
 902   @Override
 903   public double[] getE()
 904   {
 905     return e;
 906   }
 907
 908   @Override
 909   public int height() {
 910     return rows;
 911   }
 912
 913   @Override
 914   public int width()
 915   {
 916     return cols;
 917   }
 918
 919   @Override
 920   public double[] getRow(int i)
 921   {
 922     double[] row = new double[cols];
 923     System.arraycopy(value[i], 0, row, 0, cols);
 924     return row;
 925   }
 926 }