X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=src%2Fjalview%2Fanalysis%2FPCA.java;h=3ec7995f9f67aa84c2553fcd154df0b4c934a223;hb=136c0793b90b72b928c4d77dc109dd5c644e00d3;hp=f9a034334a2b78e0be9fd50b7c89ebc20de95a5c;hpb=506d60f0e188723ddc91c26824b41ac7034df3fe;p=jalview.git diff --git a/src/jalview/analysis/PCA.java b/src/jalview/analysis/PCA.java index f9a0343..3ec7995 100755 --- a/src/jalview/analysis/PCA.java +++ b/src/jalview/analysis/PCA.java @@ -1,148 +1,93 @@ /* - * Jalview - A Sequence Alignment Editor and Viewer (Version 2.4) - * Copyright (C) 2008 AM Waterhouse, J Procter, G Barton, M Clamp, S Searle + * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$) + * Copyright (C) $$Year-Rel$$ The Jalview Authors * - * This program is free software; you can redistribute it and/or - * modify it under the terms of the GNU General Public License - * as published by the Free Software Foundation; either version 2 - * of the License, or (at your option) any later version. + * This file is part of Jalview. * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. + * Jalview is free software: you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation, either version 3 + * of the License, or (at your option) any later version. + * + * Jalview is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR + * PURPOSE. See the GNU General Public License for more details. * * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA + * along with Jalview. If not, see . + * The Jalview Authors are detailed in the 'AUTHORS' file. */ package jalview.analysis; -import java.io.*; +import jalview.api.analysis.ScoreModelI; +import jalview.api.analysis.SimilarityParamsI; +import jalview.datamodel.AlignmentView; +import jalview.math.MatrixI; -import jalview.datamodel.*; -import jalview.math.*; +import java.io.PrintStream; /** * Performs Principal Component Analysis on given sequences - * - * @author $author$ - * @version $Revision$ */ public class PCA implements Runnable { - Matrix m; - - Matrix symm; - - Matrix m2; + MatrixI symm; double[] eigenvalue; - Matrix eigenvector; - - StringBuffer details = new StringBuffer(); + MatrixI eigenvector; - /** - * Creates a new PCA object. - * - * @param s - * Set of sequences to perform PCA on - */ - public PCA(String[] s) - { + StringBuilder details = new StringBuilder(1024); - BinarySequence[] bs = new BinarySequence[s.length]; - int ii = 0; + final private AlignmentView seqs; - while ((ii < s.length) && (s[ii] != null)) - { - bs[ii] = new BinarySequence(s[ii]); - bs[ii].encode(); - ii++; - } + private ScoreModelI scoreModel; + + private SimilarityParamsI similarityParams; - BinarySequence[] bs2 = new BinarySequence[s.length]; - ii = 0; - - while ((ii < s.length) && (s[ii] != null)) - { - bs2[ii] = new BinarySequence(s[ii]); - bs2[ii].blosumEncode(); - ii++; - } - - // System.out.println("Created binary encoding"); - // printMemory(rt); - int count = 0; - - while ((count < bs.length) && (bs[count] != null)) - { - count++; - } - - double[][] seqmat = new double[count][bs[0].getDBinary().length]; - double[][] seqmat2 = new double[count][bs2[0].getDBinary().length]; - int i = 0; - - while (i < count) - { - seqmat[i] = bs[i].getDBinary(); - seqmat2[i] = bs2[i].getDBinary(); - i++; - } - - // System.out.println("Created array"); - // printMemory(rt); - // System.out.println(" --- Original matrix ---- "); - m = new Matrix(seqmat, count, bs[0].getDBinary().length); - m2 = new Matrix(seqmat2, count, bs2[0].getDBinary().length); - - } - - /** - * Returns the matrix used in PCA calculation - * - * @return java.math.Matrix object - */ - - public Matrix getM() + public PCA(AlignmentView s, ScoreModelI sm, SimilarityParamsI options) { - return m; + this.seqs = s; + this.similarityParams = options; + this.scoreModel = sm; + + details.append("PCA calculation using " + sm.getName() + + " sequence similarity matrix\n========\n\n"); } /** * Returns Eigenvalue * * @param i - * Index of diagonal within matrix + * Index of diagonal within matrix * * @return Returns value of diagonal from matrix */ public double getEigenvalue(int i) { - return eigenvector.d[i]; + return eigenvector.getD()[i]; } /** * DOCUMENT ME! * * @param l - * DOCUMENT ME! + * DOCUMENT ME! * @param n - * DOCUMENT ME! + * DOCUMENT ME! * @param mm - * DOCUMENT ME! + * DOCUMENT ME! * @param factor - * DOCUMENT ME! + * DOCUMENT ME! * * @return DOCUMENT ME! */ public float[][] getComponents(int l, int n, int mm, float factor) { - float[][] out = new float[m.rows][3]; + float[][] out = new float[getHeight()][3]; - for (int i = 0; i < m.rows; i++) + for (int i = 0; i < getHeight(); i++) { out[i][0] = (float) component(i, l) * factor; out[i][1] = (float) component(i, n) * factor; @@ -156,16 +101,16 @@ public class PCA implements Runnable * DOCUMENT ME! * * @param n - * DOCUMENT ME! + * DOCUMENT ME! * * @return DOCUMENT ME! */ public double[] component(int n) { // n = index of eigenvector - double[] out = new double[m.rows]; + double[] out = new double[getHeight()]; - for (int i = 0; i < m.rows; i++) + for (int i = 0; i < out.length; i++) { out[i] = component(i, n); } @@ -177,9 +122,9 @@ public class PCA implements Runnable * DOCUMENT ME! * * @param row - * DOCUMENT ME! + * DOCUMENT ME! * @param n - * DOCUMENT ME! + * DOCUMENT ME! * * @return DOCUMENT ME! */ @@ -187,12 +132,12 @@ public class PCA implements Runnable { double out = 0.0; - for (int i = 0; i < symm.cols; i++) + for (int i = 0; i < symm.width(); i++) { - out += (symm.value[row][i] * eigenvector.value[i][n]); + out += (symm.getValue(row, i) * eigenvector.getValue(i, n)); } - return out / eigenvector.d[n]; + return out / eigenvector.getD()[n]; } public String getDetails() @@ -203,50 +148,78 @@ public class PCA implements Runnable /** * DOCUMENT ME! */ + @Override public void run() { - Matrix mt = m.transpose(); - - details.append(" --- OrigT * Orig ---- \n"); - eigenvector = mt.preMultiply(m2); - PrintStream ps = new PrintStream(System.out) { + @Override public void print(String x) { details.append(x); } + @Override public void println() { details.append("\n"); } }; - eigenvector.print(ps); + // long now = System.currentTimeMillis(); + try + { + eigenvector = scoreModel.findSimilarities(seqs, similarityParams); + + details.append(" --- OrigT * Orig ---- \n"); + eigenvector.print(ps, "%8.2f"); - symm = eigenvector.copy(); + symm = eigenvector.copy(); - eigenvector.tred(); + eigenvector.tred(); - details.append(" ---Tridiag transform matrix ---\n"); - details.append(" --- D vector ---\n"); - eigenvector.printD(ps); - ps.println(); - details.append("--- E vector ---\n"); - eigenvector.printE(ps); - ps.println(); + details.append(" ---Tridiag transform matrix ---\n"); + details.append(" --- D vector ---\n"); + eigenvector.printD(ps, "%15.4e"); + ps.println(); + details.append("--- E vector ---\n"); + eigenvector.printE(ps, "%15.4e"); + ps.println(); - // Now produce the diagonalization matrix - eigenvector.tqli(); + // Now produce the diagonalization matrix + eigenvector.tqli(); + } catch (Exception q) + { + q.printStackTrace(); + details.append("\n*** Unexpected exception when performing PCA ***\n" + + q.getLocalizedMessage()); + details.append("*** Matrices below may not be fully diagonalised. ***\n"); + } details.append(" --- New diagonalization matrix ---\n"); + eigenvector.print(ps, "%8.2f"); details.append(" --- Eigenvalues ---\n"); - eigenvector.printD(ps); + eigenvector.printD(ps, "%15.4e"); ps.println(); - // taps.println(); - // taps.println("Transformed sequences = "); - // Matrix trans = m.preMultiply(eigenvector); - // trans.print(System.out); + /* + * for (int seq=0;seq