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.Console;
  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,
  64           SimilarityParamsI options)
  65   {
  66     this.seqs = sequences;
  67     this.scoreModel = sm;
  68     this.similarityParams = options;
  69   }
  70
  71   /**
  72    * Returns Eigenvalue
  73    *
  74    * @param i
  75    *          Index of diagonal within matrix
  76    *
  77    * @return Returns value of diagonal from matrix
  78    */
  79   public double getEigenvalue(int i)
  80   {
  81     return eigenMatrix.getD()[i];
  82   }
  83
  84   /**
  85    * DOCUMENT ME!
  86    *
  87    * @param l
  88    *          DOCUMENT ME!
  89    * @param n
  90    *          DOCUMENT ME!
  91    * @param mm
  92    *          DOCUMENT ME!
  93    * @param factor
  94    *          DOCUMENT ME!
  95    *
  96    * @return DOCUMENT ME!
  97    */
  98   public Point[] getComponents(int l, int n, int mm, float factor)
  99   {
 100     Point[] out = new Point[getHeight()];
 101
 102     for (int i = 0; i < getHeight(); i++)
 103     {
 104       float x = (float) component(i, l) * factor;
 105       float y = (float) component(i, n) * factor;
 106       float z = (float) component(i, mm) * factor;
 107       out[i] = new Point(x, y, z);
 108     }
 109
 110     return out;
 111   }
 112
 113   /**
 114    * DOCUMENT ME!
 115    *
 116    * @param n
 117    *          DOCUMENT ME!
 118    *
 119    * @return DOCUMENT ME!
 120    */
 121   public double[] component(int n)
 122   {
 123     // n = index of eigenvector
 124     double[] out = new double[getHeight()];
 125
 126     for (int i = 0; i < out.length; i++)
 127     {
 128       out[i] = component(i, n);
 129     }
 130
 131     return out;
 132   }
 133
 134   /**
 135    * DOCUMENT ME!
 136    *
 137    * @param row
 138    *          DOCUMENT ME!
 139    * @param n
 140    *          DOCUMENT ME!
 141    *
 142    * @return DOCUMENT ME!
 143    */
 144   double component(int row, int n)
 145   {
 146     double out = 0.0;
 147
 148     for (int i = 0; i < pairwiseScores.width(); i++)
 149     {
 150       out += (pairwiseScores.getValue(row, i) * eigenMatrix.getValue(i, n));
 151     }
 152
 153     return out / eigenMatrix.getD()[n];
 154   }
 155
 156   /**
 157    * Answers a formatted text report of the PCA calculation results (matrices
 158    * and eigenvalues) suitable for display
 159    *
 160    * @return
 161    */
 162   public String getDetails()
 163   {
 164     StringBuilder sb = new StringBuilder(1024);
 165     sb.append("PCA calculation using ").append(scoreModel.getName())
 166             .append(" sequence similarity matrix\n========\n\n");
 167     PrintStream ps = wrapOutputBuffer(sb);
 168
 169     /*
 170      * pairwise similarity scores
 171      */
 172     sb.append(" --- OrigT * Orig ---- \n");
 173     pairwiseScores.print(ps, "%8.2f");
 174
 175     /*
 176      * tridiagonal matrix, with D and E vectors
 177      */
 178     sb.append(" ---Tridiag transform matrix ---\n");
 179     sb.append(" --- D vector ---\n");
 180     tridiagonal.printD(ps, "%15.4e");
 181     ps.println();
 182     sb.append("--- E vector ---\n");
 183     tridiagonal.printE(ps, "%15.4e");
 184     ps.println();
 185
 186     /*
 187      * eigenvalues matrix, with D vector
 188      */
 189     sb.append(" --- New diagonalization matrix ---\n");
 190     eigenMatrix.print(ps, "%8.2f");
 191     sb.append(" --- Eigenvalues ---\n");
 192     eigenMatrix.printD(ps, "%15.4e");
 193     ps.println();
 194
 195     return sb.toString();
 196   }
 197
 198   /**
 199    * Performs the PCA calculation
 200    */
 201   @Override
 202   public void run()
 203   {
 204     try
 205     {
 206       /*
 207        * sequence pairwise similarity scores
 208        */
 209       pairwiseScores = scoreModel.findSimilarities(seqs, similarityParams);
 210
 211       /*
 212        * tridiagonal matrix
 213        */
 214       tridiagonal = pairwiseScores.copy();
 215       tridiagonal.tred();
 216
 217       /*
 218        * the diagonalization matrix
 219        */
 220       eigenMatrix = tridiagonal.copy();
 221       eigenMatrix.tqli();
 222     } catch (Exception q)
 223     {
 224       Console.error("Error computing PCA:  " + q.getMessage());
 225       q.printStackTrace();
 226     }
 227   }
 228
 229   /**
 230    * Returns a PrintStream that wraps (appends its output to) the given
 231    * StringBuilder
 232    *
 233    * @param sb
 234    * @return
 235    */
 236   protected PrintStream wrapOutputBuffer(StringBuilder sb)
 237   {
 238     PrintStream ps = new PrintStream(System.out)
 239     {
 240       @Override
 241       public void print(String x)
 242       {
 243         sb.append(x);
 244       }
 245
 246       @Override
 247       public void println()
 248       {
 249         sb.append("\n");
 250       }
 251     };
 252     return ps;
 253   }
 254
 255   /**
 256    * Answers the N dimensions of the NxN PCA matrix. This is the number of
 257    * sequences involved in the pairwise score calculation.
 258    *
 259    * @return
 260    */
 261   public int getHeight()
 262   {
 263     // TODO can any of seqs[] be null?
 264     return pairwiseScores.height();// seqs.getSequences().length;
 265   }
 266
 267   /**
 268    * Answers the sequence pairwise similarity scores which were the first step
 269    * of the PCA calculation
 270    *
 271    * @return
 272    */
 273   public MatrixI getPairwiseScores()
 274   {
 275     return pairwiseScores;
 276   }
 277
 278   public void setPairwiseScores(MatrixI m)
 279   {
 280     pairwiseScores = m;
 281   }
 282
 283   public MatrixI getEigenmatrix()
 284   {
 285     return eigenMatrix;
 286   }
 287
 288   public void setEigenmatrix(MatrixI m)
 289   {
 290     eigenMatrix = m;
 291   }
 292
 293   public MatrixI getTridiagonal()
 294   {
 295     return tridiagonal;
 296   }
 297
 298   public void setTridiagonal(MatrixI tridiagonal)
 299   {
 300     this.tridiagonal = tridiagonal;
 301   }
 302 }