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 \r
  26 /**\r
  27  * Performs Principal Component Analysis on given sequences\r
  28  *\r
  29  * @author $author$\r
  30  * @version $Revision$\r
  31  */\r
  32 public class PCA implements Runnable\r
  33 {\r
  34     Matrix m;\r
  35     Matrix symm;\r
  36     Matrix m2;\r
  37     double[] eigenvalue;\r
  38     Matrix eigenvector;\r
  39 \r
  40 \r
  41     /**\r
  42      * Creates a new PCA object.\r
  43      *\r
  44      * @param s Set of sequences to perform PCA on\r
  45      */\r
  46     public PCA(SequenceI[] s)\r
  47     {\r
  48 \r
  49         BinarySequence[] bs = new BinarySequence[s.length];\r
  50         int ii = 0;\r
  51 \r
  52         while ((ii < s.length) && (s[ii] != null))\r
  53         {\r
  54             bs[ii] = new BinarySequence(s[ii]);\r
  55             bs[ii].encode();\r
  56             ii++;\r
  57         }\r
  58 \r
  59         BinarySequence[] bs2 = new BinarySequence[s.length];\r
  60         ii = 0;\r
  61 \r
  62         while ((ii < s.length) && (s[ii] != null))\r
  63         {\r
  64             bs2[ii] = new BinarySequence(s[ii]);\r
  65             bs2[ii].blosumEncode();\r
  66             ii++;\r
  67         }\r
  68 \r
  69         //System.out.println("Created binary encoding");\r
  70         //printMemory(rt);\r
  71         int count = 0;\r
  72 \r
  73         while ((count < bs.length) && (bs[count] != null))\r
  74         {\r
  75             count++;\r
  76         }\r
  77 \r
  78         double[][] seqmat = new double[count][bs[0].getDBinary().length];\r
  79         double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];\r
  80         int i = 0;\r
  81 \r
  82         while (i < count)\r
  83         {\r
  84             seqmat[i] = bs[i].getDBinary();\r
  85             seqmat2[i] = bs2[i].getDBinary();\r
  86             i++;\r
  87         }\r
  88 \r
  89         //System.out.println("Created array");\r
  90         //printMemory(rt);\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
  94 \r
  95       }\r
  96 \r
  97       /**\r
  98        * Returns the matrix used in PCA calculation\r
  99        *\r
 100        * @return java.math.Matrix object\r
 101        */\r
 102 \r
 103       public Matrix getM()\r
 104       {\r
 105         return m;\r
 106       }\r
 107 \r
 108     /**\r
 109      * Returns Eigenvalue\r
 110      *\r
 111      * @param i Index of diagonal within matrix\r
 112      *\r
 113      * @return Returns value of diagonal from matrix\r
 114      */\r
 115     public double getEigenvalue(int i)\r
 116     {\r
 117         return eigenvector.d[i];\r
 118     }\r
 119 \r
 120      /**\r
 121      * DOCUMENT ME!\r
 122      *\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
 127      *\r
 128      * @return DOCUMENT ME!\r
 129      */\r
 130     public float[][] getComponents(int l, int n, int mm, float factor)\r
 131     {\r
 132         float[][] out = new float[m.rows][3];\r
 133 \r
 134         for (int i = 0; i < m.rows; i++)\r
 135         {\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
 139         }\r
 140 \r
 141         return out;\r
 142     }\r
 143 \r
 144     /**\r
 145      * DOCUMENT ME!\r
 146      *\r
 147      * @param n DOCUMENT ME!\r
 148      *\r
 149      * @return DOCUMENT ME!\r
 150      */\r
 151     public double[] component(int n)\r
 152     {\r
 153         // n = index of eigenvector\r
 154         double[] out = new double[m.rows];\r
 155 \r
 156         for (int i = 0; i < m.rows; i++)\r
 157         {\r
 158             out[i] = component(i, n);\r
 159         }\r
 160 \r
 161         return out;\r
 162     }\r
 163 \r
 164     /**\r
 165      * DOCUMENT ME!\r
 166      *\r
 167      * @param row DOCUMENT ME!\r
 168      * @param n DOCUMENT ME!\r
 169      *\r
 170      * @return DOCUMENT ME!\r
 171      */\r
 172     double component(int row, int n)\r
 173     {\r
 174         double out = 0.0;\r
 175 \r
 176         for (int i = 0; i < symm.cols; i++)\r
 177         {\r
 178             out += (symm.value[row][i] * eigenvector.value[i][n]);\r
 179         }\r
 180 \r
 181         return out / eigenvector.d[n];\r
 182     }\r
 183 \r
 184 \r
 185     /**\r
 186      * DOCUMENT ME!\r
 187      */\r
 188     public void run()\r
 189     {\r
 190         Matrix mt = m.transpose();\r
 191 \r
 192         //    System.out.println(" --- OrigT * Orig ---- ");\r
 193         eigenvector = mt.preMultiply(m2);\r
 194 \r
 195         //  eigenvector.print(System.out);\r
 196         symm = eigenvector.copy();\r
 197 \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
 204         //f.show();\r
 205         //symm.print(taps);\r
 206         long tstart = System.currentTimeMillis();\r
 207         eigenvector.tred();\r
 208 \r
 209         long tend = System.currentTimeMillis();\r
 210 \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
 215         //taps.println();\r
 216         //taps.println(" --- E vector ---");\r
 217         //    eigenvector.printE(taps);\r
 218         //taps.println();\r
 219         // Now produce the diagonalization matrix\r
 220         tstart = System.currentTimeMillis();\r
 221         eigenvector.tqli();\r
 222         tend = System.currentTimeMillis();\r
 223 \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
 231         // taps.println();\r
 232         // }\r
 233         //  taps.println();\r
 234         //  taps.println("Transformed sequences = ");\r
 235         // Matrix trans =  m.preMultiply(eigenvector);\r
 236         //  trans.print(System.out);\r
 237     }\r
 238 }\r