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