src/jalview/analysis/PCA.java

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