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