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 package jalview.analysis;\r
  20 \r
  21 import jalview.datamodel.*;\r
  22 \r
  23 import jalview.math.*;\r
  24 \r
  25 import jalview.util.*;\r
  26 \r
  27 import java.awt.*;\r
  28 \r
  29 import java.io.*;\r
  30 \r
  31 \r
  32 public class PCA implements Runnable {\r
  33     Matrix m;\r
  34     Matrix symm;\r
  35     Matrix m2;\r
  36     double[] eigenvalue;\r
  37     Matrix eigenvector;\r
  38 \r
  39     public PCA(Matrix m) {\r
  40         this.m = m;\r
  41     }\r
  42 \r
  43     public PCA(SequenceI[] s) {\r
  44         Runtime rt = Runtime.getRuntime();\r
  45 \r
  46         BinarySequence[] bs = new BinarySequence[s.length];\r
  47         int ii = 0;\r
  48 \r
  49         while ((ii < s.length) && (s[ii] != null)) {\r
  50             bs[ii] = new BinarySequence(s[ii]);\r
  51             bs[ii].encode();\r
  52             ii++;\r
  53         }\r
  54 \r
  55         BinarySequence[] bs2 = new BinarySequence[s.length];\r
  56         ii = 0;\r
  57 \r
  58         while ((ii < s.length) && (s[ii] != null)) {\r
  59             bs2[ii] = new BinarySequence(s[ii]);\r
  60             bs2[ii].blosumEncode();\r
  61             ii++;\r
  62         }\r
  63 \r
  64         //System.out.println("Created binary encoding");\r
  65         //printMemory(rt);\r
  66         int count = 0;\r
  67 \r
  68         while ((count < bs.length) && (bs[count] != null)) {\r
  69             count++;\r
  70         }\r
  71 \r
  72         double[][] seqmat = new double[count][bs[0].getDBinary().length];\r
  73         double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];\r
  74         int i = 0;\r
  75 \r
  76         while (i < count) {\r
  77             seqmat[i] = bs[i].getDBinary();\r
  78             seqmat2[i] = bs2[i].getDBinary();\r
  79             i++;\r
  80         }\r
  81 \r
  82         //System.out.println("Created array");\r
  83         //printMemory(rt);\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
  87 \r
  88         //System.out.println("Created matrix");\r
  89         printMemory(rt);\r
  90     }\r
  91 \r
  92     public static void printMemory(Runtime rt) {\r
  93         System.out.println("PCA:Free memory = " + rt.freeMemory());\r
  94     }\r
  95 \r
  96     public Matrix getM() {\r
  97         return m;\r
  98     }\r
  99 \r
 100     public double[] getEigenvector(int i) {\r
 101         return eigenvector.getColumn(i);\r
 102     }\r
 103 \r
 104     public double getEigenvalue(int i) {\r
 105         return eigenvector.d[i];\r
 106     }\r
 107 \r
 108     public float[][] getComponents(int l, int n, int mm) {\r
 109         return getComponents(l, n, mm, 1);\r
 110     }\r
 111 \r
 112     public float[][] getComponents(int l, int n, int mm, float factor) {\r
 113         float[][] out = new float[m.rows][3];\r
 114 \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
 119         }\r
 120 \r
 121         return out;\r
 122     }\r
 123 \r
 124     public double[] component(int n) {\r
 125         // n = index of eigenvector\r
 126         double[] out = new double[m.rows];\r
 127 \r
 128         for (int i = 0; i < m.rows; i++) {\r
 129             out[i] = component(i, n);\r
 130         }\r
 131 \r
 132         return out;\r
 133     }\r
 134 \r
 135     public double component(int row, int n) {\r
 136         double out = 0.0;\r
 137 \r
 138         for (int i = 0; i < symm.cols; i++) {\r
 139             out += (symm.value[row][i] * eigenvector.value[i][n]);\r
 140         }\r
 141 \r
 142         return out / eigenvector.d[n];\r
 143     }\r
 144 \r
 145     public void checkEigenvector(int n, PrintStream ps) {\r
 146         ps.println(" --- Eigenvector " + n + " --- ");\r
 147 \r
 148         double[] eigenv = eigenvector.getColumn(n);\r
 149 \r
 150         for (int i = 0; i < eigenv.length; i++) {\r
 151             Format.print(ps, "%15.4f", eigenv[i]);\r
 152         }\r
 153 \r
 154         System.out.println();\r
 155 \r
 156         double[] neigenv = symm.vectorPostMultiply(eigenv);\r
 157         System.out.println(" --- symmat * eigenv / lambda --- ");\r
 158 \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
 162             }\r
 163         }\r
 164 \r
 165         System.out.println();\r
 166     }\r
 167 \r
 168     public void run() {\r
 169         Matrix mt = m.transpose();\r
 170 \r
 171         //    System.out.println(" --- OrigT * Orig ---- ");\r
 172         eigenvector = mt.preMultiply(m2);\r
 173 \r
 174         //  eigenvector.print(System.out);\r
 175         symm = eigenvector.copy();\r
 176 \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
 183         //f.show();\r
 184         //symm.print(taps);\r
 185         long tstart = System.currentTimeMillis();\r
 186         eigenvector.tred();\r
 187 \r
 188         long tend = System.currentTimeMillis();\r
 189 \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
 194         //taps.println();\r
 195         //taps.println(" --- E vector ---");\r
 196         //    eigenvector.printE(taps);\r
 197         //taps.println();\r
 198         // Now produce the diagonalization matrix\r
 199         tstart = System.currentTimeMillis();\r
 200         eigenvector.tqli();\r
 201         tend = System.currentTimeMillis();\r
 202 \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
 210         // taps.println();\r
 211         // }\r
 212         //  taps.println();\r
 213         //  taps.println("Transformed sequences = ");\r
 214         // Matrix trans =  m.preMultiply(eigenvector);\r
 215         //  trans.print(System.out);\r
 216     }\r
 217 }\r