src/jalview/analysis/PCA.java

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