2 * Jalview - A Sequence Alignment Editor and Viewer
\r
3 * Copyright (C) 2005 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle
\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
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
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
19 package jalview.analysis;
\r
21 import jalview.datamodel.*;
\r
23 import jalview.math.*;
\r
27 * Performs Principal Component Analysis on given sequences
\r
30 * @version $Revision$
\r
32 public class PCA implements Runnable
\r
37 double[] eigenvalue;
\r
42 * Creates a new PCA object.
\r
44 * @param s Set of sequences to perform PCA on
\r
46 public PCA(SequenceI[] s)
\r
49 BinarySequence[] bs = new BinarySequence[s.length];
\r
52 while ((ii < s.length) && (s[ii] != null))
\r
54 bs[ii] = new BinarySequence(s[ii]);
\r
59 BinarySequence[] bs2 = new BinarySequence[s.length];
\r
62 while ((ii < s.length) && (s[ii] != null))
\r
64 bs2[ii] = new BinarySequence(s[ii]);
\r
65 bs2[ii].blosumEncode();
\r
69 //System.out.println("Created binary encoding");
\r
73 while ((count < bs.length) && (bs[count] != null))
\r
78 double[][] seqmat = new double[count][bs[0].getDBinary().length];
\r
79 double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
\r
84 seqmat[i] = bs[i].getDBinary();
\r
85 seqmat2[i] = bs2[i].getDBinary();
\r
89 //System.out.println("Created array");
\r
91 // System.out.println(" --- Original matrix ---- ");
\r
92 m = new Matrix(seqmat, count, bs[0].getDBinary().length);
\r
93 m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length);
\r
98 * Returns the matrix used in PCA calculation
\r
100 * @return java.math.Matrix object
\r
103 public Matrix getM()
\r
109 * Returns Eigenvalue
\r
111 * @param i Index of diagonal within matrix
\r
113 * @return Returns value of diagonal from matrix
\r
115 public double getEigenvalue(int i)
\r
117 return eigenvector.d[i];
\r
123 * @param l DOCUMENT ME!
\r
124 * @param n DOCUMENT ME!
\r
125 * @param mm DOCUMENT ME!
\r
126 * @param factor DOCUMENT ME!
\r
128 * @return DOCUMENT ME!
\r
130 public float[][] getComponents(int l, int n, int mm, float factor)
\r
132 float[][] out = new float[m.rows][3];
\r
134 for (int i = 0; i < m.rows; i++)
\r
136 out[i][0] = (float) component(i, l) * factor;
\r
137 out[i][1] = (float) component(i, n) * factor;
\r
138 out[i][2] = (float) component(i, mm) * factor;
\r
147 * @param n DOCUMENT ME!
\r
149 * @return DOCUMENT ME!
\r
151 public double[] component(int n)
\r
153 // n = index of eigenvector
\r
154 double[] out = new double[m.rows];
\r
156 for (int i = 0; i < m.rows; i++)
\r
158 out[i] = component(i, n);
\r
167 * @param row DOCUMENT ME!
\r
168 * @param n DOCUMENT ME!
\r
170 * @return DOCUMENT ME!
\r
172 double component(int row, int n)
\r
176 for (int i = 0; i < symm.cols; i++)
\r
178 out += (symm.value[row][i] * eigenvector.value[i][n]);
\r
181 return out / eigenvector.d[n];
\r
190 Matrix mt = m.transpose();
\r
192 // System.out.println(" --- OrigT * Orig ---- ");
\r
193 eigenvector = mt.preMultiply(m2);
\r
195 // eigenvector.print(System.out);
\r
196 symm = eigenvector.copy();
\r
198 //TextArea ta = new TextArea(25,72);
\r
199 //TextAreaPrintStream taps = new TextAreaPrintStream(System.out,ta);
\r
200 //Frame f = new Frame("PCA output");
\r
201 //f.resize(500,500);
\r
202 //f.setLayout(new BorderLayout());
\r
203 //f.add("Center",ta);
\r
205 //symm.print(taps);
\r
206 long tstart = System.currentTimeMillis();
\r
207 eigenvector.tred();
\r
209 long tend = System.currentTimeMillis();
\r
211 //taps.println("Time take for tred = " + (tend-tstart) + "ms");
\r
212 //taps.println(" ---Tridiag transform matrix ---");
\r
213 //taps.println(" --- D vector ---");
\r
214 //eigenvector.printD(taps);
\r
216 //taps.println(" --- E vector ---");
\r
217 // eigenvector.printE(taps);
\r
219 // Now produce the diagonalization matrix
\r
220 tstart = System.currentTimeMillis();
\r
221 eigenvector.tqli();
\r
222 tend = System.currentTimeMillis();
\r
224 //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");
\r
225 //System.out.println(" --- New diagonalization matrix ---");
\r
226 //System.out.println(" --- Eigenvalues ---");
\r
227 //eigenvector.printD(taps);
\r
228 //System.out.println();
\r
229 // for (int i=0; i < eigenvector.cols; i++) {
\r
230 // checkEigenvector(i,taps);
\r
234 // taps.println("Transformed sequences = ");
\r
235 // Matrix trans = m.preMultiply(eigenvector);
\r
236 // trans.print(System.out);
\r