JAL-1159 Merging 2.8.1 patches and Anne Menards work from 2012 into development branch
[jalview.git] / src / jalview / analysis / PCA.java
1 /*
2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.8)
3  * Copyright (C) 2012 J Procter, AM Waterhouse, LM Lui, J Engelhardt, 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.datamodel.BinarySequence.InvalidSequenceTypeException;
24 import jalview.math.*;
25 import jalview.schemes.ResidueProperties;
26 import jalview.schemes.ScoreMatrix;
27
28 /**
29  * Performs Principal Component Analysis on given sequences
30  * 
31  * @author $author$
32  * @version $Revision$
33  */
34 public class PCA implements Runnable
35 {
36   Matrix m;
37
38   Matrix symm;
39
40   Matrix m2;
41
42   double[] eigenvalue;
43
44   Matrix eigenvector;
45
46   StringBuffer details = new StringBuffer();
47
48   /**
49    * Creates a new PCA object. By default, uses blosum62 matrix to generate
50    * sequence similarity matrices
51    * 
52    * @param s
53    *          Set of amino acid sequences to perform PCA on
54    */
55   public PCA(String[] s)
56   {
57     this(s, false);
58   }
59
60   /**
61    * Creates a new PCA object. By default, uses blosum62 matrix to generate
62    * sequence similarity matrices
63    * 
64    * @param s
65    *          Set of sequences to perform PCA on
66    * @param nucleotides
67    *          if true, uses standard DNA/RNA matrix for sequence similarity
68    *          calculation.
69    */
70   public PCA(String[] s, boolean nucleotides)
71   {
72
73     BinarySequence[] bs = new BinarySequence[s.length];
74     int ii = 0;
75
76     while ((ii < s.length) && (s[ii] != null))
77     {
78       bs[ii] = new BinarySequence(s[ii], nucleotides);
79       bs[ii].encode();
80       ii++;
81     }
82
83     BinarySequence[] bs2 = new BinarySequence[s.length];
84     ii = 0;
85
86     String sm = nucleotides ? "DNA" : "BLOSUM62";
87     ScoreMatrix smtrx = ResidueProperties.getScoreMatrix(sm);
88     details.append("PCA calculation using " + sm
89             + " sequence similarity matrix\n========\n\n");
90     while ((ii < s.length) && (s[ii] != null))
91     {
92       bs2[ii] = new BinarySequence(s[ii], nucleotides);
93       if (smtrx != null)
94       {
95         try
96         {
97           bs2[ii].matrixEncode(smtrx);
98         } catch (InvalidSequenceTypeException x)
99         {
100           details.append("Unexpected mismatch of sequence type and score matrix. Calculation will not be valid!\n\n");
101         }
102       }
103       ii++;
104     }
105
106     // System.out.println("Created binary encoding");
107     // printMemory(rt);
108     int count = 0;
109
110     while ((count < bs.length) && (bs[count] != null))
111     {
112       count++;
113     }
114
115     double[][] seqmat = new double[count][bs[0].getDBinary().length];
116     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
117     int i = 0;
118
119     while (i < count)
120     {
121       seqmat[i] = bs[i].getDBinary();
122       seqmat2[i] = bs2[i].getDBinary();
123       i++;
124     }
125
126     // System.out.println("Created array");
127     // printMemory(rt);
128     // System.out.println(" --- Original matrix ---- ");
129     m = new Matrix(seqmat, count, bs[0].getDBinary().length);
130     m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length);
131
132   }
133
134   /**
135    * Returns the matrix used in PCA calculation
136    * 
137    * @return java.math.Matrix object
138    */
139
140   public Matrix getM()
141   {
142     return m;
143   }
144
145   /**
146    * Returns Eigenvalue
147    * 
148    * @param i
149    *          Index of diagonal within matrix
150    * 
151    * @return Returns value of diagonal from matrix
152    */
153   public double getEigenvalue(int i)
154   {
155     return eigenvector.d[i];
156   }
157
158   /**
159    * DOCUMENT ME!
160    * 
161    * @param l
162    *          DOCUMENT ME!
163    * @param n
164    *          DOCUMENT ME!
165    * @param mm
166    *          DOCUMENT ME!
167    * @param factor
168    *          DOCUMENT ME!
169    * 
170    * @return DOCUMENT ME!
171    */
172   public float[][] getComponents(int l, int n, int mm, float factor)
173   {
174     float[][] out = new float[m.rows][3];
175
176     for (int i = 0; i < m.rows; i++)
177     {
178       out[i][0] = (float) component(i, l) * factor;
179       out[i][1] = (float) component(i, n) * factor;
180       out[i][2] = (float) component(i, mm) * factor;
181     }
182
183     return out;
184   }
185
186   /**
187    * DOCUMENT ME!
188    * 
189    * @param n
190    *          DOCUMENT ME!
191    * 
192    * @return DOCUMENT ME!
193    */
194   public double[] component(int n)
195   {
196     // n = index of eigenvector
197     double[] out = new double[m.rows];
198
199     for (int i = 0; i < m.rows; i++)
200     {
201       out[i] = component(i, n);
202     }
203
204     return out;
205   }
206
207   /**
208    * DOCUMENT ME!
209    * 
210    * @param row
211    *          DOCUMENT ME!
212    * @param n
213    *          DOCUMENT ME!
214    * 
215    * @return DOCUMENT ME!
216    */
217   double component(int row, int n)
218   {
219     double out = 0.0;
220
221     for (int i = 0; i < symm.cols; i++)
222     {
223       out += (symm.value[row][i] * eigenvector.value[i][n]);
224     }
225
226     return out / eigenvector.d[n];
227   }
228
229   public String getDetails()
230   {
231     return details.toString();
232   }
233
234   /**
235    * DOCUMENT ME!
236    */
237   public void run()
238   {
239     details.append("PCA Calculation Mode is "
240             + (jvCalcMode ? "Jalview variant" : "Original SeqSpace") + "\n");
241     Matrix mt = m.transpose();
242
243     details.append(" --- OrigT * Orig ---- \n");
244     if (!jvCalcMode)
245     {
246       eigenvector = mt.preMultiply(m); // standard seqspace comparison matrix
247     }
248     else
249     {
250       eigenvector = mt.preMultiply(m2); // jalview variation on seqsmace method
251     }
252
253     PrintStream ps = new PrintStream(System.out)
254     {
255       public void print(String x)
256       {
257         details.append(x);
258       }
259
260       public void println()
261       {
262         details.append("\n");
263       }
264     };
265
266     eigenvector.print(ps);
267
268     symm = eigenvector.copy();
269
270     eigenvector.tred();
271
272     details.append(" ---Tridiag transform matrix ---\n");
273     details.append(" --- D vector ---\n");
274     eigenvector.printD(ps);
275     ps.println();
276     details.append("--- E vector ---\n");
277     eigenvector.printE(ps);
278     ps.println();
279
280     // Now produce the diagonalization matrix
281     eigenvector.tqli();
282
283     details.append(" --- New diagonalization matrix ---\n");
284     eigenvector.print(ps);
285     details.append(" --- Eigenvalues ---\n");
286     eigenvector.printD(ps);
287     ps.println();
288     /*
289      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
290      * (int ev=0;ev<symm.rows; ev++) {
291      * 
292      * ps.print(","+component(seq, ev)); } ps.println(); }
293      */
294   }
295
296   boolean jvCalcMode = true;
297
298   public void setJvCalcMode(boolean calcMode)
299   {
300     this.jvCalcMode = calcMode;
301   }
302 }