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