src/jalview/math/Matrix.java

   1 /*\r
   2 * Jalview - A Sequence Alignment Editor and Viewer\r
   3 * Copyright (C) 2005 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle\r
   4 *\r
   5 * This program is free software; you can redistribute it and/or\r
   6 * modify it under the terms of the GNU General Public License\r
   7 * as published by the Free Software Foundation; either version 2\r
   8 * of the License, or (at your option) any later version.\r
   9 *\r
  10 * This program is distributed in the hope that it will be useful,\r
  11 * but WITHOUT ANY WARRANTY; without even the implied warranty of\r
  12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
  13 * GNU General Public License for more details.\r
  14 *\r
  15 * You should have received a copy of the GNU General Public License\r
  16 * along with this program; if not, write to the Free Software\r
  17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA\r
  18 */\r
  19 \r
  20 package jalview.math;\r
  21 \r
  22 import jalview.util.*;\r
  23 \r
  24 import java.io.*;\r
  25 \r
  26 public class Matrix {\r
  27 \r
  28   /**\r
  29    * SMJSPUBLIC\r
  30    */\r
  31   public double[][] value;\r
  32   public int rows;\r
  33   public int cols;\r
  34   public double[] d;  // Diagonal\r
  35   public double[] e;  // off diagonal\r
  36 \r
  37   public Matrix(double[][] value, int rows, int cols) {\r
  38     this.rows = rows;\r
  39     this.cols = cols;\r
  40     this.value = value;\r
  41   }\r
  42 \r
  43   public Matrix transpose() {\r
  44     double[][] out = new double[cols][rows];\r
  45 \r
  46     for (int i = 0; i < cols; i++) {\r
  47       for (int j = 0; j < rows ; j++) {\r
  48         out[i][j] = value[j][i];\r
  49       }\r
  50     }\r
  51     return new Matrix(out,cols,rows);\r
  52   }\r
  53 \r
  54   public void print(PrintStream ps) {\r
  55 \r
  56     for (int i = 0; i < rows; i++) {\r
  57       for (int j = 0; j < cols; j++) {\r
  58         Format.print(ps,"%8.2f",value[i][j]);\r
  59       }\r
  60       ps.println();\r
  61     }\r
  62   }\r
  63 \r
  64 \r
  65   public Matrix  preMultiply(Matrix in) {\r
  66     double[][] tmp = new double[in.rows][this.cols];\r
  67 \r
  68     for (int i = 0; i < in.rows; i++) {\r
  69       for (int j = 0; j < this.cols; j++ ) {\r
  70         tmp[i][j] = 0.0;\r
  71 \r
  72         for (int k = 0; k < in.cols; k++) {\r
  73           tmp[i][j] += in.value[i][k]*this.value[k][j];\r
  74         }\r
  75 \r
  76       }\r
  77     }\r
  78 \r
  79     return new Matrix(tmp,in.rows,this.cols);\r
  80   }\r
  81 \r
  82   public double[] vectorPostMultiply(double[] in) {\r
  83     double[] out = new double[in.length];\r
  84     for (int i = 0; i < in.length; i++) {\r
  85       out[i] = 0.0;\r
  86       for (int k=0; k < in.length; k++) {\r
  87         out[i] += value[i][k] * in[k];\r
  88       }\r
  89     }\r
  90     return out;\r
  91   }\r
  92   public Matrix postMultiply(Matrix in) {\r
  93 \r
  94     double[][] out = new double[this.rows][in.cols];\r
  95     for (int i = 0; i < this.rows; i++) {\r
  96       for (int j = 0; j < in.cols; j++ ) {\r
  97 \r
  98         out[i][j] = 0.0;\r
  99 \r
 100         for (int k = 0; k < rows; k++) {\r
 101           out[i][j] = out[i][j] + value[i][k]*in.value[k][j];\r
 102         }\r
 103 \r
 104       }\r
 105     }\r
 106     return new Matrix(out,this.cols,in.rows);\r
 107   }\r
 108 \r
 109   public Matrix copy() {\r
 110     double[][] newmat = new double[rows][cols];\r
 111 \r
 112     for (int i = 0; i < rows; i++) {\r
 113       for (int j = 0; j < cols; j++) {\r
 114         newmat[i][j] = value[i][j];\r
 115       }\r
 116     }\r
 117     return new Matrix(newmat,rows,cols);\r
 118   }\r
 119 \r
 120   public void tred() {\r
 121     int n = rows;\r
 122     int l;\r
 123     int k;\r
 124     int j;\r
 125     int i;\r
 126 \r
 127     double scale;\r
 128     double hh;\r
 129     double h;\r
 130     double g;\r
 131     double f;\r
 132 \r
 133     this.d = new double[rows];\r
 134     this.e = new double[rows];\r
 135 \r
 136     for (i=n; i >= 2;i--) {\r
 137       l=i-1;\r
 138       h = 0.0;\r
 139       scale = 0.0;\r
 140 \r
 141       if (l > 1) {\r
 142         for (k=1;k<=l;k++) {\r
 143           scale += Math.abs(value[i-1][k-1]);\r
 144         }\r
 145         if (scale == 0.0) {\r
 146           e[i-1] = value[i-1][l-1];\r
 147         } else {\r
 148           for (k=1; k <= l; k++) {\r
 149             value[i-1][k-1] /= scale;\r
 150             h += value[i-1][k-1]*value[i-1][k-1];\r
 151           }\r
 152           f = value[i-1][l-1];\r
 153           if (f>0) {\r
 154             g = -1.0*Math.sqrt(h);\r
 155           } else {\r
 156             g = Math.sqrt(h);\r
 157           }\r
 158           e[i-1] = scale*g;\r
 159           h -= f*g;\r
 160           value[i-1][l-1] = f-g;\r
 161           f=0.0;\r
 162           for (j=1; j <= l; j++) {\r
 163             value[j-1][i-1] = value[i-1][j-1]/h;\r
 164             g=0.0;\r
 165             for (k= 1; k <= j; k++) {\r
 166               g += value[j-1][k-1]*value[i-1][k-1];\r
 167             }\r
 168             for (k=j+1; k<=l;k++) {\r
 169               g+= value[k-1][j-1]*value[i-1][k-1];\r
 170             }\r
 171             e[j-1] = g/h;\r
 172             f+=e[j-1]*value[i-1][j-1];\r
 173           }\r
 174           hh=f/(h+h);\r
 175           for (j=1;j<=l;j++) {\r
 176             f=value[i-1][j-1];\r
 177             g = e[j-1] - hh*f;\r
 178             e[j-1] = g;\r
 179             for (k=1;k<=j;k++) {\r
 180               value[j-1][k-1] -= (f*e[k-1]+g*value[i-1][k-1]);\r
 181             }\r
 182           }\r
 183         }\r
 184       } else {\r
 185         e[i-1] = value[i-1][l-1];\r
 186       }\r
 187       d[i-1] = h;\r
 188     }\r
 189     d[0] = 0.0;\r
 190     e[0] = 0.0;\r
 191     for (i=1;i<=n;i++) {\r
 192       l=i-1;\r
 193       if (d[i-1] != 0.0) {\r
 194         for (j=1;j<=l;j++) {\r
 195           g=0.0;\r
 196           for (k=1;k<=l;k++) {\r
 197             g+= value[i-1][k-1]*value[k-1][j-1];\r
 198           }\r
 199           for (k=1;k<=l;k++) {\r
 200             value[k-1][j-1] -= g*value[k-1][i-1];\r
 201           }\r
 202         }\r
 203       }\r
 204       d[i-1] = value[i-1][i-1];\r
 205       value[i-1][i-1] = 1.0;\r
 206       for (j=1;j<=l;j++) {\r
 207         value[j-1][i-1] = 0.0;\r
 208         value[i-1][j-1] = 0.0;\r
 209       }\r
 210     }\r
 211   }\r
 212 \r
 213   public void tqli() {\r
 214     int n = rows;\r
 215 \r
 216     int m;\r
 217     int l;\r
 218     int iter;\r
 219     int i;\r
 220     int k;\r
 221     double s;\r
 222     double r;\r
 223     double p;\r
 224     ;\r
 225     double g;\r
 226     double f;\r
 227     double dd;\r
 228     double c;\r
 229     double b;\r
 230 \r
 231     for (i=2;i<=n;i++) {\r
 232       e[i-2] = e[i-1];\r
 233     }\r
 234     e[n-1] = 0.0;\r
 235     for (l=1;l<=n;l++) {\r
 236       iter=0;\r
 237       do {\r
 238         for (m=l;m<=(n-1);m++) {\r
 239           dd=Math.abs(d[m-1]) + Math.abs(d[m]);\r
 240           if (Math.abs(e[m-1]) + dd == dd)\r
 241             break;\r
 242         }\r
 243         if (m != l) {\r
 244           iter++;\r
 245           if (iter == 30) {\r
 246             System.err.print("Too many iterations in tqli");\r
 247             System.exit(0); // JBPNote - should this really be here ???\r
 248           } else {\r
 249             //      System.out.println("Iteration " + iter);\r
 250           }\r
 251           g=(d[l]-d[l-1])/(2.0*e[l-1]);\r
 252           r = Math.sqrt((g*g) + 1.0);\r
 253           g=d[m-1]-d[l-1]+e[l-1]/(g + sign(r,g));\r
 254           c = 1.0;\r
 255           s = c;\r
 256           p=0.0;\r
 257           for (i=m-1;i>=l;i--) {\r
 258             f = s*e[i-1];\r
 259             b = c*e[i-1];\r
 260             if (Math.abs(f) >= Math.abs(g)) {\r
 261               c=g/f;\r
 262               r = Math.sqrt((c*c)+1.0);\r
 263               e[i] = f*r;\r
 264               s = 1.0/r;\r
 265               c *= s;\r
 266             } else {\r
 267               s=f/g;\r
 268               r = Math.sqrt((s*s)+1.0);\r
 269               e[i] = g*r;\r
 270               c = 1.0/r;\r
 271               s *= c;\r
 272             }\r
 273             g=d[i] -p;\r
 274             r=(d[i-1]-g)*s + 2.0*c*b;\r
 275             p=s*r;\r
 276             d[i] = g + p;\r
 277             g = c * r - b;\r
 278             for (k=1; k <= n; k++) {\r
 279               f=value[k-1][i];\r
 280               value[k-1][i] = s*value[k-1][i-1] + c*f;\r
 281               value[k-1][i-1] = c*value[k-1][i-1] - s*f;\r
 282             }\r
 283           }\r
 284           d[l-1] = d[l-1] - p;\r
 285           e[l-1] = g;\r
 286           e[m-1] = 0.0;\r
 287         }\r
 288       } while ( m != l);\r
 289     }\r
 290   }\r
 291   public void tred2() {\r
 292     int n = rows;\r
 293     int l;\r
 294     int k;\r
 295     int j;\r
 296     int i;\r
 297 \r
 298     double scale;\r
 299     double hh;\r
 300     double h;\r
 301     double g;\r
 302     double f;\r
 303 \r
 304     this.d = new double[rows];\r
 305     this.e = new double[rows];\r
 306 \r
 307     for (i=n-1; i >= 1;i--) {\r
 308       l=i-1;\r
 309       h = 0.0;\r
 310       scale = 0.0;\r
 311 \r
 312       if (l > 0) {\r
 313         for (k=0;k<l;k++) {\r
 314           scale += Math.abs(value[i][k]);\r
 315         }\r
 316         if (scale == 0.0) {\r
 317           e[i] = value[i][l];\r
 318         } else {\r
 319           for (k=0; k < l; k++) {\r
 320             value[i][k] /= scale;\r
 321             h += value[i][k]*value[i][k];\r
 322           }\r
 323           f = value[i][l];\r
 324           if (f>0) {\r
 325             g = -1.0*Math.sqrt(h);\r
 326           } else {\r
 327             g = Math.sqrt(h);\r
 328           }\r
 329           e[i] = scale*g;\r
 330           h -= f*g;\r
 331           value[i][l] = f-g;\r
 332           f=0.0;\r
 333           for (j=0; j < l; j++) {\r
 334             value[j][i] = value[i][j]/h;\r
 335             g=0.0;\r
 336             for (k= 0; k < j; k++) {\r
 337               g += value[j][k]*value[i][k];\r
 338             }\r
 339             for (k=j; k<l;k++) {\r
 340               g+= value[k][j]*value[i][k];\r
 341             }\r
 342             e[j] = g/h;\r
 343             f+=e[j]*value[i][j];\r
 344           }\r
 345           hh=f/(h+h);\r
 346           for (j=0;j<l;j++) {\r
 347             f=value[i][j];\r
 348             g = e[j] - hh*f;\r
 349             e[j] = g;\r
 350             for (k=0;k<j;k++) {\r
 351               value[j][k] -= (f*e[k]+g*value[i][k]);\r
 352             }\r
 353           }\r
 354         }\r
 355       } else {\r
 356         e[i] = value[i][l];\r
 357       }\r
 358       d[i] = h;\r
 359     }\r
 360     d[0] = 0.0;\r
 361     e[0] = 0.0;\r
 362     for (i=0;i<n;i++) {\r
 363       l=i-1;\r
 364       if (d[i] != 0.0) {\r
 365         for (j=0;j<l;j++) {\r
 366           g=0.0;\r
 367           for (k=0;k<l;k++) {\r
 368             g+= value[i][k]*value[k][j];\r
 369           }\r
 370           for (k=0;k<l;k++) {\r
 371             value[k][j] -= g*value[k][i];\r
 372           }\r
 373         }\r
 374       }\r
 375       d[i] = value[i][i];\r
 376       value[i][i] = 1.0;\r
 377       for (j=0;j<l;j++) {\r
 378         value[j][i] = 0.0;\r
 379         value[i][j] = 0.0;\r
 380       }\r
 381     }\r
 382   }\r
 383 \r
 384   public void tqli2() {\r
 385     int n = rows;\r
 386 \r
 387     int m;\r
 388     int l;\r
 389     int iter;\r
 390     int i;\r
 391     int k;\r
 392     double s;\r
 393     double r;\r
 394     double p;\r
 395     ;\r
 396     double g;\r
 397     double f;\r
 398     double dd;\r
 399     double c;\r
 400     double b;\r
 401 \r
 402     for (i=2;i<=n;i++) {\r
 403       e[i-2] = e[i-1];\r
 404     }\r
 405     e[n-1] = 0.0;\r
 406     for (l=1;l<=n;l++) {\r
 407       iter=0;\r
 408       do {\r
 409         for (m=l;m<=(n-1);m++) {\r
 410           dd=Math.abs(d[m-1]) + Math.abs(d[m]);\r
 411           if (Math.abs(e[m-1]) + dd == dd)\r
 412             break;\r
 413         }\r
 414         if (m != l) {\r
 415           iter++;\r
 416           if (iter == 30) {\r
 417             System.err.print("Too many iterations in tqli");\r
 418             System.exit(0); // JBPNote - same as above - not a graceful exit!\r
 419           } else {\r
 420             //      System.out.println("Iteration " + iter);\r
 421           }\r
 422           g=(d[l]-d[l-1])/(2.0*e[l-1]);\r
 423           r = Math.sqrt((g*g) + 1.0);\r
 424           g=d[m-1]-d[l-1]+e[l-1]/(g + sign(r,g));\r
 425           c = 1.0;\r
 426           s = c;\r
 427           p=0.0;\r
 428           for (i=m-1;i>=l;i--) {\r
 429             f = s*e[i-1];\r
 430             b = c*e[i-1];\r
 431             if (Math.abs(f) >= Math.abs(g)) {\r
 432               c=g/f;\r
 433               r = Math.sqrt((c*c)+1.0);\r
 434               e[i] = f*r;\r
 435               s = 1.0/r;\r
 436               c *= s;\r
 437             } else {\r
 438               s=f/g;\r
 439               r = Math.sqrt((s*s)+1.0);\r
 440               e[i] = g*r;\r
 441               c = 1.0/r;\r
 442               s *= c;\r
 443             }\r
 444             g=d[i] -p;\r
 445             r=(d[i-1]-g)*s + 2.0*c*b;\r
 446             p=s*r;\r
 447             d[i] = g + p;\r
 448             g = c * r - b;\r
 449             for (k=1; k <= n; k++) {\r
 450               f=value[k-1][i];\r
 451               value[k-1][i] = s*value[k-1][i-1] + c*f;\r
 452               value[k-1][i-1] = c*value[k-1][i-1] - s*f;\r
 453             }\r
 454           }\r
 455           d[l-1] = d[l-1] - p;\r
 456           e[l-1] = g;\r
 457           e[m-1] = 0.0;\r
 458         }\r
 459       } while ( m != l);\r
 460     }\r
 461   }\r
 462 \r
 463   public double sign(double a, double b) {\r
 464     if (b < 0) {\r
 465       return -Math.abs(a);\r
 466     } else {\r
 467       return Math.abs(a);\r
 468     }\r
 469   }\r
 470 \r
 471   public double[] getColumn(int n) {\r
 472     double[] out  = new double[rows];\r
 473     for (int i=0;i<rows;i++) {\r
 474       out[i] = value[i][n];\r
 475     }\r
 476     return out;\r
 477   }\r
 478 \r
 479 \r
 480   public void printD(PrintStream ps) {\r
 481 \r
 482     for (int j = 0; j < rows;j++) {\r
 483       Format.print(ps,"%15.4e",d[j]);\r
 484     }\r
 485   }\r
 486   public void printE(PrintStream ps) {\r
 487 \r
 488     for (int j = 0; j < rows;j++) {\r
 489       Format.print(ps,"%15.4e",e[j]);\r
 490     }\r
 491   }\r
 492 \r
 493   public static void main(String[] args) {\r
 494     int n = Integer.parseInt(args[0]);\r
 495     double[][] in = new double[n][n];\r
 496 \r
 497     for (int i = 0;i < n;i++) {\r
 498       for (int j = 0; j < n; j++) {\r
 499         in[i][j] = (double)Math.random();\r
 500       }\r
 501     }\r
 502 \r
 503     Matrix origmat = new Matrix(in,n,n);\r
 504     //    System.out.println(" --- Original matrix ---- ");\r
 505     ///    origmat.print(System.out);\r
 506     //System.out.println();\r
 507 \r
 508     //System.out.println(" --- transpose matrix ---- ");\r
 509     Matrix trans = origmat.transpose();\r
 510     //trans.print(System.out);\r
 511     //System.out.println();\r
 512 \r
 513     //System.out.println(" --- OrigT * Orig ---- ");\r
 514 \r
 515     Matrix symm = trans.postMultiply(origmat);\r
 516     //symm.print(System.out);\r
 517     //System.out.println();\r
 518 \r
 519     // Copy the symmetric matrix for later\r
 520     Matrix origsymm = symm.copy();\r
 521 \r
 522 \r
 523     // This produces the tridiagonal transformation matrix\r
 524     long tstart = System.currentTimeMillis();\r
 525     symm.tred();\r
 526     long tend = System.currentTimeMillis();\r
 527     //System.out.println("Time take for tred = " + (tend-tstart) + "ms");\r
 528     //System.out.println(" ---Tridiag transform matrix ---");\r
 529     //symm.print(System.out);\r
 530     //System.out.println();\r
 531 \r
 532     //System.out.println(" --- D vector ---");\r
 533     //symm.printD(System.out);\r
 534     //System.out.println();\r
 535     //System.out.println(" --- E vector ---");\r
 536     //symm.printE(System.out);\r
 537     //System.out.println();\r
 538 \r
 539 \r
 540     // Now produce the diagonalization matrix\r
 541     tstart = System.currentTimeMillis();\r
 542     symm.tqli();\r
 543     tend = System.currentTimeMillis();\r
 544     //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");\r
 545 \r
 546     //System.out.println(" --- New diagonalization matrix ---");\r
 547     //symm.print(System.out);\r
 548     //System.out.println();\r
 549 \r
 550     //System.out.println(" --- D vector ---");\r
 551     //symm.printD(System.out);\r
 552     //System.out.println();\r
 553     //System.out.println(" --- E vector ---");\r
 554     //symm.printE(System.out);\r
 555     //System.out.println();\r
 556 \r
 557     //System.out.println(" --- First eigenvector --- ");\r
 558     //double[] eigenv = symm.getColumn(0);\r
 559     //for (int i=0; i < eigenv.length;i++) {\r
 560     //  Format.print(System.out,"%15.4f",eigenv[i]);\r
 561     // }\r
 562     //System.out.println();\r
 563 \r
 564     //double[] neigenv = origsymm.vectorPostMultiply(eigenv);\r
 565 \r
 566     //for (int i=0; i < neigenv.length;i++) {\r
 567     //  Format.print(System.out,"%15.4f",neigenv[i]/symm.d[0]);\r
 568     //}\r
 569 \r
 570     //System.out.println();\r
 571   }\r
 572 \r
 573 }\r
 574 \r
 575 \r
 576 \r
 577 \r
 578 \r
 579 \r
 580 \r
 581 \r
 582 \r
 583 \r
 584 \r
 585 \r
 586 \r
 587 \r
 588 \r
 589 \r
 590 \r
 591 \r