src/jalview/math/Matrix.java

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