1 // PCA.java with an irrelevant comment
\r
3 package jalview.analysis;
\r
5 import jalview.math.*;
\r
6 import jalview.datamodel.*;
\r
7 import jalview.util.*;
\r
12 public class PCA implements Runnable {
\r
17 double[] eigenvalue;
\r
20 public PCA(Matrix m) {
\r
24 public PCA(SequenceI[] s) {
\r
25 Runtime rt = Runtime.getRuntime();
\r
27 BinarySequence[] bs = new BinarySequence[s.length];
\r
29 while (ii < s.length && s[ii] != null) {
\r
31 bs[ii] = new BinarySequence(s[ii]);
\r
36 BinarySequence[] bs2 = new BinarySequence[s.length];
\r
38 while (ii < s.length && s[ii] != null) {
\r
40 bs2[ii] = new BinarySequence(s[ii]);
\r
41 bs2[ii].blosumEncode();
\r
46 //System.out.println("Created binary encoding");
\r
50 while (count < bs.length && bs[count] != null) {
\r
53 double[][] seqmat = new double[count][bs[0].getDBinary().length];
\r
54 double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
\r
57 seqmat[i] = bs[i].getDBinary();
\r
58 seqmat2[i] = bs2[i].getDBinary();
\r
61 //System.out.println("Created array");
\r
63 // System.out.println(" --- Original matrix ---- ");
\r
64 m = new Matrix(seqmat,count,bs[0].getDBinary().length);
\r
65 m2 = new Matrix(seqmat2,count,bs2[0].getDBinary().length);
\r
67 //System.out.println("Created matrix");
\r
71 public static void printMemory(Runtime rt) {
\r
72 System.out.println("Free memory = " + rt.freeMemory());
\r
75 public Matrix getM() {
\r
79 public double[] getEigenvector(int i) {
\r
80 return eigenvector.getColumn(i);
\r
83 public double getEigenvalue(int i) {
\r
84 return eigenvector.d[i];
\r
86 public float[][] getComponents(int l, int n, int mm) {
\r
87 return getComponents(l,n,mm,1);
\r
89 public float[][] getComponents(int l, int n, int mm, float factor) {
\r
90 float[][] out = new float[m.rows][3];
\r
92 for (int i = 0; i < m.rows;i++) {
\r
93 out[i][0] = (float)component(i,l)*factor;
\r
94 out[i][1] = (float)component(i,n)*factor;
\r
95 out[i][2] = (float)component(i,mm)*factor;
\r
100 public double[] component(int n) {
\r
101 // n = index of eigenvector
\r
102 double[] out = new double[m.rows];
\r
104 for (int i=0; i < m.rows; i++) {
\r
105 out[i] = component(i,n);
\r
109 public double component(int row, int n) {
\r
112 for (int i = 0; i < symm.cols; i++) {
\r
113 out += symm.value[row][i] * eigenvector.value[i][n];
\r
115 return out/eigenvector.d[n];
\r
118 public void checkEigenvector(int n,PrintStream ps) {
\r
119 ps.println(" --- Eigenvector " + n + " --- ");
\r
121 double[] eigenv = eigenvector.getColumn(n);
\r
123 for (int i=0; i < eigenv.length;i++) {
\r
124 Format.print(ps,"%15.4f",eigenv[i]);
\r
127 System.out.println();
\r
129 double[] neigenv = symm.vectorPostMultiply(eigenv);
\r
130 System.out.println(" --- symmat * eigenv / lambda --- ");
\r
131 if (eigenvector.d[n] > 1e-4) {
\r
132 for (int i=0; i < neigenv.length;i++) {
\r
133 Format.print(System.out,"%15.4f",neigenv[i]/eigenvector.d[n]);
\r
136 System.out.println();
\r
139 public void run() {
\r
140 Matrix mt = m.transpose();
\r
141 // System.out.println(" --- OrigT * Orig ---- ");
\r
142 eigenvector = mt.preMultiply(m2);
\r
143 // eigenvector.print(System.out);
\r
144 symm = eigenvector.copy();
\r
146 //TextArea ta = new TextArea(25,72);
\r
147 //TextAreaPrintStream taps = new TextAreaPrintStream(System.out,ta);
\r
148 //Frame f = new Frame("PCA output");
\r
149 //f.resize(500,500);
\r
150 //f.setLayout(new BorderLayout());
\r
151 //f.add("Center",ta);
\r
153 //symm.print(taps);
\r
154 long tstart = System.currentTimeMillis();
\r
155 eigenvector.tred();
\r
156 long tend = System.currentTimeMillis();
\r
157 //taps.println("Time take for tred = " + (tend-tstart) + "ms");
\r
158 //taps.println(" ---Tridiag transform matrix ---");
\r
160 //taps.println(" --- D vector ---");
\r
161 //eigenvector.printD(taps);
\r
163 //taps.println(" --- E vector ---");
\r
164 // eigenvector.printE(taps);
\r
167 // Now produce the diagonalization matrix
\r
168 tstart = System.currentTimeMillis();
\r
169 eigenvector.tqli();
\r
170 tend = System.currentTimeMillis();
\r
171 //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");
\r
173 //System.out.println(" --- New diagonalization matrix ---");
\r
175 //System.out.println(" --- Eigenvalues ---");
\r
176 //eigenvector.printD(taps);
\r
178 //System.out.println();
\r
180 // for (int i=0; i < eigenvector.cols; i++) {
\r
181 // checkEigenvector(i,taps);
\r
186 // taps.println("Transformed sequences = ");
\r
187 // Matrix trans = m.preMultiply(eigenvector);
\r
188 // trans.print(System.out);
\r