src/jalview/analysis/PCA.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.analysis;
  22
  23 import jalview.api.analysis.ScoreModelI;
  24 import jalview.api.analysis.SimilarityParamsI;
  25 import jalview.bin.Cache;
  26 import jalview.datamodel.AlignmentView;
  27 import jalview.datamodel.Point;
  28 import jalview.math.MatrixI;
  29
  30 import java.io.PrintStream;
  31
  32 /**
  33  * Performs Principal Component Analysis on given sequences
  34  */
  35 public class PCA implements Runnable
  36 {
  37   /*
  38    * inputs
  39    */
  40   final private AlignmentView seqs;
  41
  42   final private ScoreModelI scoreModel;
  43
  44   final private SimilarityParamsI similarityParams;
  45
  46   /*
  47    * outputs
  48    */
  49   private MatrixI pairwiseScores;
  50
  51   private MatrixI tridiagonal;
  52
  53   private MatrixI eigenMatrix;
  54
  55   /**
  56    * Constructor given the sequences to compute for, the similarity model to
  57    * use, and a set of parameters for sequence comparison
  58    *
  59    * @param sequences
  60    * @param sm
  61    * @param options
  62    */
  63   public PCA(AlignmentView sequences, ScoreModelI sm, SimilarityParamsI options)
  64   {
  65     this.seqs = sequences;
  66     this.scoreModel = sm;
  67     this.similarityParams = options;
  68   }
  69
  70   /**
  71    * Returns Eigenvalue
  72    *
  73    * @param i
  74    *          Index of diagonal within matrix
  75    *
  76    * @return Returns value of diagonal from matrix
  77    */
  78   public double getEigenvalue(int i)
  79   {
  80     return eigenMatrix.getD()[i];
  81   }
  82
  83   /**
  84    * DOCUMENT ME!
  85    *
  86    * @param l
  87    *          DOCUMENT ME!
  88    * @param n
  89    *          DOCUMENT ME!
  90    * @param mm
  91    *          DOCUMENT ME!
  92    * @param factor
  93    *          DOCUMENT ME!
  94    *
  95    * @return DOCUMENT ME!
  96    */
  97   public Point[] getComponents(int l, int n, int mm, float factor)
  98   {
  99     Point[] out = new Point[getHeight()];
 100
 101     for (int i = 0; i < getHeight(); i++)
 102     {
 103       float x = (float) component(i, l) * factor;
 104       float y = (float) component(i, n) * factor;
 105       float z = (float) component(i, mm) * factor;
 106       out[i] = new Point(x, y, z);
 107     }
 108
 109     return out;
 110   }
 111
 112   /**
 113    * DOCUMENT ME!
 114    *
 115    * @param n
 116    *          DOCUMENT ME!
 117    *
 118    * @return DOCUMENT ME!
 119    */
 120   public double[] component(int n)
 121   {
 122     // n = index of eigenvector
 123     double[] out = new double[getHeight()];
 124
 125     for (int i = 0; i < out.length; i++)
 126     {
 127       out[i] = component(i, n);
 128     }
 129
 130     return out;
 131   }
 132
 133   /**
 134    * DOCUMENT ME!
 135    *
 136    * @param row
 137    *          DOCUMENT ME!
 138    * @param n
 139    *          DOCUMENT ME!
 140    *
 141    * @return DOCUMENT ME!
 142    */
 143   double component(int row, int n)
 144   {
 145     double out = 0.0;
 146
 147     for (int i = 0; i < pairwiseScores.width(); i++)
 148     {
 149       out += (pairwiseScores.getValue(row, i) * eigenMatrix.getValue(i, n));
 150     }
 151
 152     return out / eigenMatrix.getD()[n];
 153   }
 154
 155   /**
 156    * Answers a formatted text report of the PCA calculation results (matrices
 157    * and eigenvalues) suitable for display
 158    *
 159    * @return
 160    */
 161   public String getDetails()
 162   {
 163     StringBuilder sb = new StringBuilder(1024);
 164     sb.append("PCA calculation using ").append(scoreModel.getName())
 165             .append(" sequence similarity matrix\n========\n\n");
 166     PrintStream ps = wrapOutputBuffer(sb);
 167
 168     /*
 169      * pairwise similarity scores
 170      */
 171     sb.append(" --- OrigT * Orig ---- \n");
 172     pairwiseScores.print(ps, "%8.2f");
 173
 174     /*
 175      * tridiagonal matrix, with D and E vectors
 176      */
 177     sb.append(" ---Tridiag transform matrix ---\n");
 178     sb.append(" --- D vector ---\n");
 179     tridiagonal.printD(ps, "%15.4e");
 180     ps.println();
 181     sb.append("--- E vector ---\n");
 182     tridiagonal.printE(ps, "%15.4e");
 183     ps.println();
 184
 185     /*
 186      * eigenvalues matrix, with D vector
 187      */
 188     sb.append(" --- New diagonalization matrix ---\n");
 189     eigenMatrix.print(ps, "%8.2f");
 190     sb.append(" --- Eigenvalues ---\n");
 191     eigenMatrix.printD(ps, "%15.4e");
 192     ps.println();
 193
 194     return sb.toString();
 195   }
 196
 197   /**
 198    * Performs the PCA calculation
 199    */
 200   @Override
 201   public void run()
 202   {
 203     try
 204     {
 205       /*
 206        * sequence pairwise similarity scores
 207        */
 208       pairwiseScores = scoreModel.findSimilarities(seqs, similarityParams);
 209
 210       /*
 211        * tridiagonal matrix
 212        */
 213       tridiagonal = pairwiseScores.copy();
 214       tridiagonal.tred();
 215
 216       /*
 217        * the diagonalization matrix
 218        */
 219       eigenMatrix = tridiagonal.copy();
 220       eigenMatrix.tqli();
 221     } catch (Exception q)
 222     {
 223       Cache.error("Error computing PCA:  " + q.getMessage());
 224       q.printStackTrace();
 225     }
 226   }
 227
 228   /**
 229    * Returns a PrintStream that wraps (appends its output to) the given
 230    * StringBuilder
 231    *
 232    * @param sb
 233    * @return
 234    */
 235   protected PrintStream wrapOutputBuffer(StringBuilder sb)
 236   {
 237     PrintStream ps = new PrintStream(System.out)
 238     {
 239       @Override
 240       public void print(String x)
 241       {
 242         sb.append(x);
 243       }
 244
 245       @Override
 246       public void println()
 247       {
 248         sb.append("\n");
 249       }
 250     };
 251     return ps;
 252   }
 253
 254   /**
 255    * Answers the N dimensions of the NxN PCA matrix. This is the number of
 256    * sequences involved in the pairwise score calculation.
 257    *
 258    * @return
 259    */
 260   public int getHeight()
 261   {
 262     // TODO can any of seqs[] be null?
 263     return pairwiseScores.height();// seqs.getSequences().length;
 264   }
 265
 266   /**
 267    * Answers the sequence pairwise similarity scores which were the first step
 268    * of the PCA calculation
 269    *
 270    * @return
 271    */
 272   public MatrixI getPairwiseScores()
 273   {
 274     return pairwiseScores;
 275   }
 276
 277   public void setPairwiseScores(MatrixI m)
 278   {
 279     pairwiseScores = m;
 280   }
 281
 282   public MatrixI getEigenmatrix()
 283   {
 284     return eigenMatrix;
 285   }
 286
 287   public void setEigenmatrix(MatrixI m)
 288   {
 289     eigenMatrix = m;
 290   }
 291
 292   public MatrixI getTridiagonal()
 293   {
 294     return tridiagonal;
 295   }
 296
 297   public void setTridiagonal(MatrixI tridiagonal)
 298   {
 299     this.tridiagonal = tridiagonal;
 300   }
 301 }