src/jalview/analysis/PCA.java

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