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