added help and documentation for annotation editing.
[jalview.git] / src / jalview / analysis / PCA.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 package jalview.analysis;\r
20 \r
21 import jalview.datamodel.*;\r
22 \r
23 import jalview.math.*;\r
24 \r
25 import jalview.util.*;\r
26 \r
27 import java.io.*;\r
28 \r
29 \r
30 /**\r
31  * Performs Principal Component Analysis on given sequences\r
32  *\r
33  * @author $author$\r
34  * @version $Revision$\r
35  */\r
36 public class PCA implements Runnable\r
37 {\r
38     Matrix m;\r
39     Matrix symm;\r
40     Matrix m2;\r
41     double[] eigenvalue;\r
42     Matrix eigenvector;\r
43 \r
44 \r
45     /**\r
46      * Creates a new PCA object.\r
47      *\r
48      * @param s Set of sequences to perform PCA on\r
49      */\r
50     public PCA(SequenceI[] s)\r
51     {\r
52 \r
53         BinarySequence[] bs = new BinarySequence[s.length];\r
54         int ii = 0;\r
55 \r
56         while ((ii < s.length) && (s[ii] != null))\r
57         {\r
58             bs[ii] = new BinarySequence(s[ii]);\r
59             bs[ii].encode();\r
60             ii++;\r
61         }\r
62 \r
63         BinarySequence[] bs2 = new BinarySequence[s.length];\r
64         ii = 0;\r
65 \r
66         while ((ii < s.length) && (s[ii] != null))\r
67         {\r
68             bs2[ii] = new BinarySequence(s[ii]);\r
69             bs2[ii].blosumEncode();\r
70             ii++;\r
71         }\r
72 \r
73         //System.out.println("Created binary encoding");\r
74         //printMemory(rt);\r
75         int count = 0;\r
76 \r
77         while ((count < bs.length) && (bs[count] != null))\r
78         {\r
79             count++;\r
80         }\r
81 \r
82         double[][] seqmat = new double[count][bs[0].getDBinary().length];\r
83         double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];\r
84         int i = 0;\r
85 \r
86         while (i < count)\r
87         {\r
88             seqmat[i] = bs[i].getDBinary();\r
89             seqmat2[i] = bs2[i].getDBinary();\r
90             i++;\r
91         }\r
92 \r
93         //System.out.println("Created array");\r
94         //printMemory(rt);\r
95         //    System.out.println(" --- Original matrix ---- ");\r
96         m = new Matrix(seqmat, count, bs[0].getDBinary().length);\r
97         m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length);\r
98 \r
99       }\r
100 \r
101       /**\r
102        * Returns the matrix used in PCA calculation\r
103        *\r
104        * @return java.math.Matrix object\r
105        */\r
106 \r
107       public Matrix getM()\r
108       {\r
109         return m;\r
110       }\r
111 \r
112     /**\r
113      * Returns Eigenvalue\r
114      *\r
115      * @param i Index of diagonal within matrix\r
116      *\r
117      * @return Returns value of diagonal from matrix\r
118      */\r
119     public double getEigenvalue(int i)\r
120     {\r
121         return eigenvector.d[i];\r
122     }\r
123 \r
124      /**\r
125      * DOCUMENT ME!\r
126      *\r
127      * @param l DOCUMENT ME!\r
128      * @param n DOCUMENT ME!\r
129      * @param mm DOCUMENT ME!\r
130      * @param factor DOCUMENT ME!\r
131      *\r
132      * @return DOCUMENT ME!\r
133      */\r
134     public float[][] getComponents(int l, int n, int mm, float factor)\r
135     {\r
136         float[][] out = new float[m.rows][3];\r
137 \r
138         for (int i = 0; i < m.rows; i++)\r
139         {\r
140             out[i][0] = (float) component(i, l) * factor;\r
141             out[i][1] = (float) component(i, n) * factor;\r
142             out[i][2] = (float) component(i, mm) * factor;\r
143         }\r
144 \r
145         return out;\r
146     }\r
147 \r
148     /**\r
149      * DOCUMENT ME!\r
150      *\r
151      * @param n DOCUMENT ME!\r
152      *\r
153      * @return DOCUMENT ME!\r
154      */\r
155     public double[] component(int n)\r
156     {\r
157         // n = index of eigenvector\r
158         double[] out = new double[m.rows];\r
159 \r
160         for (int i = 0; i < m.rows; i++)\r
161         {\r
162             out[i] = component(i, n);\r
163         }\r
164 \r
165         return out;\r
166     }\r
167 \r
168     /**\r
169      * DOCUMENT ME!\r
170      *\r
171      * @param row DOCUMENT ME!\r
172      * @param n DOCUMENT ME!\r
173      *\r
174      * @return DOCUMENT ME!\r
175      */\r
176     double component(int row, int n)\r
177     {\r
178         double out = 0.0;\r
179 \r
180         for (int i = 0; i < symm.cols; i++)\r
181         {\r
182             out += (symm.value[row][i] * eigenvector.value[i][n]);\r
183         }\r
184 \r
185         return out / eigenvector.d[n];\r
186     }\r
187 \r
188 \r
189     /**\r
190      * DOCUMENT ME!\r
191      */\r
192     public void run()\r
193     {\r
194         Matrix mt = m.transpose();\r
195 \r
196         //    System.out.println(" --- OrigT * Orig ---- ");\r
197         eigenvector = mt.preMultiply(m2);\r
198 \r
199         //  eigenvector.print(System.out);\r
200         symm = eigenvector.copy();\r
201 \r
202         //TextArea ta = new TextArea(25,72);\r
203         //TextAreaPrintStream taps = new TextAreaPrintStream(System.out,ta);\r
204         //Frame f = new Frame("PCA output");\r
205         //f.resize(500,500);\r
206         //f.setLayout(new BorderLayout());\r
207         //f.add("Center",ta);\r
208         //f.show();\r
209         //symm.print(taps);\r
210         long tstart = System.currentTimeMillis();\r
211         eigenvector.tred();\r
212 \r
213         long tend = System.currentTimeMillis();\r
214 \r
215         //taps.println("Time take for tred = " + (tend-tstart) + "ms");\r
216         //taps.println(" ---Tridiag transform matrix ---");\r
217         //taps.println(" --- D vector ---");\r
218         //eigenvector.printD(taps);\r
219         //taps.println();\r
220         //taps.println(" --- E vector ---");\r
221         //    eigenvector.printE(taps);\r
222         //taps.println();\r
223         // Now produce the diagonalization matrix\r
224         tstart = System.currentTimeMillis();\r
225         eigenvector.tqli();\r
226         tend = System.currentTimeMillis();\r
227 \r
228         //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");\r
229         //System.out.println(" --- New diagonalization matrix ---");\r
230         //System.out.println(" --- Eigenvalues ---");\r
231         //eigenvector.printD(taps);\r
232         //System.out.println();\r
233         // for (int i=0; i < eigenvector.cols; i++) {\r
234         // checkEigenvector(i,taps);\r
235         // taps.println();\r
236         // }\r
237         //  taps.println();\r
238         //  taps.println("Transformed sequences = ");\r
239         // Matrix trans =  m.preMultiply(eigenvector);\r
240         //  trans.print(System.out);\r
241     }\r
242 }\r