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
25 import jalview.util.*;
\r
32 public class PCA implements Runnable {
\r
36 double[] eigenvalue;
\r
39 public PCA(Matrix m) {
\r
43 public PCA(SequenceI[] s) {
\r
44 Runtime rt = Runtime.getRuntime();
\r
46 BinarySequence[] bs = new BinarySequence[s.length];
\r
49 while ((ii < s.length) && (s[ii] != null)) {
\r
50 bs[ii] = new BinarySequence(s[ii]);
\r
55 BinarySequence[] bs2 = new BinarySequence[s.length];
\r
58 while ((ii < s.length) && (s[ii] != null)) {
\r
59 bs2[ii] = new BinarySequence(s[ii]);
\r
60 bs2[ii].blosumEncode();
\r
64 //System.out.println("Created binary encoding");
\r
68 while ((count < bs.length) && (bs[count] != null)) {
\r
72 double[][] seqmat = new double[count][bs[0].getDBinary().length];
\r
73 double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
\r
77 seqmat[i] = bs[i].getDBinary();
\r
78 seqmat2[i] = bs2[i].getDBinary();
\r
82 //System.out.println("Created array");
\r
84 // System.out.println(" --- Original matrix ---- ");
\r
85 m = new Matrix(seqmat, count, bs[0].getDBinary().length);
\r
86 m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length);
\r
88 //System.out.println("Created matrix");
\r
92 public static void printMemory(Runtime rt) {
\r
93 System.out.println("PCA:Free memory = " + rt.freeMemory());
\r
96 public Matrix getM() {
\r
100 public double[] getEigenvector(int i) {
\r
101 return eigenvector.getColumn(i);
\r
104 public double getEigenvalue(int i) {
\r
105 return eigenvector.d[i];
\r
108 public float[][] getComponents(int l, int n, int mm) {
\r
109 return getComponents(l, n, mm, 1);
\r
112 public float[][] getComponents(int l, int n, int mm, float factor) {
\r
113 float[][] out = new float[m.rows][3];
\r
115 for (int i = 0; i < m.rows; i++) {
\r
116 out[i][0] = (float) component(i, l) * factor;
\r
117 out[i][1] = (float) component(i, n) * factor;
\r
118 out[i][2] = (float) component(i, mm) * factor;
\r
124 public double[] component(int n) {
\r
125 // n = index of eigenvector
\r
126 double[] out = new double[m.rows];
\r
128 for (int i = 0; i < m.rows; i++) {
\r
129 out[i] = component(i, n);
\r
135 public double component(int row, int n) {
\r
138 for (int i = 0; i < symm.cols; i++) {
\r
139 out += (symm.value[row][i] * eigenvector.value[i][n]);
\r
142 return out / eigenvector.d[n];
\r
145 public void checkEigenvector(int n, PrintStream ps) {
\r
146 ps.println(" --- Eigenvector " + n + " --- ");
\r
148 double[] eigenv = eigenvector.getColumn(n);
\r
150 for (int i = 0; i < eigenv.length; i++) {
\r
151 Format.print(ps, "%15.4f", eigenv[i]);
\r
154 System.out.println();
\r
156 double[] neigenv = symm.vectorPostMultiply(eigenv);
\r
157 System.out.println(" --- symmat * eigenv / lambda --- ");
\r
159 if (eigenvector.d[n] > 1e-4) {
\r
160 for (int i = 0; i < neigenv.length; i++) {
\r
161 Format.print(System.out, "%15.4f", neigenv[i] / eigenvector.d[n]);
\r
165 System.out.println();
\r
168 public void run() {
\r
169 Matrix mt = m.transpose();
\r
171 // System.out.println(" --- OrigT * Orig ---- ");
\r
172 eigenvector = mt.preMultiply(m2);
\r
174 // eigenvector.print(System.out);
\r
175 symm = eigenvector.copy();
\r
177 //TextArea ta = new TextArea(25,72);
\r
178 //TextAreaPrintStream taps = new TextAreaPrintStream(System.out,ta);
\r
179 //Frame f = new Frame("PCA output");
\r
180 //f.resize(500,500);
\r
181 //f.setLayout(new BorderLayout());
\r
182 //f.add("Center",ta);
\r
184 //symm.print(taps);
\r
185 long tstart = System.currentTimeMillis();
\r
186 eigenvector.tred();
\r
188 long tend = System.currentTimeMillis();
\r
190 //taps.println("Time take for tred = " + (tend-tstart) + "ms");
\r
191 //taps.println(" ---Tridiag transform matrix ---");
\r
192 //taps.println(" --- D vector ---");
\r
193 //eigenvector.printD(taps);
\r
195 //taps.println(" --- E vector ---");
\r
196 // eigenvector.printE(taps);
\r
198 // Now produce the diagonalization matrix
\r
199 tstart = System.currentTimeMillis();
\r
200 eigenvector.tqli();
\r
201 tend = System.currentTimeMillis();
\r
203 //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");
\r
204 //System.out.println(" --- New diagonalization matrix ---");
\r
205 //System.out.println(" --- Eigenvalues ---");
\r
206 //eigenvector.printD(taps);
\r
207 //System.out.println();
\r
208 // for (int i=0; i < eigenvector.cols; i++) {
\r
209 // checkEigenvector(i,taps);
\r
213 // taps.println("Transformed sequences = ");
\r
214 // Matrix trans = m.preMultiply(eigenvector);
\r
215 // trans.print(System.out);
\r