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