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