src/jalview/analysis/PCA.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer (Version 2.4)
   3  * Copyright (C) 2008 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle
   4  *
   5  * This program is free software; you can redistribute it and/or
   6  * modify it under the terms of the GNU General Public License
   7  * as published by the Free Software Foundation; either version 2
   8  * of the License, or (at your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13  * GNU General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  18  */
  19 package jalview.analysis;
  20
  21 import java.io.*;
  22
  23 import jalview.datamodel.*;
  24 import jalview.math.*;
  25
  26 /**
  27  * Performs Principal Component Analysis on given sequences
  28  *
  29  * @author $author$
  30  * @version $Revision$
  31  */
  32 public class PCA
  33     implements Runnable
  34 {
  35   Matrix m;
  36   Matrix symm;
  37   Matrix m2;
  38   double[] eigenvalue;
  39   Matrix eigenvector;
  40   StringBuffer details = new StringBuffer();
  41
  42   /**
  43    * Creates a new PCA object.
  44    *
  45    * @param s Set of sequences to perform PCA on
  46    */
  47   public PCA(String[] s)
  48   {
  49
  50     BinarySequence[] bs = new BinarySequence[s.length];
  51     int ii = 0;
  52
  53     while ( (ii < s.length) && (s[ii] != null))
  54     {
  55       bs[ii] = new BinarySequence(s[ii]);
  56       bs[ii].encode();
  57       ii++;
  58     }
  59
  60     BinarySequence[] bs2 = new BinarySequence[s.length];
  61     ii = 0;
  62
  63     while ( (ii < s.length) && (s[ii] != null))
  64     {
  65       bs2[ii] = new BinarySequence(s[ii]);
  66       bs2[ii].blosumEncode();
  67       ii++;
  68     }
  69
  70     //System.out.println("Created binary encoding");
  71     //printMemory(rt);
  72     int count = 0;
  73
  74     while ( (count < bs.length) && (bs[count] != null))
  75     {
  76       count++;
  77     }
  78
  79     double[][] seqmat = new double[count][bs[0].getDBinary().length];
  80     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];
  81     int i = 0;
  82
  83     while (i < count)
  84     {
  85       seqmat[i] = bs[i].getDBinary();
  86       seqmat2[i] = bs2[i].getDBinary();
  87       i++;
  88     }
  89
  90     //System.out.println("Created array");
  91     //printMemory(rt);
  92     //    System.out.println(" --- Original matrix ---- ");
  93     m = new Matrix(seqmat, count, bs[0].getDBinary().length);
  94     m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length);
  95
  96   }
  97
  98   /**
  99    * Returns the matrix used in PCA calculation
 100    *
 101    * @return java.math.Matrix object
 102    */
 103
 104   public Matrix getM()
 105   {
 106     return m;
 107   }
 108
 109   /**
 110    * Returns Eigenvalue
 111    *
 112    * @param i Index of diagonal within matrix
 113    *
 114    * @return Returns value of diagonal from matrix
 115    */
 116   public double getEigenvalue(int i)
 117   {
 118     return eigenvector.d[i];
 119   }
 120
 121   /**
 122    * DOCUMENT ME!
 123    *
 124    * @param l DOCUMENT ME!
 125    * @param n DOCUMENT ME!
 126    * @param mm DOCUMENT ME!
 127    * @param factor DOCUMENT ME!
 128    *
 129    * @return DOCUMENT ME!
 130    */
 131   public float[][] getComponents(int l, int n, int mm, float factor)
 132   {
 133     float[][] out = new float[m.rows][3];
 134
 135     for (int i = 0; i < m.rows; i++)
 136     {
 137       out[i][0] = (float) component(i, l) * factor;
 138       out[i][1] = (float) component(i, n) * factor;
 139       out[i][2] = (float) component(i, mm) * factor;
 140     }
 141
 142     return out;
 143   }
 144
 145   /**
 146    * DOCUMENT ME!
 147    *
 148    * @param n DOCUMENT ME!
 149    *
 150    * @return DOCUMENT ME!
 151    */
 152   public double[] component(int n)
 153   {
 154     // n = index of eigenvector
 155     double[] out = new double[m.rows];
 156
 157     for (int i = 0; i < m.rows; i++)
 158     {
 159       out[i] = component(i, n);
 160     }
 161
 162     return out;
 163   }
 164
 165   /**
 166    * DOCUMENT ME!
 167    *
 168    * @param row DOCUMENT ME!
 169    * @param n DOCUMENT ME!
 170    *
 171    * @return DOCUMENT ME!
 172    */
 173   double component(int row, int n)
 174   {
 175     double out = 0.0;
 176
 177     for (int i = 0; i < symm.cols; i++)
 178     {
 179       out += (symm.value[row][i] * eigenvector.value[i][n]);
 180     }
 181
 182     return out / eigenvector.d[n];
 183   }
 184
 185   public String getDetails()
 186   {
 187     return details.toString();
 188   }
 189
 190   /**
 191    * DOCUMENT ME!
 192    */
 193   public void run()
 194   {
 195     Matrix mt = m.transpose();
 196
 197     details.append(" --- OrigT * Orig ---- \n");
 198     eigenvector = mt.preMultiply(m2);
 199
 200     PrintStream ps = new PrintStream(System.out)
 201     {
 202       public void print(String x)
 203       {
 204         details.append(x);
 205       }
 206
 207       public void println()
 208       {
 209         details.append("\n");
 210       }
 211     };
 212
 213     eigenvector.print(ps);
 214
 215     symm = eigenvector.copy();
 216
 217     eigenvector.tred();
 218
 219     details.append(" ---Tridiag transform matrix ---\n");
 220     details.append(" --- D vector ---\n");
 221     eigenvector.printD(ps);
 222     ps.println();
 223     details.append("--- E vector ---\n");
 224     eigenvector.printE(ps);
 225     ps.println();
 226
 227     // Now produce the diagonalization matrix
 228     eigenvector.tqli();
 229
 230     details.append(" --- New diagonalization matrix ---\n");
 231     details.append(" --- Eigenvalues ---\n");
 232     eigenvector.printD(ps);
 233     ps.println();
 234     //  taps.println();
 235     //  taps.println("Transformed sequences = ");
 236     // Matrix trans =  m.preMultiply(eigenvector);
 237     //  trans.print(System.out);
 238   }
 239 }