src/jalview/analysis/PCA.java

   1 package jalview.analysis;\r
   2 \r
   3 import jalview.math.*;\r
   4 import jalview.datamodel.*;\r
   5 import jalview.util.*;\r
   6 \r
   7 import java.awt.*;\r
   8 import java.io.*;\r
   9 \r
  10 public class PCA implements Runnable {\r
  11   Matrix m;\r
  12   Matrix symm;\r
  13   Matrix m2;\r
  14 \r
  15   double[] eigenvalue;\r
  16   Matrix eigenvector;\r
  17 \r
  18   public PCA(Matrix m) {\r
  19     this.m = m;\r
  20   }\r
  21 \r
  22   public PCA(SequenceI[] s) {\r
  23     Runtime rt = Runtime.getRuntime();\r
  24 \r
  25     BinarySequence[] bs = new BinarySequence[s.length];\r
  26     int ii = 0;\r
  27     while (ii < s.length && s[ii] != null) {\r
  28 \r
  29       bs[ii] = new BinarySequence(s[ii]);\r
  30       bs[ii].encode();\r
  31       ii++;\r
  32     }\r
  33 \r
  34     BinarySequence[] bs2 = new BinarySequence[s.length];\r
  35     ii = 0;\r
  36     while (ii < s.length && s[ii] != null) {\r
  37 \r
  38       bs2[ii] = new BinarySequence(s[ii]);\r
  39       bs2[ii].blosumEncode();\r
  40       ii++;\r
  41     }\r
  42 \r
  43 \r
  44     //System.out.println("Created binary encoding");\r
  45     //printMemory(rt);\r
  46 \r
  47     int count=0;\r
  48     while (count < bs.length && bs[count] != null) {\r
  49       count++;\r
  50     }\r
  51     double[][] seqmat = new double[count][bs[0].getDBinary().length];\r
  52     double[][] seqmat2 = new double[count][bs2[0].getDBinary().length];\r
  53     int i=0;\r
  54     while (i < count) {\r
  55       seqmat[i] = bs[i].getDBinary();\r
  56       seqmat2[i] = bs2[i].getDBinary();\r
  57       i++;\r
  58     }\r
  59     //System.out.println("Created array");\r
  60     //printMemory(rt);\r
  61     //    System.out.println(" --- Original matrix ---- ");\r
  62     m = new Matrix(seqmat,count,bs[0].getDBinary().length);\r
  63     m2 = new Matrix(seqmat2,count,bs2[0].getDBinary().length);\r
  64 \r
  65     //System.out.println("Created matrix");\r
  66     printMemory(rt);\r
  67   }\r
  68 \r
  69   public static void printMemory(Runtime rt) {\r
  70     System.out.println("Free memory = " + rt.freeMemory());\r
  71   }\r
  72 \r
  73   public Matrix getM() {\r
  74     return m;\r
  75   }\r
  76 \r
  77   public double[] getEigenvector(int i) {\r
  78     return eigenvector.getColumn(i);\r
  79   }\r
  80 \r
  81   public double getEigenvalue(int i) {\r
  82     return eigenvector.d[i];\r
  83   }\r
  84   public float[][] getComponents(int l, int n, int mm) {\r
  85     return getComponents(l,n,mm,1);\r
  86   }\r
  87   public float[][] getComponents(int l, int n, int mm, float factor) {\r
  88     float[][] out = new float[m.rows][3];\r
  89 \r
  90     for (int i = 0; i < m.rows;i++) {\r
  91       out[i][0] = (float)component(i,l)*factor;\r
  92       out[i][1] = (float)component(i,n)*factor;\r
  93       out[i][2] = (float)component(i,mm)*factor;\r
  94     }\r
  95     return out;\r
  96   }\r
  97 \r
  98   public double[] component(int n) {\r
  99     // n = index of eigenvector\r
 100     double[] out = new double[m.rows];\r
 101 \r
 102     for (int i=0; i < m.rows; i++) {\r
 103       out[i] = component(i,n);\r
 104     }\r
 105     return out;\r
 106   }\r
 107   public double component(int row, int n) {\r
 108     double out = 0.0;\r
 109 \r
 110     for (int i = 0; i < symm.cols; i++) {\r
 111       out += symm.value[row][i] * eigenvector.value[i][n];\r
 112     }\r
 113     return out/eigenvector.d[n];\r
 114   }\r
 115 \r
 116   public void checkEigenvector(int n,PrintStream ps) {\r
 117     ps.println(" --- Eigenvector " + n  + " --- ");\r
 118 \r
 119     double[] eigenv = eigenvector.getColumn(n);\r
 120 \r
 121     for (int i=0; i < eigenv.length;i++) {\r
 122       Format.print(ps,"%15.4f",eigenv[i]);\r
 123     }\r
 124 \r
 125     System.out.println();\r
 126 \r
 127     double[] neigenv = symm.vectorPostMultiply(eigenv);\r
 128     System.out.println(" --- symmat * eigenv / lambda --- ");\r
 129     if (eigenvector.d[n] > 1e-4) {\r
 130       for (int i=0; i < neigenv.length;i++) {\r
 131         Format.print(System.out,"%15.4f",neigenv[i]/eigenvector.d[n]);\r
 132       }\r
 133     }\r
 134     System.out.println();\r
 135   }\r
 136 \r
 137   public void run() {\r
 138     Matrix mt = m.transpose();\r
 139     //    System.out.println(" --- OrigT * Orig ---- ");\r
 140     eigenvector = mt.preMultiply(m2);\r
 141     //  eigenvector.print(System.out);\r
 142     symm = eigenvector.copy();\r
 143 \r
 144     //TextArea ta = new TextArea(25,72);\r
 145     //TextAreaPrintStream taps = new TextAreaPrintStream(System.out,ta);\r
 146     //Frame f = new Frame("PCA output");\r
 147     //f.resize(500,500);\r
 148     //f.setLayout(new BorderLayout());\r
 149     //f.add("Center",ta);\r
 150     //f.show();\r
 151     //symm.print(taps);\r
 152     long tstart = System.currentTimeMillis();\r
 153     eigenvector.tred();\r
 154     long tend = System.currentTimeMillis();\r
 155     //taps.println("Time take for tred = " + (tend-tstart) + "ms");\r
 156     //taps.println(" ---Tridiag transform matrix ---");\r
 157 \r
 158     //taps.println(" --- D vector ---");\r
 159     //eigenvector.printD(taps);\r
 160     //taps.println();\r
 161     //taps.println(" --- E vector ---");\r
 162     //    eigenvector.printE(taps);\r
 163     //taps.println();\r
 164 \r
 165     // Now produce the diagonalization matrix\r
 166     tstart = System.currentTimeMillis();\r
 167     eigenvector.tqli();\r
 168     tend = System.currentTimeMillis();\r
 169     //System.out.println("Time take for tqli = " + (tend-tstart) + " ms");\r
 170 \r
 171     //System.out.println(" --- New diagonalization matrix ---");\r
 172 \r
 173     //System.out.println(" --- Eigenvalues ---");\r
 174     //eigenvector.printD(taps);\r
 175 \r
 176     //System.out.println();\r
 177 \r
 178     // for (int i=0; i < eigenvector.cols; i++) {\r
 179     // checkEigenvector(i,taps);\r
 180     // taps.println();\r
 181     // }\r
 182 \r
 183     //  taps.println();\r
 184     //  taps.println("Transformed sequences = ");\r
 185     // Matrix trans =  m.preMultiply(eigenvector);\r
 186     //  trans.print(System.out);\r
 187   }\r
 188 \r
 189 }\r