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(getValue(i - 1, k - 1));
 271           scale += v;
 272         }
 273
 274         if (scale == 0.0)
 275         {
 276           e[i - 1] = getValue(i - 1, l - 1);
 277         }
 278         else
 279         {
 280           for (k = 1; k <= l; k++)
 281           {
 282             double v = divideValue(i - 1, k - 1, scale);
 283             h += v * v;
 284           }
 285
 286           f = getValue(i - 1, l - 1);
 287
 288           if (f > 0)
 289           {
 290             g = -1.0 * Math.sqrt(h);
 291           }
 292           else
 293           {
 294             g = Math.sqrt(h);
 295           }
 296
 297           e[i - 1] = scale * g;
 298           h -= (f * g);
 299           setValue(i - 1, l - 1, f - g);
 300           f = 0.0;
 301
 302           for (j = 1; j <= l; j++)
 303           {
 304             double val = getValue(i - 1, j - 1) / h;
 305             setValue(j - 1, i - 1, val);
 306             g = 0.0;
 307
 308             for (k = 1; k <= j; k++)
 309             {
 310               g += (getValue(j - 1, k - 1) * getValue(i - 1, k - 1));
 311             }
 312
 313             for (k = j + 1; k <= l; k++)
 314             {
 315               g += (getValue(k - 1, j - 1) * getValue(i - 1, k - 1));
 316             }
 317
 318             e[j - 1] = g / h;
 319             f += (e[j - 1] * getValue(i - 1, j - 1));
 320           }
 321
 322           hh = f / (h + h);
 323
 324           for (j = 1; j <= l; j++)
 325           {
 326             f = getValue(i - 1, j - 1);
 327             g = e[j - 1] - (hh * f);
 328             e[j - 1] = g;
 329
 330             for (k = 1; k <= j; k++)
 331             {
 332               double val = (f * e[k - 1]) + (g * getValue(i - 1, k - 1));
 333               addValue(j - 1, k - 1, -val);
 334             }
 335           }
 336         }
 337       }
 338       else
 339       {
 340         e[i - 1] = getValue(i - 1, l - 1);
 341       }
 342
 343       d[i - 1] = h;
 344     }
 345
 346     d[0] = 0.0;
 347     e[0] = 0.0;
 348
 349     for (i = 1; i <= n; i++)
 350     {
 351       final int l = i - 1;
 352
 353       if (d[i - 1] != 0.0)
 354       {
 355         for (j = 1; j <= l; j++)
 356         {
 357           g = 0.0;
 358
 359           for (k = 1; k <= l; k++)
 360           {
 361             g += (getValue(i - 1, k - 1) * getValue(k - 1, j - 1));
 362           }
 363
 364           for (k = 1; k <= l; k++)
 365           {
 366             addValue(k - 1, j - 1, -(g * getValue(k - 1, i - 1)));
 367           }
 368         }
 369       }
 370
 371       d[i - 1] = getValue(i - 1, i - 1);
 372       setValue(i - 1, i - 1, 1.0);
 373
 374       for (j = 1; j <= l; j++)
 375       {
 376         setValue(j - 1, i - 1, 0.0);
 377         setValue(i - 1, j - 1, 0.0);
 378       }
 379     }
 380   }
 381
 382   /**
 383    * Adds f to the value at [i, j] and returns the new value
 384    *
 385    * @param i
 386    * @param j
 387    * @param f
 388    */
 389   protected double addValue(int i, int j, double f)
 390   {
 391     double v = value[i][j] + f;
 392     value[i][j] = v;
 393     return v;
 394   }
 395
 396   /**
 397    * Divides the value at [i, j] by divisor and returns the new value. If d is
 398    * zero, returns the unchanged value.
 399    *
 400    * @param i
 401    * @param j
 402    * @param divisor
 403    * @return
 404    */
 405   protected double divideValue(int i, int j, double divisor)
 406   {
 407     if (divisor == 0d)
 408     {
 409       return getValue(i, j);
 410     }
 411     double v = value[i][j];
 412     v = v / divisor;
 413     value[i][j] = v;
 414     return v;
 415   }
 416
 417   /**
 418    * DOCUMENT ME!
 419    */
 420   @Override
 421   public void tqli() throws Exception
 422   {
 423     int n = rows;
 424
 425     int m;
 426     int l;
 427     int iter;
 428     int i;
 429     int k;
 430     double s;
 431     double r;
 432     double p;
 433
 434     double g;
 435     double f;
 436     double dd;
 437     double c;
 438     double b;
 439
 440     for (i = 2; i <= n; i++)
 441     {
 442       e[i - 2] = e[i - 1];
 443     }
 444
 445     e[n - 1] = 0.0;
 446
 447     for (l = 1; l <= n; l++)
 448     {
 449       iter = 0;
 450
 451       do
 452       {
 453         for (m = l; m <= (n - 1); m++)
 454         {
 455           dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
 456
 457           if ((Math.abs(e[m - 1]) + dd) == dd)
 458           {
 459             break;
 460           }
 461         }
 462
 463         if (m != l)
 464         {
 465           iter++;
 466
 467           if (iter == maxIter)
 468           {
 469             throw new Exception(MessageManager.formatMessage(
 470                     "exception.matrix_too_many_iteration", new String[] {
 471                         "tqli", Integer.valueOf(maxIter).toString() }));
 472           }
 473           else
 474           {
 475             // System.out.println("Iteration " + iter);
 476           }
 477
 478           g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
 479           r = Math.sqrt((g * g) + 1.0);
 480           g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
 481           c = 1.0;
 482           s = c;
 483           p = 0.0;
 484
 485           for (i = m - 1; i >= l; i--)
 486           {
 487             f = s * e[i - 1];
 488             b = c * e[i - 1];
 489
 490             if (Math.abs(f) >= Math.abs(g))
 491             {
 492               c = g / f;
 493               r = Math.sqrt((c * c) + 1.0);
 494               e[i] = f * r;
 495               s = 1.0 / r;
 496               c *= s;
 497             }
 498             else
 499             {
 500               s = f / g;
 501               r = Math.sqrt((s * s) + 1.0);
 502               e[i] = g * r;
 503               c = 1.0 / r;
 504               s *= c;
 505             }
 506
 507             g = d[i] - p;
 508             r = ((d[i - 1] - g) * s) + (2.0 * c * b);
 509             p = s * r;
 510             d[i] = g + p;
 511             g = (c * r) - b;
 512
 513             for (k = 1; k <= n; k++)
 514             {
 515               f = getValue(k - 1, i);
 516               setValue(k - 1, i, (s * getValue(k - 1, i - 1)) + (c * f));
 517               setValue(k - 1, i - 1, (c * getValue(k - 1, i - 1)) - (s * f));
 518             }
 519           }
 520
 521           d[l - 1] = d[l - 1] - p;
 522           e[l - 1] = g;
 523           e[m - 1] = 0.0;
 524         }
 525       } while (m != l);
 526     }
 527   }
 528
 529   @Override
 530   public double getValue(int i, int j)
 531   {
 532     return value[i][j];
 533   }
 534
 535   public void setValue(int i, int j, double val)
 536   {
 537     value[i][j] = val;
 538   }
 539
 540   /**
 541    * DOCUMENT ME!
 542    */
 543   public void tred2()
 544   {
 545     int n = rows;
 546     int l;
 547     int k;
 548     int j;
 549     int i;
 550
 551     double scale;
 552     double hh;
 553     double h;
 554     double g;
 555     double f;
 556
 557     this.d = new double[rows];
 558     this.e = new double[rows];
 559
 560     for (i = n - 1; i >= 1; i--)
 561     {
 562       l = i - 1;
 563       h = 0.0;
 564       scale = 0.0;
 565
 566       if (l > 0)
 567       {
 568         for (k = 0; k < l; k++)
 569         {
 570           scale += Math.abs(value[i][k]);
 571         }
 572
 573         if (scale == 0.0)
 574         {
 575           e[i] = value[i][l];
 576         }
 577         else
 578         {
 579           for (k = 0; k < l; k++)
 580           {
 581             value[i][k] /= scale;
 582             h += (value[i][k] * value[i][k]);
 583           }
 584
 585           f = value[i][l];
 586
 587           if (f > 0)
 588           {
 589             g = -1.0 * Math.sqrt(h);
 590           }
 591           else
 592           {
 593             g = Math.sqrt(h);
 594           }
 595
 596           e[i] = scale * g;
 597           h -= (f * g);
 598           value[i][l] = f - g;
 599           f = 0.0;
 600
 601           for (j = 0; j < l; j++)
 602           {
 603             value[j][i] = value[i][j] / h;
 604             g = 0.0;
 605
 606             for (k = 0; k < j; k++)
 607             {
 608               g += (value[j][k] * value[i][k]);
 609             }
 610
 611             for (k = j; k < l; k++)
 612             {
 613               g += (value[k][j] * value[i][k]);
 614             }
 615
 616             e[j] = g / h;
 617             f += (e[j] * value[i][j]);
 618           }
 619
 620           hh = f / (h + h);
 621
 622           for (j = 0; j < l; j++)
 623           {
 624             f = value[i][j];
 625             g = e[j] - (hh * f);
 626             e[j] = g;
 627
 628             for (k = 0; k < j; k++)
 629             {
 630               value[j][k] -= ((f * e[k]) + (g * value[i][k]));
 631             }
 632           }
 633         }
 634       }
 635       else
 636       {
 637         e[i] = value[i][l];
 638       }
 639
 640       d[i] = h;
 641     }
 642
 643     d[0] = 0.0;
 644     e[0] = 0.0;
 645
 646     for (i = 0; i < n; i++)
 647     {
 648       l = i - 1;
 649
 650       if (d[i] != 0.0)
 651       {
 652         for (j = 0; j < l; j++)
 653         {
 654           g = 0.0;
 655
 656           for (k = 0; k < l; k++)
 657           {
 658             g += (value[i][k] * value[k][j]);
 659           }
 660
 661           for (k = 0; k < l; k++)
 662           {
 663             value[k][j] -= (g * value[k][i]);
 664           }
 665         }
 666       }
 667
 668       d[i] = value[i][i];
 669       value[i][i] = 1.0;
 670
 671       for (j = 0; j < l; j++)
 672       {
 673         value[j][i] = 0.0;
 674         value[i][j] = 0.0;
 675       }
 676     }
 677   }
 678
 679   /**
 680    * DOCUMENT ME!
 681    */
 682   public void tqli2() throws Exception
 683   {
 684     int n = rows;
 685
 686     int m;
 687     int l;
 688     int iter;
 689     int i;
 690     int k;
 691     double s;
 692     double r;
 693     double p;
 694     ;
 695
 696     double g;
 697     double f;
 698     double dd;
 699     double c;
 700     double b;
 701
 702     for (i = 2; i <= n; i++)
 703     {
 704       e[i - 2] = e[i - 1];
 705     }
 706
 707     e[n - 1] = 0.0;
 708
 709     for (l = 1; l <= n; l++)
 710     {
 711       iter = 0;
 712
 713       do
 714       {
 715         for (m = l; m <= (n - 1); m++)
 716         {
 717           dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
 718
 719           if ((Math.abs(e[m - 1]) + dd) == dd)
 720           {
 721             break;
 722           }
 723         }
 724
 725         if (m != l)
 726         {
 727           iter++;
 728
 729           if (iter == maxIter)
 730           {
 731             throw new Exception(MessageManager.formatMessage(
 732                     "exception.matrix_too_many_iteration", new String[] {
 733                         "tqli2", Integer.valueOf(maxIter).toString() }));
 734           }
 735           else
 736           {
 737             // System.out.println("Iteration " + iter);
 738           }
 739
 740           g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
 741           r = Math.sqrt((g * g) + 1.0);
 742           g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
 743           c = 1.0;
 744           s = c;
 745           p = 0.0;
 746
 747           for (i = m - 1; i >= l; i--)
 748           {
 749             f = s * e[i - 1];
 750             b = c * e[i - 1];
 751
 752             if (Math.abs(f) >= Math.abs(g))
 753             {
 754               c = g / f;
 755               r = Math.sqrt((c * c) + 1.0);
 756               e[i] = f * r;
 757               s = 1.0 / r;
 758               c *= s;
 759             }
 760             else
 761             {
 762               s = f / g;
 763               r = Math.sqrt((s * s) + 1.0);
 764               e[i] = g * r;
 765               c = 1.0 / r;
 766               s *= c;
 767             }
 768
 769             g = d[i] - p;
 770             r = ((d[i - 1] - g) * s) + (2.0 * c * b);
 771             p = s * r;
 772             d[i] = g + p;
 773             g = (c * r) - b;
 774
 775             for (k = 1; k <= n; k++)
 776             {
 777               f = value[k - 1][i];
 778               value[k - 1][i] = (s * value[k - 1][i - 1]) + (c * f);
 779               value[k - 1][i - 1] = (c * value[k - 1][i - 1]) - (s * f);
 780             }
 781           }
 782
 783           d[l - 1] = d[l - 1] - p;
 784           e[l - 1] = g;
 785           e[m - 1] = 0.0;
 786         }
 787       } while (m != l);
 788     }
 789   }
 790
 791   /**
 792    * Answers the first argument with the sign of the second argument
 793    *
 794    * @param a
 795    * @param b
 796    *
 797    * @return
 798    */
 799   static double sign(double a, double b)
 800   {
 801     if (b < 0)
 802     {
 803       return -Math.abs(a);
 804     }
 805     else
 806     {
 807       return Math.abs(a);
 808     }
 809   }
 810
 811   /**
 812    * Returns an array containing the values in the specified column
 813    *
 814    * @param col
 815    *
 816    * @return
 817    */
 818   public double[] getColumn(int col)
 819   {
 820     double[] out = new double[rows];
 821
 822     for (int i = 0; i < rows; i++)
 823     {
 824       out[i] = value[i][col];
 825     }
 826
 827     return out;
 828   }
 829
 830   /**
 831    * DOCUMENT ME!
 832    *
 833    * @param ps
 834    *          DOCUMENT ME!
 835    * @param format
 836    */
 837   @Override
 838   public void printD(PrintStream ps, String format)
 839   {
 840     for (int j = 0; j < rows; j++)
 841     {
 842       Format.print(ps, format, d[j]);
 843     }
 844   }
 845
 846   /**
 847    * DOCUMENT ME!
 848    *
 849    * @param ps
 850    *          DOCUMENT ME!
 851    * @param format TODO
 852    */
 853   @Override
 854   public void printE(PrintStream ps, String format)
 855   {
 856     for (int j = 0; j < rows; j++)
 857     {
 858       Format.print(ps, format, e[j]);
 859     }
 860   }
 861
 862   @Override
 863   public double[] getD()
 864   {
 865     return d;
 866   }
 867
 868   @Override
 869   public double[] getE()
 870   {
 871     return e;
 872   }
 873
 874   @Override
 875   public int height() {
 876     return rows;
 877   }
 878
 879   @Override
 880   public int width()
 881   {
 882     return cols;
 883   }
 884
 885   @Override
 886   public double[] getRow(int i)
 887   {
 888     double[] row = new double[cols];
 889     System.arraycopy(value[i], 0, row, 0, cols);
 890     return row;
 891   }
 892 }