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