src/jalview/analysis/PCA.java

   1 /*\r
   2 * Jalview - A Sequence Alignment Editor and Viewer\r
   3 * Copyright (C) 2005 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle\r
   4 *\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
   9 *\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
  14 *\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
  18 */\r
  19 \r
  20 package jalview.analysis;\r
  21 \r
  22 import jalview.math.*;\r
  23 import jalview.datamodel.*;\r
  24 import jalview.util.*;\r
  25 \r
  26 import java.awt.*;\r
  27 import java.io.*;\r
  28 \r
  29 public class PCA implements Runnable {\r
  30   Matrix m;\r
  31   Matrix symm;\r
  32   Matrix m2;\r
  33 \r
  34   double[] eigenvalue;\r
  35   Matrix eigenvector;\r
  36 \r
  37   public PCA(Matrix m) {\r
  38     this.m = m;\r
  39   }\r
  40 \r
  41   public PCA(SequenceI[] s) {\r
  42     Runtime rt = Runtime.getRuntime();\r
  43 \r
  44     BinarySequence[] bs = new BinarySequence[s.length];\r
  45     int ii = 0;\r
  46     while (ii < s.length && s[ii] != null) {\r
  47 \r
  48       bs[ii] = new BinarySequence(s[ii]);\r
  49       bs[ii].encode();\r
  50       ii++;\r
  51     }\r
  52 \r
  53     BinarySequence[] bs2 = new BinarySequence[s.length];\r
  54     ii = 0;\r
  55     while (ii < s.length && s[ii] != null) {\r
  56 \r
  57       bs2[ii] = new BinarySequence(s[ii]);\r
  58       bs2[ii].blosumEncode();\r
  59       ii++;\r
  60     }\r
  61 \r
  62 \r
  63     //System.out.println("Created binary encoding");\r
  64     //printMemory(rt);\r
  65 \r
  66     int count=0;\r
  67     while (count < bs.length && bs[count] != null) {\r
  68       count++;\r
  69     }\r
  70     double[][] seqmat = new double[count][bs[0].getDBinary().length];\r
  71     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];\r
  72     int i=0;\r
  73     while (i < count) {\r
  74       seqmat[i] = bs[i].getDBinary();\r
  75       seqmat2[i] = bs2[i].getDBinary();\r
  76       i++;\r
  77     }\r
  78     //System.out.println("Created array");\r
  79     //printMemory(rt);\r
  80     //    System.out.println(" --- Original matrix ---- ");\r
  81     m = new Matrix(seqmat,count,bs[0].getDBinary().length);\r
  82     m2 = new Matrix(seqmat2,count,bs2[0].getDBinary().length);\r
  83 \r
  84     //System.out.println("Created matrix");\r
  85     printMemory(rt);\r
  86   }\r
  87 \r
  88   public static void printMemory(Runtime rt) {\r
  89     System.out.println("PCA:Free memory = " + rt.freeMemory());\r
  90   }\r
  91 \r
  92   public Matrix getM() {\r
  93     return m;\r
  94   }\r
  95 \r
  96   public double[] getEigenvector(int i) {\r
  97     return eigenvector.getColumn(i);\r
  98   }\r
  99 \r
 100   public double getEigenvalue(int i) {\r
 101     return eigenvector.d[i];\r
 102   }\r
 103   public float[][] getComponents(int l, int n, int mm) {\r
 104     return getComponents(l,n,mm,1);\r
 105   }\r
 106   public float[][] getComponents(int l, int n, int mm, float factor) {\r
 107     float[][] out = new float[m.rows][3];\r
 108 \r
 109     for (int i = 0; i < m.rows;i++) {\r
 110       out[i][0] = (float)component(i,l)*factor;\r
 111       out[i][1] = (float)component(i,n)*factor;\r
 112       out[i][2] = (float)component(i,mm)*factor;\r
 113     }\r
 114     return out;\r
 115   }\r
 116 \r
 117   public double[] component(int n) {\r
 118     // n = index of eigenvector\r
 119     double[] out = new double[m.rows];\r
 120 \r
 121     for (int i=0; i < m.rows; i++) {\r
 122       out[i] = component(i,n);\r
 123     }\r
 124     return out;\r
 125   }\r
 126   public double component(int row, int n) {\r
 127     double out = 0.0;\r
 128 \r
 129     for (int i = 0; i < symm.cols; i++) {\r
 130       out += symm.value[row][i] * eigenvector.value[i][n];\r
 131     }\r
 132     return out/eigenvector.d[n];\r
 133   }\r
 134 \r
 135   public void checkEigenvector(int n,PrintStream ps) {\r
 136     ps.println(" --- Eigenvector " + n  + " --- ");\r
 137 \r
 138     double[] eigenv = eigenvector.getColumn(n);\r
 139 \r
 140     for (int i=0; i < eigenv.length;i++) {\r
 141       Format.print(ps,"%15.4f",eigenv[i]);\r
 142     }\r
 143 \r
 144     System.out.println();\r
 145 \r
 146     double[] neigenv = symm.vectorPostMultiply(eigenv);\r
 147     System.out.println(" --- symmat * eigenv / lambda --- ");\r
 148     if (eigenvector.d[n] > 1e-4) {\r
 149       for (int i=0; i < neigenv.length;i++) {\r
 150         Format.print(System.out,"%15.4f",neigenv[i]/eigenvector.d[n]);\r
 151       }\r
 152     }\r
 153     System.out.println();\r
 154   }\r
 155 \r
 156   public void run() {\r
 157     Matrix mt = m.transpose();\r
 158     //    System.out.println(" --- OrigT * Orig ---- ");\r
 159     eigenvector = mt.preMultiply(m2);\r
 160     //  eigenvector.print(System.out);\r
 161     symm = eigenvector.copy();\r
 162 \r
 163     //TextArea ta = new TextArea(25,72);\r
 164     //TextAreaPrintStream taps = new TextAreaPrintStream(System.out,ta);\r
 165     //Frame f = new Frame("PCA output");\r
 166     //f.resize(500,500);\r
 167     //f.setLayout(new BorderLayout());\r
 168     //f.add("Center",ta);\r
 169     //f.show();\r
 170     //symm.print(taps);\r
 171     long tstart = System.currentTimeMillis();\r
 172     eigenvector.tred();\r
 173     long tend = System.currentTimeMillis();\r
 174     //taps.println("Time take for tred = " + (tend-tstart) + "ms");\r
 175     //taps.println(" ---Tridiag transform matrix ---");\r
 176 \r
 177     //taps.println(" --- D vector ---");\r
 178     //eigenvector.printD(taps);\r
 179     //taps.println();\r
 180     //taps.println(" --- E vector ---");\r
 181     //    eigenvector.printE(taps);\r
 182     //taps.println();\r
 183 \r
 184     // Now produce the diagonalization matrix\r
 185     tstart = System.currentTimeMillis();\r
 186     eigenvector.tqli();\r
 187     tend = System.currentTimeMillis();\r
 188     //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");\r
 189 \r
 190     //System.out.println(" --- New diagonalization matrix ---");\r
 191 \r
 192     //System.out.println(" --- Eigenvalues ---");\r
 193     //eigenvector.printD(taps);\r
 194 \r
 195     //System.out.println();\r
 196 \r
 197     // for (int i=0; i < eigenvector.cols; i++) {\r
 198     // checkEigenvector(i,taps);\r
 199     // taps.println();\r
 200     // }\r
 201 \r
 202     //  taps.println();\r
 203     //  taps.println("Transformed sequences = ");\r
 204     // Matrix trans =  m.preMultiply(eigenvector);\r
 205     //  trans.print(System.out);\r
 206   }\r
 207 \r
 208 }\r