src/jalview/analysis/PCA.java

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