src/jalview/math/Matrix.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.9)
   3  * Copyright (C) 2015 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.
  61    *
  62    * @param value
  63    *          DOCUMENT ME!
  64    * @param rows
  65    *          DOCUMENT ME!
  66    * @param cols
  67    *          DOCUMENT ME!
  68    */
  69   public Matrix(double[][] value, int rows, int cols)
  70   {
  71     this.rows = rows;
  72     this.cols = cols;
  73     this.value = value;
  74   }
  75
  76   /**
  77    * DOCUMENT ME!
  78    *
  79    * @return DOCUMENT ME!
  80    */
  81   public Matrix transpose()
  82   {
  83     double[][] out = new double[cols][rows];
  84
  85     for (int i = 0; i < cols; i++)
  86     {
  87       for (int j = 0; j < rows; j++)
  88       {
  89         out[i][j] = value[j][i];
  90       }
  91     }
  92
  93     return new Matrix(out, cols, rows);
  94   }
  95
  96   /**
  97    * DOCUMENT ME!
  98    *
  99    * @param ps
 100    *          DOCUMENT ME!
 101    */
 102   public void print(PrintStream ps)
 103   {
 104     for (int i = 0; i < rows; i++)
 105     {
 106       for (int j = 0; j < cols; j++)
 107       {
 108         Format.print(ps, "%8.2f", value[i][j]);
 109       }
 110
 111       ps.println();
 112     }
 113   }
 114
 115   /**
 116    * DOCUMENT ME!
 117    *
 118    * @param in
 119    *          DOCUMENT ME!
 120    *
 121    * @return DOCUMENT ME!
 122    */
 123   public Matrix preMultiply(Matrix in)
 124   {
 125     double[][] tmp = new double[in.rows][this.cols];
 126
 127     for (int i = 0; i < in.rows; i++)
 128     {
 129       for (int j = 0; j < this.cols; j++)
 130       {
 131         tmp[i][j] = 0.0;
 132
 133         for (int k = 0; k < in.cols; k++)
 134         {
 135           tmp[i][j] += (in.value[i][k] * this.value[k][j]);
 136         }
 137       }
 138     }
 139
 140     return new Matrix(tmp, in.rows, this.cols);
 141   }
 142
 143   /**
 144    * DOCUMENT ME!
 145    *
 146    * @param in
 147    *          DOCUMENT ME!
 148    *
 149    * @return DOCUMENT ME!
 150    */
 151   public double[] vectorPostMultiply(double[] in)
 152   {
 153     double[] out = new double[in.length];
 154
 155     for (int i = 0; i < in.length; i++)
 156     {
 157       out[i] = 0.0;
 158
 159       for (int k = 0; k < in.length; k++)
 160       {
 161         out[i] += (value[i][k] * in[k]);
 162       }
 163     }
 164
 165     return out;
 166   }
 167
 168   /**
 169    * DOCUMENT ME!
 170    *
 171    * @param in
 172    *          DOCUMENT ME!
 173    *
 174    * @return DOCUMENT ME!
 175    */
 176   public Matrix postMultiply(Matrix in)
 177   {
 178     double[][] out = new double[this.rows][in.cols];
 179
 180     for (int i = 0; i < this.rows; i++)
 181     {
 182       for (int j = 0; j < in.cols; j++)
 183       {
 184         out[i][j] = 0.0;
 185
 186         for (int k = 0; k < rows; k++)
 187         {
 188           out[i][j] = out[i][j] + (value[i][k] * in.value[k][j]);
 189         }
 190       }
 191     }
 192
 193     return new Matrix(out, this.cols, in.rows);
 194   }
 195
 196   /**
 197    * DOCUMENT ME!
 198    *
 199    * @return DOCUMENT ME!
 200    */
 201   public Matrix copy()
 202   {
 203     double[][] newmat = new double[rows][cols];
 204
 205     for (int i = 0; i < rows; i++)
 206     {
 207       for (int j = 0; j < cols; j++)
 208       {
 209         newmat[i][j] = value[i][j];
 210       }
 211     }
 212
 213     return new Matrix(newmat, rows, cols);
 214   }
 215
 216   /**
 217    * DOCUMENT ME!
 218    */
 219   public void tred()
 220   {
 221     int n = rows;
 222     int l;
 223     int k;
 224     int j;
 225     int i;
 226
 227     double scale;
 228     double hh;
 229     double h;
 230     double g;
 231     double f;
 232
 233     this.d = new double[rows];
 234     this.e = new double[rows];
 235
 236     for (i = n; i >= 2; i--)
 237     {
 238       l = i - 1;
 239       h = 0.0;
 240       scale = 0.0;
 241
 242       if (l > 1)
 243       {
 244         for (k = 1; k <= l; k++)
 245         {
 246           scale += Math.abs(value[i - 1][k - 1]);
 247         }
 248
 249         if (scale == 0.0)
 250         {
 251           e[i - 1] = value[i - 1][l - 1];
 252         }
 253         else
 254         {
 255           for (k = 1; k <= l; k++)
 256           {
 257             value[i - 1][k - 1] /= scale;
 258             h += (value[i - 1][k - 1] * value[i - 1][k - 1]);
 259           }
 260
 261           f = value[i - 1][l - 1];
 262
 263           if (f > 0)
 264           {
 265             g = -1.0 * Math.sqrt(h);
 266           }
 267           else
 268           {
 269             g = Math.sqrt(h);
 270           }
 271
 272           e[i - 1] = scale * g;
 273           h -= (f * g);
 274           value[i - 1][l - 1] = f - g;
 275           f = 0.0;
 276
 277           for (j = 1; j <= l; j++)
 278           {
 279             value[j - 1][i - 1] = value[i - 1][j - 1] / h;
 280             g = 0.0;
 281
 282             for (k = 1; k <= j; k++)
 283             {
 284               g += (value[j - 1][k - 1] * value[i - 1][k - 1]);
 285             }
 286
 287             for (k = j + 1; k <= l; k++)
 288             {
 289               g += (value[k - 1][j - 1] * value[i - 1][k - 1]);
 290             }
 291
 292             e[j - 1] = g / h;
 293             f += (e[j - 1] * value[i - 1][j - 1]);
 294           }
 295
 296           hh = f / (h + h);
 297
 298           for (j = 1; j <= l; j++)
 299           {
 300             f = value[i - 1][j - 1];
 301             g = e[j - 1] - (hh * f);
 302             e[j - 1] = g;
 303
 304             for (k = 1; k <= j; k++)
 305             {
 306               value[j - 1][k - 1] -= ((f * e[k - 1]) + (g * value[i - 1][k - 1]));
 307             }
 308           }
 309         }
 310       }
 311       else
 312       {
 313         e[i - 1] = value[i - 1][l - 1];
 314       }
 315
 316       d[i - 1] = h;
 317     }
 318
 319     d[0] = 0.0;
 320     e[0] = 0.0;
 321
 322     for (i = 1; i <= n; i++)
 323     {
 324       l = i - 1;
 325
 326       if (d[i - 1] != 0.0)
 327       {
 328         for (j = 1; j <= l; j++)
 329         {
 330           g = 0.0;
 331
 332           for (k = 1; k <= l; k++)
 333           {
 334             g += (value[i - 1][k - 1] * value[k - 1][j - 1]);
 335           }
 336
 337           for (k = 1; k <= l; k++)
 338           {
 339             value[k - 1][j - 1] -= (g * value[k - 1][i - 1]);
 340           }
 341         }
 342       }
 343
 344       d[i - 1] = value[i - 1][i - 1];
 345       value[i - 1][i - 1] = 1.0;
 346
 347       for (j = 1; j <= l; j++)
 348       {
 349         value[j - 1][i - 1] = 0.0;
 350         value[i - 1][j - 1] = 0.0;
 351       }
 352     }
 353   }
 354
 355   /**
 356    * DOCUMENT ME!
 357    */
 358   public void tqli() throws Exception
 359   {
 360     int n = rows;
 361
 362     int m;
 363     int l;
 364     int iter;
 365     int i;
 366     int k;
 367     double s;
 368     double r;
 369     double p;
 370     ;
 371
 372     double g;
 373     double f;
 374     double dd;
 375     double c;
 376     double b;
 377
 378     for (i = 2; i <= n; i++)
 379     {
 380       e[i - 2] = e[i - 1];
 381     }
 382
 383     e[n - 1] = 0.0;
 384
 385     for (l = 1; l <= n; l++)
 386     {
 387       iter = 0;
 388
 389       do
 390       {
 391         for (m = l; m <= (n - 1); m++)
 392         {
 393           dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
 394
 395           if ((Math.abs(e[m - 1]) + dd) == dd)
 396           {
 397             break;
 398           }
 399         }
 400
 401         if (m != l)
 402         {
 403           iter++;
 404
 405           if (iter == maxIter)
 406           {
 407             throw new Exception(MessageManager.formatMessage(
 408                     "exception.matrix_too_many_iteration", new String[] {
 409                         "tqli", Integer.valueOf(maxIter).toString() }));
 410           }
 411           else
 412           {
 413             // System.out.println("Iteration " + iter);
 414           }
 415
 416           g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
 417           r = Math.sqrt((g * g) + 1.0);
 418           g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
 419           c = 1.0;
 420           s = c;
 421           p = 0.0;
 422
 423           for (i = m - 1; i >= l; i--)
 424           {
 425             f = s * e[i - 1];
 426             b = c * e[i - 1];
 427
 428             if (Math.abs(f) >= Math.abs(g))
 429             {
 430               c = g / f;
 431               r = Math.sqrt((c * c) + 1.0);
 432               e[i] = f * r;
 433               s = 1.0 / r;
 434               c *= s;
 435             }
 436             else
 437             {
 438               s = f / g;
 439               r = Math.sqrt((s * s) + 1.0);
 440               e[i] = g * r;
 441               c = 1.0 / r;
 442               s *= c;
 443             }
 444
 445             g = d[i] - p;
 446             r = ((d[i - 1] - g) * s) + (2.0 * c * b);
 447             p = s * r;
 448             d[i] = g + p;
 449             g = (c * r) - b;
 450
 451             for (k = 1; k <= n; k++)
 452             {
 453               f = value[k - 1][i];
 454               value[k - 1][i] = (s * value[k - 1][i - 1]) + (c * f);
 455               value[k - 1][i - 1] = (c * value[k - 1][i - 1]) - (s * f);
 456             }
 457           }
 458
 459           d[l - 1] = d[l - 1] - p;
 460           e[l - 1] = g;
 461           e[m - 1] = 0.0;
 462         }
 463       } while (m != l);
 464     }
 465   }
 466
 467   /**
 468    * DOCUMENT ME!
 469    */
 470   public void tred2()
 471   {
 472     int n = rows;
 473     int l;
 474     int k;
 475     int j;
 476     int i;
 477
 478     double scale;
 479     double hh;
 480     double h;
 481     double g;
 482     double f;
 483
 484     this.d = new double[rows];
 485     this.e = new double[rows];
 486
 487     for (i = n - 1; i >= 1; i--)
 488     {
 489       l = i - 1;
 490       h = 0.0;
 491       scale = 0.0;
 492
 493       if (l > 0)
 494       {
 495         for (k = 0; k < l; k++)
 496         {
 497           scale += Math.abs(value[i][k]);
 498         }
 499
 500         if (scale == 0.0)
 501         {
 502           e[i] = value[i][l];
 503         }
 504         else
 505         {
 506           for (k = 0; k < l; k++)
 507           {
 508             value[i][k] /= scale;
 509             h += (value[i][k] * value[i][k]);
 510           }
 511
 512           f = value[i][l];
 513
 514           if (f > 0)
 515           {
 516             g = -1.0 * Math.sqrt(h);
 517           }
 518           else
 519           {
 520             g = Math.sqrt(h);
 521           }
 522
 523           e[i] = scale * g;
 524           h -= (f * g);
 525           value[i][l] = f - g;
 526           f = 0.0;
 527
 528           for (j = 0; j < l; j++)
 529           {
 530             value[j][i] = value[i][j] / h;
 531             g = 0.0;
 532
 533             for (k = 0; k < j; k++)
 534             {
 535               g += (value[j][k] * value[i][k]);
 536             }
 537
 538             for (k = j; k < l; k++)
 539             {
 540               g += (value[k][j] * value[i][k]);
 541             }
 542
 543             e[j] = g / h;
 544             f += (e[j] * value[i][j]);
 545           }
 546
 547           hh = f / (h + h);
 548
 549           for (j = 0; j < l; j++)
 550           {
 551             f = value[i][j];
 552             g = e[j] - (hh * f);
 553             e[j] = g;
 554
 555             for (k = 0; k < j; k++)
 556             {
 557               value[j][k] -= ((f * e[k]) + (g * value[i][k]));
 558             }
 559           }
 560         }
 561       }
 562       else
 563       {
 564         e[i] = value[i][l];
 565       }
 566
 567       d[i] = h;
 568     }
 569
 570     d[0] = 0.0;
 571     e[0] = 0.0;
 572
 573     for (i = 0; i < n; i++)
 574     {
 575       l = i - 1;
 576
 577       if (d[i] != 0.0)
 578       {
 579         for (j = 0; j < l; j++)
 580         {
 581           g = 0.0;
 582
 583           for (k = 0; k < l; k++)
 584           {
 585             g += (value[i][k] * value[k][j]);
 586           }
 587
 588           for (k = 0; k < l; k++)
 589           {
 590             value[k][j] -= (g * value[k][i]);
 591           }
 592         }
 593       }
 594
 595       d[i] = value[i][i];
 596       value[i][i] = 1.0;
 597
 598       for (j = 0; j < l; j++)
 599       {
 600         value[j][i] = 0.0;
 601         value[i][j] = 0.0;
 602       }
 603     }
 604   }
 605
 606   /**
 607    * DOCUMENT ME!
 608    */
 609   public void tqli2() throws Exception
 610   {
 611     int n = rows;
 612
 613     int m;
 614     int l;
 615     int iter;
 616     int i;
 617     int k;
 618     double s;
 619     double r;
 620     double p;
 621     ;
 622
 623     double g;
 624     double f;
 625     double dd;
 626     double c;
 627     double b;
 628
 629     for (i = 2; i <= n; i++)
 630     {
 631       e[i - 2] = e[i - 1];
 632     }
 633
 634     e[n - 1] = 0.0;
 635
 636     for (l = 1; l <= n; l++)
 637     {
 638       iter = 0;
 639
 640       do
 641       {
 642         for (m = l; m <= (n - 1); m++)
 643         {
 644           dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
 645
 646           if ((Math.abs(e[m - 1]) + dd) == dd)
 647           {
 648             break;
 649           }
 650         }
 651
 652         if (m != l)
 653         {
 654           iter++;
 655
 656           if (iter == maxIter)
 657           {
 658             throw new Exception(MessageManager.formatMessage(
 659                     "exception.matrix_too_many_iteration", new String[] {
 660                         "tqli2", Integer.valueOf(maxIter).toString() }));
 661           }
 662           else
 663           {
 664             // System.out.println("Iteration " + iter);
 665           }
 666
 667           g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
 668           r = Math.sqrt((g * g) + 1.0);
 669           g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
 670           c = 1.0;
 671           s = c;
 672           p = 0.0;
 673
 674           for (i = m - 1; i >= l; i--)
 675           {
 676             f = s * e[i - 1];
 677             b = c * e[i - 1];
 678
 679             if (Math.abs(f) >= Math.abs(g))
 680             {
 681               c = g / f;
 682               r = Math.sqrt((c * c) + 1.0);
 683               e[i] = f * r;
 684               s = 1.0 / r;
 685               c *= s;
 686             }
 687             else
 688             {
 689               s = f / g;
 690               r = Math.sqrt((s * s) + 1.0);
 691               e[i] = g * r;
 692               c = 1.0 / r;
 693               s *= c;
 694             }
 695
 696             g = d[i] - p;
 697             r = ((d[i - 1] - g) * s) + (2.0 * c * b);
 698             p = s * r;
 699             d[i] = g + p;
 700             g = (c * r) - b;
 701
 702             for (k = 1; k <= n; k++)
 703             {
 704               f = value[k - 1][i];
 705               value[k - 1][i] = (s * value[k - 1][i - 1]) + (c * f);
 706               value[k - 1][i - 1] = (c * value[k - 1][i - 1]) - (s * f);
 707             }
 708           }
 709
 710           d[l - 1] = d[l - 1] - p;
 711           e[l - 1] = g;
 712           e[m - 1] = 0.0;
 713         }
 714       } while (m != l);
 715     }
 716   }
 717
 718   /**
 719    * DOCUMENT ME!
 720    *
 721    * @param a
 722    *          DOCUMENT ME!
 723    * @param b
 724    *          DOCUMENT ME!
 725    *
 726    * @return DOCUMENT ME!
 727    */
 728   public double sign(double a, double b)
 729   {
 730     if (b < 0)
 731     {
 732       return -Math.abs(a);
 733     }
 734     else
 735     {
 736       return Math.abs(a);
 737     }
 738   }
 739
 740   /**
 741    * DOCUMENT ME!
 742    *
 743    * @param n
 744    *          DOCUMENT ME!
 745    *
 746    * @return DOCUMENT ME!
 747    */
 748   public double[] getColumn(int n)
 749   {
 750     double[] out = new double[rows];
 751
 752     for (int i = 0; i < rows; i++)
 753     {
 754       out[i] = value[i][n];
 755     }
 756
 757     return out;
 758   }
 759
 760   /**
 761    * DOCUMENT ME!
 762    *
 763    * @param ps
 764    *          DOCUMENT ME!
 765    */
 766   public void printD(PrintStream ps)
 767   {
 768     for (int j = 0; j < rows; j++)
 769     {
 770       Format.print(ps, "%15.4e", d[j]);
 771     }
 772   }
 773
 774   /**
 775    * DOCUMENT ME!
 776    *
 777    * @param ps
 778    *          DOCUMENT ME!
 779    */
 780   public void printE(PrintStream ps)
 781   {
 782     for (int j = 0; j < rows; j++)
 783     {
 784       Format.print(ps, "%15.4e", e[j]);
 785     }
 786   }
 787
 788   /**
 789    * DOCUMENT ME!
 790    *
 791    * @param args
 792    *          DOCUMENT ME!
 793    */
 794   public static void main(String[] args) throws Exception
 795   {
 796     int n = Integer.parseInt(args[0]);
 797     double[][] in = new double[n][n];
 798
 799     for (int i = 0; i < n; i++)
 800     {
 801       for (int j = 0; j < n; j++)
 802       {
 803         in[i][j] = (double) Math.random();
 804       }
 805     }
 806
 807     Matrix origmat = new Matrix(in, n, n);
 808
 809     // System.out.println(" --- Original matrix ---- ");
 810     // / origmat.print(System.out);
 811     // System.out.println();
 812     // System.out.println(" --- transpose matrix ---- ");
 813     Matrix trans = origmat.transpose();
 814
 815     // trans.print(System.out);
 816     // System.out.println();
 817     // System.out.println(" --- OrigT * Orig ---- ");
 818     Matrix symm = trans.postMultiply(origmat);
 819
 820     // symm.print(System.out);
 821     // System.out.println();
 822     // Copy the symmetric matrix for later
 823     // Matrix origsymm = symm.copy();
 824
 825     // This produces the tridiagonal transformation matrix
 826     // long tstart = System.currentTimeMillis();
 827     symm.tred();
 828
 829     // long tend = System.currentTimeMillis();
 830
 831     // System.out.println("Time take for tred = " + (tend-tstart) + "ms");
 832     // System.out.println(" ---Tridiag transform matrix ---");
 833     // symm.print(System.out);
 834     // System.out.println();
 835     // System.out.println(" --- D vector ---");
 836     // symm.printD(System.out);
 837     // System.out.println();
 838     // System.out.println(" --- E vector ---");
 839     // symm.printE(System.out);
 840     // System.out.println();
 841     // Now produce the diagonalization matrix
 842     // tstart = System.currentTimeMillis();
 843     symm.tqli();
 844     // tend = System.currentTimeMillis();
 845
 846     // System.out.println("Time take for tqli = " + (tend-tstart) + " ms");
 847     // System.out.println(" --- New diagonalization matrix ---");
 848     // symm.print(System.out);
 849     // System.out.println();
 850     // System.out.println(" --- D vector ---");
 851     // symm.printD(System.out);
 852     // System.out.println();
 853     // System.out.println(" --- E vector ---");
 854     // symm.printE(System.out);
 855     // System.out.println();
 856     // System.out.println(" --- First eigenvector --- ");
 857     // double[] eigenv = symm.getColumn(0);
 858     // for (int i=0; i < eigenv.length;i++) {
 859     // Format.print(System.out,"%15.4f",eigenv[i]);
 860     // }
 861     // System.out.println();
 862     // double[] neigenv = origsymm.vectorPostMultiply(eigenv);
 863     // for (int i=0; i < neigenv.length;i++) {
 864     // Format.print(System.out,"%15.4f",neigenv[i]/symm.d[0]);
 865     // }
 866     // System.out.println();
 867   }
 868 }