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   @Override
 536   public void setValue(int i, int j, double val)
 537   {
 538     value[i][j] = val;
 539   }
 540
 541   /**
 542    * DOCUMENT ME!
 543    */
 544   public void tred2()
 545   {
 546     int n = rows;
 547     int l;
 548     int k;
 549     int j;
 550     int i;
 551
 552     double scale;
 553     double hh;
 554     double h;
 555     double g;
 556     double f;
 557
 558     this.d = new double[rows];
 559     this.e = new double[rows];
 560
 561     for (i = n - 1; i >= 1; i--)
 562     {
 563       l = i - 1;
 564       h = 0.0;
 565       scale = 0.0;
 566
 567       if (l > 0)
 568       {
 569         for (k = 0; k < l; k++)
 570         {
 571           scale += Math.abs(value[i][k]);
 572         }
 573
 574         if (scale == 0.0)
 575         {
 576           e[i] = value[i][l];
 577         }
 578         else
 579         {
 580           for (k = 0; k < l; k++)
 581           {
 582             value[i][k] /= scale;
 583             h += (value[i][k] * value[i][k]);
 584           }
 585
 586           f = value[i][l];
 587
 588           if (f > 0)
 589           {
 590             g = -1.0 * Math.sqrt(h);
 591           }
 592           else
 593           {
 594             g = Math.sqrt(h);
 595           }
 596
 597           e[i] = scale * g;
 598           h -= (f * g);
 599           value[i][l] = f - g;
 600           f = 0.0;
 601
 602           for (j = 0; j < l; j++)
 603           {
 604             value[j][i] = value[i][j] / h;
 605             g = 0.0;
 606
 607             for (k = 0; k < j; k++)
 608             {
 609               g += (value[j][k] * value[i][k]);
 610             }
 611
 612             for (k = j; k < l; k++)
 613             {
 614               g += (value[k][j] * value[i][k]);
 615             }
 616
 617             e[j] = g / h;
 618             f += (e[j] * value[i][j]);
 619           }
 620
 621           hh = f / (h + h);
 622
 623           for (j = 0; j < l; j++)
 624           {
 625             f = value[i][j];
 626             g = e[j] - (hh * f);
 627             e[j] = g;
 628
 629             for (k = 0; k < j; k++)
 630             {
 631               value[j][k] -= ((f * e[k]) + (g * value[i][k]));
 632             }
 633           }
 634         }
 635       }
 636       else
 637       {
 638         e[i] = value[i][l];
 639       }
 640
 641       d[i] = h;
 642     }
 643
 644     d[0] = 0.0;
 645     e[0] = 0.0;
 646
 647     for (i = 0; i < n; i++)
 648     {
 649       l = i - 1;
 650
 651       if (d[i] != 0.0)
 652       {
 653         for (j = 0; j < l; j++)
 654         {
 655           g = 0.0;
 656
 657           for (k = 0; k < l; k++)
 658           {
 659             g += (value[i][k] * value[k][j]);
 660           }
 661
 662           for (k = 0; k < l; k++)
 663           {
 664             value[k][j] -= (g * value[k][i]);
 665           }
 666         }
 667       }
 668
 669       d[i] = value[i][i];
 670       value[i][i] = 1.0;
 671
 672       for (j = 0; j < l; j++)
 673       {
 674         value[j][i] = 0.0;
 675         value[i][j] = 0.0;
 676       }
 677     }
 678   }
 679
 680   /**
 681    * DOCUMENT ME!
 682    */
 683   public void tqli2() throws Exception
 684   {
 685     int n = rows;
 686
 687     int m;
 688     int l;
 689     int iter;
 690     int i;
 691     int k;
 692     double s;
 693     double r;
 694     double p;
 695     ;
 696
 697     double g;
 698     double f;
 699     double dd;
 700     double c;
 701     double b;
 702
 703     for (i = 2; i <= n; i++)
 704     {
 705       e[i - 2] = e[i - 1];
 706     }
 707
 708     e[n - 1] = 0.0;
 709
 710     for (l = 1; l <= n; l++)
 711     {
 712       iter = 0;
 713
 714       do
 715       {
 716         for (m = l; m <= (n - 1); m++)
 717         {
 718           dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
 719
 720           if ((Math.abs(e[m - 1]) + dd) == dd)
 721           {
 722             break;
 723           }
 724         }
 725
 726         if (m != l)
 727         {
 728           iter++;
 729
 730           if (iter == maxIter)
 731           {
 732             throw new Exception(MessageManager.formatMessage(
 733                     "exception.matrix_too_many_iteration", new String[] {
 734                         "tqli2", Integer.valueOf(maxIter).toString() }));
 735           }
 736           else
 737           {
 738             // System.out.println("Iteration " + iter);
 739           }
 740
 741           g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
 742           r = Math.sqrt((g * g) + 1.0);
 743           g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
 744           c = 1.0;
 745           s = c;
 746           p = 0.0;
 747
 748           for (i = m - 1; i >= l; i--)
 749           {
 750             f = s * e[i - 1];
 751             b = c * e[i - 1];
 752
 753             if (Math.abs(f) >= Math.abs(g))
 754             {
 755               c = g / f;
 756               r = Math.sqrt((c * c) + 1.0);
 757               e[i] = f * r;
 758               s = 1.0 / r;
 759               c *= s;
 760             }
 761             else
 762             {
 763               s = f / g;
 764               r = Math.sqrt((s * s) + 1.0);
 765               e[i] = g * r;
 766               c = 1.0 / r;
 767               s *= c;
 768             }
 769
 770             g = d[i] - p;
 771             r = ((d[i - 1] - g) * s) + (2.0 * c * b);
 772             p = s * r;
 773             d[i] = g + p;
 774             g = (c * r) - b;
 775
 776             for (k = 1; k <= n; k++)
 777             {
 778               f = value[k - 1][i];
 779               value[k - 1][i] = (s * value[k - 1][i - 1]) + (c * f);
 780               value[k - 1][i - 1] = (c * value[k - 1][i - 1]) - (s * f);
 781             }
 782           }
 783
 784           d[l - 1] = d[l - 1] - p;
 785           e[l - 1] = g;
 786           e[m - 1] = 0.0;
 787         }
 788       } while (m != l);
 789     }
 790   }
 791
 792   /**
 793    * Answers the first argument with the sign of the second argument
 794    *
 795    * @param a
 796    * @param b
 797    *
 798    * @return
 799    */
 800   static double sign(double a, double b)
 801   {
 802     if (b < 0)
 803     {
 804       return -Math.abs(a);
 805     }
 806     else
 807     {
 808       return Math.abs(a);
 809     }
 810   }
 811
 812   /**
 813    * Returns an array containing the values in the specified column
 814    *
 815    * @param col
 816    *
 817    * @return
 818    */
 819   public double[] getColumn(int col)
 820   {
 821     double[] out = new double[rows];
 822
 823     for (int i = 0; i < rows; i++)
 824     {
 825       out[i] = value[i][col];
 826     }
 827
 828     return out;
 829   }
 830
 831   /**
 832    * DOCUMENT ME!
 833    *
 834    * @param ps
 835    *          DOCUMENT ME!
 836    * @param format
 837    */
 838   @Override
 839   public void printD(PrintStream ps, String format)
 840   {
 841     for (int j = 0; j < rows; j++)
 842     {
 843       Format.print(ps, format, d[j]);
 844     }
 845   }
 846
 847   /**
 848    * DOCUMENT ME!
 849    *
 850    * @param ps
 851    *          DOCUMENT ME!
 852    * @param format TODO
 853    */
 854   @Override
 855   public void printE(PrintStream ps, String format)
 856   {
 857     for (int j = 0; j < rows; j++)
 858     {
 859       Format.print(ps, format, e[j]);
 860     }
 861   }
 862
 863   @Override
 864   public double[] getD()
 865   {
 866     return d;
 867   }
 868
 869   @Override
 870   public double[] getE()
 871   {
 872     return e;
 873   }
 874
 875   @Override
 876   public int height() {
 877     return rows;
 878   }
 879
 880   @Override
 881   public int width()
 882   {
 883     return cols;
 884   }
 885
 886   @Override
 887   public double[] getRow(int i)
 888   {
 889     double[] row = new double[cols];
 890     System.arraycopy(value[i], 0, row, 0, cols);
 891     return row;
 892   }
 893
 894   /**
 895    * Returns a length 2 array of {minValue, maxValue} of all values in the
 896    * matrix. Returns null if the matrix is null or empty.
 897    *
 898    * @return
 899    */
 900   double[] findMinMax()
 901   {
 902     if (value == null)
 903     {
 904       return null;
 905     }
 906     double min = Double.MAX_VALUE;
 907     double max = -Double.MAX_VALUE;
 908     boolean empty = true;
 909     for (double[] row : value)
 910     {
 911       if (row != null)
 912       {
 913         for (double x : row)
 914         {
 915           empty = false;
 916           if (x > max)
 917           {
 918             max = x;
 919           }
 920           if (x < min)
 921           {
 922             min = x;
 923           }
 924         }
 925       }
 926     }
 927     return empty ? null : new double[] { min, max };
 928   }
 929
 930   /**
 931    * {@inheritDoc}
 932    */
 933   @Override
 934   public void reverseRange(boolean maxToZero)
 935   {
 936     if (value == null)
 937     {
 938       return;
 939     }
 940     double[] minMax = findMinMax();
 941     if (minMax == null)
 942     {
 943       return; // empty matrix
 944     }
 945     double subtractFrom = maxToZero ? minMax[1] : minMax[0] + minMax[1];
 946
 947     for (double[] row : value)
 948     {
 949       if (row != null)
 950       {
 951         int j = 0;
 952         for (double x : row)
 953         {
 954           row[j] = subtractFrom - x;
 955           j++;
 956         }
 957       }
 958     }
 959   }
 960
 961   /**
 962    * Multiply every entry in the matrix by the given value. This method is not
 963    * thread-safe.
 964    */
 965   @Override
 966   public void multiply(double d)
 967   {
 968     for (double[] row : value)
 969     {
 970       if (row != null)
 971       {
 972         for (int i = 0; i < row.length; i++)
 973         {
 974           row[i] *= d;
 975         }
 976       }
 977     }
 978   }
 979 }