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.datamodel.AlignmentView;
  26 import jalview.math.MatrixI;
  27
  28 import java.io.PrintStream;
  29
  30 /**
  31  * Performs Principal Component Analysis on given sequences
  32  */
  33 public class PCA implements Runnable
  34 {
  35   MatrixI symm;
  36
  37   double[] eigenvalue;
  38
  39   MatrixI eigenvector;
  40
  41   StringBuilder details = new StringBuilder(1024);
  42
  43   final private AlignmentView seqs;
  44
  45   private ScoreModelI scoreModel;
  46
  47   private SimilarityParamsI similarityParams;
  48
  49   public PCA(AlignmentView s, ScoreModelI sm, SimilarityParamsI options)
  50   {
  51     this.seqs = s;
  52     this.similarityParams = options;
  53     this.scoreModel = sm;
  54
  55     details.append("PCA calculation using " + sm.getName()
  56             + " sequence similarity matrix\n========\n\n");
  57   }
  58
  59   /**
  60    * Returns Eigenvalue
  61    *
  62    * @param i
  63    *          Index of diagonal within matrix
  64    *
  65    * @return Returns value of diagonal from matrix
  66    */
  67   public double getEigenvalue(int i)
  68   {
  69     return eigenvector.getD()[i];
  70   }
  71
  72   /**
  73    * DOCUMENT ME!
  74    *
  75    * @param l
  76    *          DOCUMENT ME!
  77    * @param n
  78    *          DOCUMENT ME!
  79    * @param mm
  80    *          DOCUMENT ME!
  81    * @param factor
  82    *          DOCUMENT ME!
  83    *
  84    * @return DOCUMENT ME!
  85    */
  86   public float[][] getComponents(int l, int n, int mm, float factor)
  87   {
  88     float[][] out = new float[getHeight()][3];
  89
  90     for (int i = 0; i < getHeight(); i++)
  91     {
  92       out[i][0] = (float) component(i, l) * factor;
  93       out[i][1] = (float) component(i, n) * factor;
  94       out[i][2] = (float) component(i, mm) * factor;
  95     }
  96
  97     return out;
  98   }
  99
 100   /**
 101    * DOCUMENT ME!
 102    *
 103    * @param n
 104    *          DOCUMENT ME!
 105    *
 106    * @return DOCUMENT ME!
 107    */
 108   public double[] component(int n)
 109   {
 110     // n = index of eigenvector
 111     double[] out = new double[getHeight()];
 112
 113     for (int i = 0; i < out.length; i++)
 114     {
 115       out[i] = component(i, n);
 116     }
 117
 118     return out;
 119   }
 120
 121   /**
 122    * DOCUMENT ME!
 123    *
 124    * @param row
 125    *          DOCUMENT ME!
 126    * @param n
 127    *          DOCUMENT ME!
 128    *
 129    * @return DOCUMENT ME!
 130    */
 131   double component(int row, int n)
 132   {
 133     double out = 0.0;
 134
 135     for (int i = 0; i < symm.width(); i++)
 136     {
 137       out += (symm.getValue(row, i) * eigenvector.getValue(i, n));
 138     }
 139
 140     return out / eigenvector.getD()[n];
 141   }
 142
 143   public String getDetails()
 144   {
 145     return details.toString();
 146   }
 147
 148   /**
 149    * DOCUMENT ME!
 150    */
 151   @Override
 152   public void run()
 153   {
 154     PrintStream ps = new PrintStream(System.out)
 155     {
 156       @Override
 157       public void print(String x)
 158       {
 159         details.append(x);
 160       }
 161
 162       @Override
 163       public void println()
 164       {
 165         details.append("\n");
 166       }
 167     };
 168
 169     // long now = System.currentTimeMillis();
 170     try
 171     {
 172       eigenvector = scoreModel.findSimilarities(seqs, similarityParams);
 173
 174       details.append(" --- OrigT * Orig ---- \n");
 175       eigenvector.print(ps, "%8.2f");
 176
 177       symm = eigenvector.copy();
 178
 179       eigenvector.tred();
 180
 181       details.append(" ---Tridiag transform matrix ---\n");
 182       details.append(" --- D vector ---\n");
 183       eigenvector.printD(ps, "%15.4e");
 184       ps.println();
 185       details.append("--- E vector ---\n");
 186       eigenvector.printE(ps, "%15.4e");
 187       ps.println();
 188
 189       // Now produce the diagonalization matrix
 190       eigenvector.tqli();
 191     } catch (Exception q)
 192     {
 193       q.printStackTrace();
 194       details.append("\n*** Unexpected exception when performing PCA ***\n"
 195               + q.getLocalizedMessage());
 196       details.append("*** Matrices below may not be fully diagonalised. ***\n");
 197     }
 198
 199     details.append(" --- New diagonalization matrix ---\n");
 200     eigenvector.print(ps, "%8.2f");
 201     details.append(" --- Eigenvalues ---\n");
 202     eigenvector.printD(ps, "%15.4e");
 203     ps.println();
 204     /*
 205      * for (int seq=0;seq<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
 206      * (int ev=0;ev<symm.rows; ev++) {
 207      *
 208      * ps.print(","+component(seq, ev)); } ps.println(); }
 209      */
 210     // System.out.println(("PCA.run() took "
 211     // + (System.currentTimeMillis() - now) + "ms"));
 212   }
 213
 214   /**
 215    * Answers the N dimensions of the NxN PCA matrix. This is the number of
 216    * sequences involved in the pairwise score calculation.
 217    *
 218    * @return
 219    */
 220   public int getHeight()
 221   {
 222     // TODO can any of seqs[] be null?
 223     return seqs.getSequences().length;
 224   }
 225 }