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