various loadTree bits for loading and saving trees.
[jalview.git] / src / jalview / analysis / PCA.java
1 // PCA.java \r
2 package jalview.analysis;\r
3 \r
4 import jalview.math.*;\r
5 import jalview.datamodel.*;\r
6 import jalview.util.*;\r
7 \r
8 import java.awt.*;\r
9 import java.io.*;\r
10 \r
11 public class PCA implements Runnable {\r
12   Matrix m;\r
13   Matrix symm;\r
14   Matrix m2;\r
15 \r
16   double[] eigenvalue;\r
17   Matrix eigenvector;\r
18 \r
19   public PCA(Matrix m) {\r
20     this.m = m;\r
21   }\r
22 \r
23   public PCA(SequenceI[] s) {\r
24     Runtime rt = Runtime.getRuntime();\r
25 \r
26     BinarySequence[] bs = new BinarySequence[s.length];\r
27     int ii = 0;\r
28     while (ii < s.length && s[ii] != null) {\r
29 \r
30       bs[ii] = new BinarySequence(s[ii]);\r
31       bs[ii].encode();\r
32       ii++;\r
33     }\r
34 \r
35     BinarySequence[] bs2 = new BinarySequence[s.length];\r
36     ii = 0;\r
37     while (ii < s.length && s[ii] != null) {\r
38 \r
39       bs2[ii] = new BinarySequence(s[ii]);\r
40       bs2[ii].blosumEncode();\r
41       ii++;\r
42     }\r
43 \r
44 \r
45     //System.out.println("Created binary encoding");\r
46     //printMemory(rt);\r
47 \r
48     int count=0;\r
49     while (count < bs.length && bs[count] != null) {\r
50       count++;\r
51     }\r
52     double[][] seqmat = new double[count][bs[0].getDBinary().length];\r
53     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];\r
54     int i=0;\r
55     while (i < count) {\r
56       seqmat[i] = bs[i].getDBinary();\r
57       seqmat2[i] = bs2[i].getDBinary();\r
58       i++;\r
59     }\r
60     //System.out.println("Created array");\r
61     //printMemory(rt);\r
62     //    System.out.println(" --- Original matrix ---- ");\r
63     m = new Matrix(seqmat,count,bs[0].getDBinary().length);\r
64     m2 = new Matrix(seqmat2,count,bs2[0].getDBinary().length);\r
65 \r
66     //System.out.println("Created matrix");\r
67     printMemory(rt);\r
68   }\r
69 \r
70   public static void printMemory(Runtime rt) {\r
71     System.out.println("Free memory = " + rt.freeMemory());\r
72   }\r
73 \r
74   public Matrix getM() {\r
75     return m;\r
76   }\r
77 \r
78   public double[] getEigenvector(int i) {\r
79     return eigenvector.getColumn(i);\r
80   }\r
81 \r
82   public double getEigenvalue(int i) {\r
83     return eigenvector.d[i];\r
84   }\r
85   public float[][] getComponents(int l, int n, int mm) {\r
86     return getComponents(l,n,mm,1);\r
87   }\r
88   public float[][] getComponents(int l, int n, int mm, float factor) {\r
89     float[][] out = new float[m.rows][3];\r
90 \r
91     for (int i = 0; i < m.rows;i++) {\r
92       out[i][0] = (float)component(i,l)*factor;\r
93       out[i][1] = (float)component(i,n)*factor;\r
94       out[i][2] = (float)component(i,mm)*factor;\r
95     }\r
96     return out;\r
97   }\r
98 \r
99   public double[] component(int n) {\r
100     // n = index of eigenvector\r
101     double[] out = new double[m.rows];\r
102 \r
103     for (int i=0; i < m.rows; i++) {\r
104       out[i] = component(i,n);\r
105     }\r
106     return out;\r
107   }\r
108   public double component(int row, int n) {\r
109     double out = 0.0;\r
110 \r
111     for (int i = 0; i < symm.cols; i++) {\r
112       out += symm.value[row][i] * eigenvector.value[i][n];\r
113     }\r
114     return out/eigenvector.d[n];\r
115   }\r
116 \r
117   public void checkEigenvector(int n,PrintStream ps) {\r
118     ps.println(" --- Eigenvector " + n  + " --- ");\r
119 \r
120     double[] eigenv = eigenvector.getColumn(n);\r
121 \r
122     for (int i=0; i < eigenv.length;i++) {\r
123       Format.print(ps,"%15.4f",eigenv[i]);\r
124     }\r
125 \r
126     System.out.println();\r
127 \r
128     double[] neigenv = symm.vectorPostMultiply(eigenv);\r
129     System.out.println(" --- symmat * eigenv / lambda --- ");\r
130     if (eigenvector.d[n] > 1e-4) {\r
131       for (int i=0; i < neigenv.length;i++) {\r
132         Format.print(System.out,"%15.4f",neigenv[i]/eigenvector.d[n]);\r
133       }\r
134     }\r
135     System.out.println();\r
136   }\r
137 \r
138   public void run() {\r
139     Matrix mt = m.transpose();\r
140     //    System.out.println(" --- OrigT * Orig ---- ");\r
141     eigenvector = mt.preMultiply(m2);\r
142     //  eigenvector.print(System.out);\r
143     symm = eigenvector.copy();\r
144 \r
145     //TextArea ta = new TextArea(25,72);\r
146     //TextAreaPrintStream taps = new TextAreaPrintStream(System.out,ta);\r
147     //Frame f = new Frame("PCA output");\r
148     //f.resize(500,500);\r
149     //f.setLayout(new BorderLayout());\r
150     //f.add("Center",ta);\r
151     //f.show();\r
152     //symm.print(taps);\r
153     long tstart = System.currentTimeMillis();\r
154     eigenvector.tred();\r
155     long tend = System.currentTimeMillis();\r
156     //taps.println("Time take for tred = " + (tend-tstart) + "ms");\r
157     //taps.println(" ---Tridiag transform matrix ---");\r
158 \r
159     //taps.println(" --- D vector ---");\r
160     //eigenvector.printD(taps);\r
161     //taps.println();\r
162     //taps.println(" --- E vector ---");\r
163     //    eigenvector.printE(taps);\r
164     //taps.println();\r
165 \r
166     // Now produce the diagonalization matrix\r
167     tstart = System.currentTimeMillis();\r
168     eigenvector.tqli();\r
169     tend = System.currentTimeMillis();\r
170     //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");\r
171 \r
172     //System.out.println(" --- New diagonalization matrix ---");\r
173 \r
174     //System.out.println(" --- Eigenvalues ---");\r
175     //eigenvector.printD(taps);\r
176 \r
177     //System.out.println();\r
178 \r
179     // for (int i=0; i < eigenvector.cols; i++) {\r
180     // checkEigenvector(i,taps);\r
181     // taps.println();\r
182     // }\r
183 \r
184     //  taps.println();\r
185     //  taps.println("Transformed sequences = ");\r
186     // Matrix trans =  m.preMultiply(eigenvector);\r
187     //  trans.print(System.out);\r
188   }\r
189 \r
190 }\r