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