mev-4.0.01/source/org/tigr/microarray/mev/cluster/algorithm/impl/PCA.java

Code
Comments
Other
Rev Date Author Line
2 26 Feb 07 jari 1 /*
2 26 Feb 07 jari 2 Copyright @ 1999-2003, The Institute for Genomic Research (TIGR).
2 26 Feb 07 jari 3 All rights reserved.
2 26 Feb 07 jari 4 */
2 26 Feb 07 jari 5 /*
2 26 Feb 07 jari 6  * $RCSfile: PCA.java,v $
2 26 Feb 07 jari 7  * $Revision: 1.3 $
2 26 Feb 07 jari 8  * $Date: 2005/03/10 15:45:20 $
2 26 Feb 07 jari 9  * $Author: braistedj $
2 26 Feb 07 jari 10  * $State: Exp $
2 26 Feb 07 jari 11  */
2 26 Feb 07 jari 12 package org.tigr.microarray.mev.cluster.algorithm.impl;
2 26 Feb 07 jari 13
2 26 Feb 07 jari 14 import java.util.Vector;
2 26 Feb 07 jari 15
2 26 Feb 07 jari 16 import org.tigr.microarray.mev.cluster.algorithm.AbortException;
2 26 Feb 07 jari 17 import org.tigr.microarray.mev.cluster.algorithm.AbstractAlgorithm;
2 26 Feb 07 jari 18 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmData;
2 26 Feb 07 jari 19 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmEvent;
2 26 Feb 07 jari 20 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmException;
2 26 Feb 07 jari 21 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmParameters;
2 26 Feb 07 jari 22 import org.tigr.util.FloatMatrix;
2 26 Feb 07 jari 23 import org.tigr.util.Maths;
2 26 Feb 07 jari 24 import org.tigr.util.QSort;
2 26 Feb 07 jari 25
2 26 Feb 07 jari 26 public class PCA extends AbstractAlgorithm {
2 26 Feb 07 jari 27     
2 26 Feb 07 jari 28     private boolean stop = false;
2 26 Feb 07 jari 29     private int numNeighbors;    
2 26 Feb 07 jari 30     private int numGenes, numExps;  
2 26 Feb 07 jari 31     private float factor;    
2 26 Feb 07 jari 32     
2 26 Feb 07 jari 33     public AlgorithmData execute(AlgorithmData data) throws AlgorithmException {
2 26 Feb 07 jari 34   FloatMatrix expMatrix = data.getMatrix("experiment");
2 26 Feb 07 jari 35   AlgorithmParameters map = data.getParams();
2 26 Feb 07 jari 36   
2 26 Feb 07 jari 37   int function = map.getInt("distance-function", COVARIANCE);
2 26 Feb 07 jari 38   factor = map.getFloat("distance-factor", 1.0f);
2 26 Feb 07 jari 39   boolean absolute = map.getBoolean("distance-absolute", false);
2 26 Feb 07 jari 40   int mode = map.getInt("pca-mode", 0);
2 26 Feb 07 jari 41   numNeighbors = map.getInt("numNeighbors", 10);  
2 26 Feb 07 jari 42         
2 26 Feb 07 jari 43   numGenes = expMatrix.getRowDimension();
2 26 Feb 07 jari 44   numExps = expMatrix.getColumnDimension();        
2 26 Feb 07 jari 45   
2 26 Feb 07 jari 46   final int numberOfGenes = expMatrix.getRowDimension();
2 26 Feb 07 jari 47   final int numberOfSamples = expMatrix.getColumnDimension();
2 26 Feb 07 jari 48   
2 26 Feb 07 jari 49   AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.PROGRESS_VALUE, 0);
2 26 Feb 07 jari 50   int eventValue = 0;
2 26 Feb 07 jari 51   event.setIntValue(eventValue);
2 26 Feb 07 jari 52   event.setDescription("Calculate covariance matrix\n");
2 26 Feb 07 jari 53   fireValueChanged(event);
2 26 Feb 07 jari 54         /*
2 26 Feb 07 jari 55   FloatMatrix An = new FloatMatrix(numberOfGenes, numberOfSamples);
2 26 Feb 07 jari 56   for (int row=0; row<numberOfGenes; row++) {
2 26 Feb 07 jari 57       for (int column=0; column<numberOfSamples; column++) {
2 26 Feb 07 jari 58     if (!Float.isNaN(expMatrix.get(row, column))) {
2 26 Feb 07 jari 59         An.set(row, column, expMatrix.get(row, column));
2 26 Feb 07 jari 60     } else {
2 26 Feb 07 jari 61         An.set(row, column, 0);
2 26 Feb 07 jari 62     }
2 26 Feb 07 jari 63       }
2 26 Feb 07 jari 64   }
2 26 Feb 07 jari 65         */
2 26 Feb 07 jari 66         FloatMatrix An = imputeKNearestMatrix(expMatrix, numNeighbors);
2 26 Feb 07 jari 67   FloatMatrix matrix = null;
2 26 Feb 07 jari 68   if (mode==0) {
2 26 Feb 07 jari 69       matrix = An;
2 26 Feb 07 jari 70   } else {
2 26 Feb 07 jari 71       matrix = new FloatMatrix(numberOfSamples, numberOfSamples);
2 26 Feb 07 jari 72       for (int column=0; column<numberOfSamples; column++) {
2 26 Feb 07 jari 73     for (int row=0; row<numberOfSamples; row++) {
2 26 Feb 07 jari 74         //matrix.set(row, column, ExperimentUtil.distance(expMatrix, row, column, function, factor, absolute));
2 26 Feb 07 jari 75                     matrix.set(row, column, ExperimentUtil.distance(An, row, column, function, factor, absolute));
2 26 Feb 07 jari 76     }
2 26 Feb 07 jari 77       }
2 26 Feb 07 jari 78   }
2 26 Feb 07 jari 79   
2 26 Feb 07 jari 80   float[][] A = matrix.getArrayCopy();
2 26 Feb 07 jari 81   int m = matrix.getRowDimension();
2 26 Feb 07 jari 82   int n = matrix.getColumnDimension();
2 26 Feb 07 jari 83   int nu = Math.min(m,n);
2 26 Feb 07 jari 84   float[] s = new float [Math.min(m+1,n)];
2 26 Feb 07 jari 85   int[] order = new int [Math.min(m+1,n)];
2 26 Feb 07 jari 86   float[][] U = new float [m][nu];
2 26 Feb 07 jari 87   float[][] V = new float [n][n];
2 26 Feb 07 jari 88   float[] e = new float [n];
2 26 Feb 07 jari 89   float[] work = new float [m];
2 26 Feb 07 jari 90   boolean wantu = true;
2 26 Feb 07 jari 91   boolean wantv = true;
2 26 Feb 07 jari 92   // Reduce A to bidiagonal form, storing the diagonal elements
2 26 Feb 07 jari 93   // in s and the super-diagonal elements in e.
2 26 Feb 07 jari 94   int nct = Math.min(m-1,n);
2 26 Feb 07 jari 95   int nrt = Math.max(0,Math.min(n-2,m));
2 26 Feb 07 jari 96   for (int i=0; i<Math.min(m+1,n); i++) {
2 26 Feb 07 jari 97       order[i]=i;
2 26 Feb 07 jari 98   }
2 26 Feb 07 jari 99   eventValue++;
2 26 Feb 07 jari 100   event.setIntValue(eventValue);
2 26 Feb 07 jari 101   event.setDescription("Reducing A to bidiagonal form\n");
2 26 Feb 07 jari 102   fireValueChanged(event);
2 26 Feb 07 jari 103   int counter = 0;
2 26 Feb 07 jari 104   //int factor=(int)Math.round(Math.max(nct,nrt)/50.0);
2 26 Feb 07 jari 105   for (int k = 0; k < Math.max(nct,nrt); k++) {
2 26 Feb 07 jari 106       counter++;
2 26 Feb 07 jari 107       if (k < nct) {
2 26 Feb 07 jari 108     // Compute the transformation for the k-th column and
2 26 Feb 07 jari 109     // place the k-th diagonal in s[k].
2 26 Feb 07 jari 110     // Compute 2-norm of k-th column without under/overflow.
2 26 Feb 07 jari 111     s[k] = 0;
2 26 Feb 07 jari 112     for (int i = k; i < m; i++) {
2 26 Feb 07 jari 113         s[k] = Maths.hypot(s[k],A[i][k]);
2 26 Feb 07 jari 114     }
2 26 Feb 07 jari 115     if (s[k] != 0.0) {
2 26 Feb 07 jari 116         if (A[k][k] < 0.0) {
2 26 Feb 07 jari 117       s[k] = -s[k];
2 26 Feb 07 jari 118         }
2 26 Feb 07 jari 119         for (int i = k; i < m; i++) {
2 26 Feb 07 jari 120       A[i][k] /= s[k];
2 26 Feb 07 jari 121         }
2 26 Feb 07 jari 122         A[k][k] += 1.0;
2 26 Feb 07 jari 123     }
2 26 Feb 07 jari 124     s[k] = -s[k];
2 26 Feb 07 jari 125       }
2 26 Feb 07 jari 126       for (int j = k+1; j < n; j++) {
2 26 Feb 07 jari 127     if ((k < nct) & (s[k] != 0.0)) {
2 26 Feb 07 jari 128         // Apply the transformation.
2 26 Feb 07 jari 129         float t = 0;
2 26 Feb 07 jari 130         for (int i = k; i < m; i++) {
2 26 Feb 07 jari 131       t += A[i][k]*A[i][j];
2 26 Feb 07 jari 132         }
2 26 Feb 07 jari 133         t = -t/A[k][k];
2 26 Feb 07 jari 134         for (int i = k; i < m; i++) {
2 26 Feb 07 jari 135       A[i][j] += t*A[i][k];
2 26 Feb 07 jari 136         }
2 26 Feb 07 jari 137     }
2 26 Feb 07 jari 138     // Place the k-th row of A into e for the
2 26 Feb 07 jari 139     // subsequent calculation of the row transformation.
2 26 Feb 07 jari 140     e[j] = A[k][j];
2 26 Feb 07 jari 141       }
2 26 Feb 07 jari 142       if (wantu & (k < nct)) {
2 26 Feb 07 jari 143     // Place the transformation in U for subsequent back
2 26 Feb 07 jari 144     // multiplication.
2 26 Feb 07 jari 145     for (int i = k; i < m; i++) {
2 26 Feb 07 jari 146         U[i][k] = A[i][k];
2 26 Feb 07 jari 147     }
2 26 Feb 07 jari 148       }
2 26 Feb 07 jari 149       if (k < nrt) {
2 26 Feb 07 jari 150     // Compute the k-th row transformation and place the
2 26 Feb 07 jari 151     // k-th super-diagonal in e[k].
2 26 Feb 07 jari 152     // Compute 2-norm without under/overflow.
2 26 Feb 07 jari 153     e[k] = 0;
2 26 Feb 07 jari 154     for (int i = k+1; i < n; i++) {
2 26 Feb 07 jari 155         e[k] = Maths.hypot(e[k],e[i]);
2 26 Feb 07 jari 156     }
2 26 Feb 07 jari 157     if (e[k] != 0.0) {
2 26 Feb 07 jari 158         if (e[k+1] < 0.0) {
2 26 Feb 07 jari 159       e[k] = -e[k];
2 26 Feb 07 jari 160         }
2 26 Feb 07 jari 161         for (int i = k+1; i < n; i++) {
2 26 Feb 07 jari 162       e[i] /= e[k];
2 26 Feb 07 jari 163         }
2 26 Feb 07 jari 164         e[k+1] += 1.0;
2 26 Feb 07 jari 165     }
2 26 Feb 07 jari 166     e[k] = -e[k];
2 26 Feb 07 jari 167     if ((k+1 < m) & (e[k] != 0.0)) {
2 26 Feb 07 jari 168         // Apply the transformation.
2 26 Feb 07 jari 169         for (int i = k+1; i < m; i++) {
2 26 Feb 07 jari 170       work[i] = 0.0f;
2 26 Feb 07 jari 171         }
2 26 Feb 07 jari 172         for (int j = k+1; j < n; j++) {
2 26 Feb 07 jari 173       for (int i = k+1; i < m; i++) {
2 26 Feb 07 jari 174           work[i] += e[j]*A[i][j];
2 26 Feb 07 jari 175       }
2 26 Feb 07 jari 176         }
2 26 Feb 07 jari 177         for (int j = k+1; j < n; j++) {
2 26 Feb 07 jari 178       float t = -e[j]/e[k+1];
2 26 Feb 07 jari 179       for (int i = k+1; i < m; i++) {
2 26 Feb 07 jari 180           A[i][j] += t*work[i];
2 26 Feb 07 jari 181       }
2 26 Feb 07 jari 182         }
2 26 Feb 07 jari 183     }
2 26 Feb 07 jari 184     if (wantv) {
2 26 Feb 07 jari 185         // Place the transformation in V for subsequent
2 26 Feb 07 jari 186         // back multiplication.
2 26 Feb 07 jari 187         for (int i = k+1; i < n; i++) {
2 26 Feb 07 jari 188       V[i][k] = e[i];
2 26 Feb 07 jari 189         }
2 26 Feb 07 jari 190     }
2 26 Feb 07 jari 191       }
2 26 Feb 07 jari 192   }
2 26 Feb 07 jari 193   // Set up the final bidiagonal matrix or order p.
2 26 Feb 07 jari 194   int p = Math.min(n,m+1);
2 26 Feb 07 jari 195   if (nct < n) {
2 26 Feb 07 jari 196       s[nct] = A[nct][nct];
2 26 Feb 07 jari 197   }
2 26 Feb 07 jari 198   if (m < p) {
2 26 Feb 07 jari 199       s[p-1] = 0.0f;
2 26 Feb 07 jari 200   }
2 26 Feb 07 jari 201   if (nrt+1 < p) {
2 26 Feb 07 jari 202       e[nrt] = A[nrt][p-1];
2 26 Feb 07 jari 203   }
2 26 Feb 07 jari 204   e[p-1] = 0.0f;
2 26 Feb 07 jari 205   // If required, generate U.
2 26 Feb 07 jari 206   if (wantu) {
2 26 Feb 07 jari 207       event.setDescription("Generating Matrix U\n");
2 26 Feb 07 jari 208       eventValue++;
2 26 Feb 07 jari 209       event.setIntValue(eventValue);
2 26 Feb 07 jari 210       fireValueChanged(event);
2 26 Feb 07 jari 211       for (int j = nct; j < nu; j++) {
2 26 Feb 07 jari 212     for (int i = 0; i < m; i++) {
2 26 Feb 07 jari 213         U[i][j] = 0.0f;
2 26 Feb 07 jari 214     }
2 26 Feb 07 jari 215     U[j][j] = 1.0f;
2 26 Feb 07 jari 216       }
2 26 Feb 07 jari 217       for (int k = nct-1; k >= 0; k--) {
2 26 Feb 07 jari 218     if (s[k] != 0.0) {
2 26 Feb 07 jari 219         for (int j = k+1; j < nu; j++) {
2 26 Feb 07 jari 220       float t = 0;
2 26 Feb 07 jari 221       for (int i = k; i < m; i++) {
2 26 Feb 07 jari 222           t += U[i][k]*U[i][j];
2 26 Feb 07 jari 223       }
2 26 Feb 07 jari 224       t = -t/U[k][k];
2 26 Feb 07 jari 225       for (int i = k; i < m; i++) {
2 26 Feb 07 jari 226           U[i][j] += t*U[i][k];
2 26 Feb 07 jari 227       }
2 26 Feb 07 jari 228         }
2 26 Feb 07 jari 229         for (int i = k; i < m; i++) {
2 26 Feb 07 jari 230       U[i][k] = -U[i][k];
2 26 Feb 07 jari 231         }
2 26 Feb 07 jari 232         U[k][k] = 1.0f + U[k][k];
2 26 Feb 07 jari 233         for (int i = 0; i < k-1; i++) {
2 26 Feb 07 jari 234       U[i][k] = 0.0f;
2 26 Feb 07 jari 235         }
2 26 Feb 07 jari 236     } else {
2 26 Feb 07 jari 237         for (int i = 0; i < m; i++) {
2 26 Feb 07 jari 238       U[i][k] = 0.0f;
2 26 Feb 07 jari 239         }
2 26 Feb 07 jari 240         U[k][k] = 1.0f;
2 26 Feb 07 jari 241     }
2 26 Feb 07 jari 242       }
2 26 Feb 07 jari 243   }
2 26 Feb 07 jari 244   // If required, generate V.
2 26 Feb 07 jari 245   if (wantv) {
2 26 Feb 07 jari 246       event.setDescription("Generating Matrix V\n");
2 26 Feb 07 jari 247       eventValue++;
2 26 Feb 07 jari 248       event.setIntValue(eventValue);
2 26 Feb 07 jari 249       fireValueChanged(event);
2 26 Feb 07 jari 250       for (int k = n-1; k >= 0; k--) {
2 26 Feb 07 jari 251     if ((k < nrt) & (e[k] != 0.0)) {
2 26 Feb 07 jari 252         for (int j = k+1; j < nu; j++) {
2 26 Feb 07 jari 253       float t = 0;
2 26 Feb 07 jari 254       for (int i = k+1; i < n; i++) {
2 26 Feb 07 jari 255           t += V[i][k]*V[i][j];
2 26 Feb 07 jari 256       }
2 26 Feb 07 jari 257       t = -t/V[k+1][k];
2 26 Feb 07 jari 258       for (int i = k+1; i < n; i++) {
2 26 Feb 07 jari 259           V[i][j] += t*V[i][k];
2 26 Feb 07 jari 260       }
2 26 Feb 07 jari 261         }
2 26 Feb 07 jari 262     }
2 26 Feb 07 jari 263     for (int i = 0; i < n; i++) {
2 26 Feb 07 jari 264         V[i][k] = 0.0f;
2 26 Feb 07 jari 265     }
2 26 Feb 07 jari 266     V[k][k] = 1.0f;
2 26 Feb 07 jari 267       }
2 26 Feb 07 jari 268   }
2 26 Feb 07 jari 269   // Main iteration loop for the singular values.
2 26 Feb 07 jari 270   int pp = p-1;
2 26 Feb 07 jari 271   int iter = 0;
2 26 Feb 07 jari 272   float eps = (float)Math.pow(2.0,-52.0);
2 26 Feb 07 jari 273   event.setDescription("Main iteration loop started...\n");
2 26 Feb 07 jari 274   eventValue++;
2 26 Feb 07 jari 275   event.setIntValue(eventValue);
2 26 Feb 07 jari 276   fireValueChanged(event);
2 26 Feb 07 jari 277   counter=0;
2 26 Feb 07 jari 278   while (p > 0) {
2 26 Feb 07 jari 279       counter++;
2 26 Feb 07 jari 280       int k,kase;
2 26 Feb 07 jari 281       if (counter==240) {
2 26 Feb 07 jari 282     if (stop) {
2 26 Feb 07 jari 283         throw new AbortException();
2 26 Feb 07 jari 284     }
2 26 Feb 07 jari 285     event.setDescription("Main iteration loop.\n");
2 26 Feb 07 jari 286     eventValue++;
2 26 Feb 07 jari 287     event.setIntValue(eventValue);
2 26 Feb 07 jari 288     fireValueChanged(event);
2 26 Feb 07 jari 289     counter=0;
2 26 Feb 07 jari 290       }
2 26 Feb 07 jari 291       // Here is where a test for too many iterations would go.
2 26 Feb 07 jari 292       // This section of the program inspects for
2 26 Feb 07 jari 293       // negligible elements in the s and e arrays.  On
2 26 Feb 07 jari 294       // completion the variables kase and k are set as follows.
2 26 Feb 07 jari 295       // kase = 1     if s(p) and e[k-1] are negligible and k<p
2 26 Feb 07 jari 296       // kase = 2     if s(k) is negligible and k<p
2 26 Feb 07 jari 297       // kase = 3     if e[k-1] is negligible, k<p, and
2 26 Feb 07 jari 298       //              s(k), ..., s(p) are not negligible (qr step).
2 26 Feb 07 jari 299       // kase = 4     if e(p-1) is negligible (convergence).
2 26 Feb 07 jari 300       for (k = p-2; k >= -1; k--) {
2 26 Feb 07 jari 301     if (k == -1) {
2 26 Feb 07 jari 302         break;
2 26 Feb 07 jari 303     }
2 26 Feb 07 jari 304     if (Math.abs(e[k]) <= eps*(Math.abs(s[k]) + Math.abs(s[k+1]))) {
2 26 Feb 07 jari 305         e[k] = 0.0f;
2 26 Feb 07 jari 306         break;
2 26 Feb 07 jari 307     }
2 26 Feb 07 jari 308       }
2 26 Feb 07 jari 309       if (k == p-2) {
2 26 Feb 07 jari 310     kase = 4;
2 26 Feb 07 jari 311       } else {
2 26 Feb 07 jari 312     int ks;
2 26 Feb 07 jari 313     for (ks = p-1; ks >= k; ks--) {
2 26 Feb 07 jari 314         if (ks == k) {
2 26 Feb 07 jari 315       break;
2 26 Feb 07 jari 316         }
2 26 Feb 07 jari 317         float t = (float)((ks != p ? Math.abs(e[ks]) : 0.) +
2 26 Feb 07 jari 318         (ks != k+1 ? Math.abs(e[ks-1]) : 0.));
2 26 Feb 07 jari 319         if (Math.abs(s[ks]) <= eps*t) {
2 26 Feb 07 jari 320       s[ks] = 0.0f;
2 26 Feb 07 jari 321       break;
2 26 Feb 07 jari 322         }
2 26 Feb 07 jari 323     }
2 26 Feb 07 jari 324     if (ks == k) {
2 26 Feb 07 jari 325         kase = 3;
2 26 Feb 07 jari 326     } else if (ks == p-1) {
2 26 Feb 07 jari 327         kase = 1;
2 26 Feb 07 jari 328     } else {
2 26 Feb 07 jari 329         kase = 2;
2 26 Feb 07 jari 330         k = ks;
2 26 Feb 07 jari 331     }
2 26 Feb 07 jari 332       }
2 26 Feb 07 jari 333       k++;
2 26 Feb 07 jari 334       // Perform the task indicated by kase.
2 26 Feb 07 jari 335       switch (kase) {
2 26 Feb 07 jari 336     // Deflate negligible s(p).
2 26 Feb 07 jari 337     case 1: {
2 26 Feb 07 jari 338         float f = e[p-2];
2 26 Feb 07 jari 339         e[p-2] = 0.0f;
2 26 Feb 07 jari 340         for (int j = p-2; j >= k; j--) {
2 26 Feb 07 jari 341       float t = Maths.hypot(s[j],f);
2 26 Feb 07 jari 342       float cs = s[j]/t;
2 26 Feb 07 jari 343       float sn = f/t;
2 26 Feb 07 jari 344       s[j] = t;
2 26 Feb 07 jari 345       if (j != k) {
2 26 Feb 07 jari 346           f = -sn*e[j-1];
2 26 Feb 07 jari 347           e[j-1] = cs*e[j-1];
2 26 Feb 07 jari 348       }
2 26 Feb 07 jari 349       if (wantv) {
2 26 Feb 07 jari 350           for (int i = 0; i < n; i++) {
2 26 Feb 07 jari 351         t = cs*V[i][j] + sn*V[i][p-1];
2 26 Feb 07 jari 352         V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1];
2 26 Feb 07 jari 353         V[i][j] = t;
2 26 Feb 07 jari 354           }
2 26 Feb 07 jari 355       }
2 26 Feb 07 jari 356         }
2 26 Feb 07 jari 357     }
2 26 Feb 07 jari 358     break;
2 26 Feb 07 jari 359     // Split at negligible s(k).
2 26 Feb 07 jari 360     case 2: {
2 26 Feb 07 jari 361         float f = e[k-1];
2 26 Feb 07 jari 362         e[k-1] = 0.0f;
2 26 Feb 07 jari 363         for (int j = k; j < p; j++) {
2 26 Feb 07 jari 364       float t = Maths.hypot(s[j],f);
2 26 Feb 07 jari 365       float cs = s[j]/t;
2 26 Feb 07 jari 366       float sn = f/t;
2 26 Feb 07 jari 367       s[j] = t;
2 26 Feb 07 jari 368       f = -sn*e[j];
2 26 Feb 07 jari 369       e[j] = cs*e[j];
2 26 Feb 07 jari 370       if (wantu) {
2 26 Feb 07 jari 371           for (int i = 0; i < m; i++) {
2 26 Feb 07 jari 372         t = cs*U[i][j] + sn*U[i][k-1];
2 26 Feb 07 jari 373         U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1];
2 26 Feb 07 jari 374         U[i][j] = t;
2 26 Feb 07 jari 375           }
2 26 Feb 07 jari 376       }
2 26 Feb 07 jari 377         }
2 26 Feb 07 jari 378     }
2 26 Feb 07 jari 379     break;
2 26 Feb 07 jari 380     // Perform one qr step.
2 26 Feb 07 jari 381     case 3: {
2 26 Feb 07 jari 382         // Calculate the shift.
2 26 Feb 07 jari 383         float scale = Math.max(Math.max(Math.max(Math.max(
2 26 Feb 07 jari 384         Math.abs(s[p-1]),Math.abs(s[p-2])),Math.abs(e[p-2])),
2 26 Feb 07 jari 385         Math.abs(s[k])),Math.abs(e[k]));
2 26 Feb 07 jari 386         float sp = s[p-1]/scale;
2 26 Feb 07 jari 387         float spm1 = s[p-2]/scale;
2 26 Feb 07 jari 388         float epm1 = e[p-2]/scale;
2 26 Feb 07 jari 389         float sk = s[k]/scale;
2 26 Feb 07 jari 390         float ek = e[k]/scale;
2 26 Feb 07 jari 391         float b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0f;
2 26 Feb 07 jari 392         float c = (sp*epm1)*(sp*epm1);
2 26 Feb 07 jari 393         float shift = 0.0f;
2 26 Feb 07 jari 394         if ((b != 0.0) | (c != 0.0)) {
2 26 Feb 07 jari 395       shift = (float)Math.sqrt(b*b + c);
2 26 Feb 07 jari 396       if (b < 0.0) {
2 26 Feb 07 jari 397           shift = -shift;
2 26 Feb 07 jari 398       }
2 26 Feb 07 jari 399       shift = c/(b + shift);
2 26 Feb 07 jari 400         }
2 26 Feb 07 jari 401         float f = (sk + sp)*(sk - sp) + shift;
2 26 Feb 07 jari 402         float g = sk*ek;
2 26 Feb 07 jari 403         // Chase zeros.
2 26 Feb 07 jari 404         for (int j = k; j < p-1; j++) {
2 26 Feb 07 jari 405       float t = Maths.hypot(f,g);
2 26 Feb 07 jari 406       float cs = f/t;
2 26 Feb 07 jari 407       float sn = g/t;
2 26 Feb 07 jari 408       if (j != k) {
2 26 Feb 07 jari 409           e[j-1] = t;
2 26 Feb 07 jari 410       }
2 26 Feb 07 jari 411       f = cs*s[j] + sn*e[j];
2 26 Feb 07 jari 412       e[j] = cs*e[j] - sn*s[j];
2 26 Feb 07 jari 413       g = sn*s[j+1];
2 26 Feb 07 jari 414       s[j+1] = cs*s[j+1];
2 26 Feb 07 jari 415       if (wantv) {
2 26 Feb 07 jari 416           for (int i = 0; i < n; i++) {
2 26 Feb 07 jari 417         t = cs*V[i][j] + sn*V[i][j+1];
2 26 Feb 07 jari 418         V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1];
2 26 Feb 07 jari 419         V[i][j] = t;
2 26 Feb 07 jari 420           }
2 26 Feb 07 jari 421       }
2 26 Feb 07 jari 422       t = Maths.hypot(f,g);
2 26 Feb 07 jari 423       cs = f/t;
2 26 Feb 07 jari 424       sn = g/t;
2 26 Feb 07 jari 425       s[j] = t;
2 26 Feb 07 jari 426       f = cs*e[j] + sn*s[j+1];
2 26 Feb 07 jari 427       s[j+1] = -sn*e[j] + cs*s[j+1];
2 26 Feb 07 jari 428       g = sn*e[j+1];
2 26 Feb 07 jari 429       e[j+1] = cs*e[j+1];
2 26 Feb 07 jari 430       if (wantu && (j < m-1)) {
2 26 Feb 07 jari 431           for (int i = 0; i < m; i++) {
2 26 Feb 07 jari 432         t = cs*U[i][j] + sn*U[i][j+1];
2 26 Feb 07 jari 433         U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1];
2 26 Feb 07 jari 434         U[i][j] = t;
2 26 Feb 07 jari 435           }
2 26 Feb 07 jari 436       }
2 26 Feb 07 jari 437         }
2 26 Feb 07 jari 438         e[p-2] = f;
2 26 Feb 07 jari 439         iter = iter + 1;
2 26 Feb 07 jari 440     }
2 26 Feb 07 jari 441     break;
2 26 Feb 07 jari 442     // Convergence.
2 26 Feb 07 jari 443     case 4: {
2 26 Feb 07 jari 444         // Make the singular values positive.
2 26 Feb 07 jari 445         if (s[k] <= 0.0) {
2 26 Feb 07 jari 446       s[k] = (s[k] < 0.0 ? -s[k] : 0.0f);
2 26 Feb 07 jari 447       if (wantv) {
2 26 Feb 07 jari 448           for (int i = 0; i <= pp; i++) {
2 26 Feb 07 jari 449         V[i][k] = -V[i][k];
2 26 Feb 07 jari 450           }
2 26 Feb 07 jari 451       }
2 26 Feb 07 jari 452         }
2 26 Feb 07 jari 453         // Order the singular values.
2 26 Feb 07 jari 454         while (k < pp) {
2 26 Feb 07 jari 455       if (s[k] >= s[k+1]) {
2 26 Feb 07 jari 456           break;
2 26 Feb 07 jari 457       }
2 26 Feb 07 jari 458       float t = s[k];
2 26 Feb 07 jari 459       s[k] = s[k+1];
2 26 Feb 07 jari 460       s[k+1] = t;
2 26 Feb 07 jari 461       
2 26 Feb 07 jari 462       int Dummy = order[k];
2 26 Feb 07 jari 463       order[k] = order[k+1];
2 26 Feb 07 jari 464       order[k+1] = Dummy;
2 26 Feb 07 jari 465       
2 26 Feb 07 jari 466       if (wantv && (k < n-1)) {
2 26 Feb 07 jari 467           for (int i = 0; i < n; i++) {
2 26 Feb 07 jari 468         t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t;
2 26 Feb 07 jari 469           }
2 26 Feb 07 jari 470       }
2 26 Feb 07 jari 471       if (wantu && (k < m-1)) {
2 26 Feb 07 jari 472           for (int i = 0; i < m; i++) {
2 26 Feb 07 jari 473         t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t;
2 26 Feb 07 jari 474           }
2 26 Feb 07 jari 475       }
2 26 Feb 07 jari 476       k++;
2 26 Feb 07 jari 477         }
2 26 Feb 07 jari 478         iter = 0;
2 26 Feb 07 jari 479         p--;
2 26 Feb 07 jari 480     }
2 26 Feb 07 jari 481     break;
2 26 Feb 07 jari 482       }
2 26 Feb 07 jari 483   }
2 26 Feb 07 jari 484   
2 26 Feb 07 jari 485   event.setDescription("End SVD calculation.\n");
2 26 Feb 07 jari 486   eventValue++;
2 26 Feb 07 jari 487   event.setIntValue(eventValue);
2 26 Feb 07 jari 488   fireValueChanged(event);
2 26 Feb 07 jari 489   
2 26 Feb 07 jari 490   // T-matrix
2 26 Feb 07 jari 491   FloatMatrix T = new FloatMatrix(U, m, Math.min(m+1, n));
2 26 Feb 07 jari 492   // V-matrix
2 26 Feb 07 jari 493   FloatMatrix Vm = new FloatMatrix(V, n, n);
2 26 Feb 07 jari 494   // S-matrix
2 26 Feb 07 jari 495   FloatMatrix S = new FloatMatrix(n, n);
2 26 Feb 07 jari 496   float[][] array = S.getArray();
2 26 Feb 07 jari 497   for (int i = 0; i < n; i++) {
2 26 Feb 07 jari 498       for (int j = 0; j < n; j++) {
2 26 Feb 07 jari 499     array[i][j] = 0.0f;
2 26 Feb 07 jari 500       }
2 26 Feb 07 jari 501       array[i][i] = s[i];
2 26 Feb 07 jari 502   }
2 26 Feb 07 jari 503   
2 26 Feb 07 jari 504   
2 26 Feb 07 jari 505   AlgorithmData result = new AlgorithmData();
2 26 Feb 07 jari 506   switch (mode) {
2 26 Feb 07 jari 507       case 0:
2 26 Feb 07 jari 508     break;
2 26 Feb 07 jari 509       case 1:
2 26 Feb 07 jari 510     result.addMatrix("U", An.times(T));
2 26 Feb 07 jari 511     break;
2 26 Feb 07 jari 512       case 2:
2 26 Feb 07 jari 513     result.addMatrix("U", An.transpose().times(T));
2 26 Feb 07 jari 514     break;
2 26 Feb 07 jari 515       case 3:
2 26 Feb 07 jari 516     FloatMatrix Q = T.copy();
2 26 Feb 07 jari 517     FloatMatrix D = S.copy();
2 26 Feb 07 jari 518     final int dim = D.getRowDimension();
2 26 Feb 07 jari 519     for (int i=0;i<dim;i++) {
2 26 Feb 07 jari 520         D.set(i,i,1.0f/(float)Math.sqrt(D.get(i,i)));
2 26 Feb 07 jari 521     }
2 26 Feb 07 jari 522     T = An.times(Q.times(D));
2 26 Feb 07 jari 523     result.addMatrix("U", An.transpose().times(T));
2 26 Feb 07 jari 524     break;
2 26 Feb 07 jari 525       default:;
2 26 Feb 07 jari 526   }
2 26 Feb 07 jari 527   
2 26 Feb 07 jari 528   result.addMatrix("T", T);
2 26 Feb 07 jari 529   result.addMatrix("S", S);
2 26 Feb 07 jari 530   result.addMatrix("V", Vm);
2 26 Feb 07 jari 531   
2 26 Feb 07 jari 532   return result;
2 26 Feb 07 jari 533     }
2 26 Feb 07 jari 534     
2 26 Feb 07 jari 535     public void abort() {
2 26 Feb 07 jari 536   stop = true;
2 26 Feb 07 jari 537     }
2 26 Feb 07 jari 538     
2 26 Feb 07 jari 539     private FloatMatrix imputeKNearestMatrix(FloatMatrix inputMatrix, int k) throws AlgorithmException {
2 26 Feb 07 jari 540         int numRows = inputMatrix.getRowDimension();
2 26 Feb 07 jari 541         int numCols = inputMatrix.getColumnDimension();
2 26 Feb 07 jari 542         FloatMatrix resultMatrix = new FloatMatrix(numRows, numCols);
2 26 Feb 07 jari 543         
2 26 Feb 07 jari 544         AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, numGenes);
2 26 Feb 07 jari 545         event.setDescription("Imputing missing values");
2 26 Feb 07 jari 546         fireValueChanged(event);
2 26 Feb 07 jari 547         event.setId(AlgorithmEvent.PROGRESS_VALUE);
2 26 Feb 07 jari 548         
2 26 Feb 07 jari 549         for (int i = 0; i < numRows; i++) {
2 26 Feb 07 jari 550             if (stop) {
2 26 Feb 07 jari 551                 throw new AbortException();
2 26 Feb 07 jari 552             }
2 26 Feb 07 jari 553             //event.setIntValue(i);
2 26 Feb 07 jari 554             //event.setDescription("Imputing missing values: Current gene = " + (i+ 1));
2 26 Feb 07 jari 555             //fireValueChanged(event);
2 26 Feb 07 jari 556             
2 26 Feb 07 jari 557             if (isMissingValues(inputMatrix, i)) {
2 26 Feb 07 jari 558                 //System.out.println("gene " + i + " is missing values");
2 26 Feb 07 jari 559                 Vector nonMissingExpts = new Vector();
2 26 Feb 07 jari 560                 for (int j = 0; j < numCols; j++) {
2 26 Feb 07 jari 561                     if (!Float.isNaN(inputMatrix.A[i][j])) {
2 26 Feb 07 jari 562                         nonMissingExpts.add(new Integer(j));
2 26 Feb 07 jari 563                     }
2 26 Feb 07 jari 564                 }
2 26 Feb 07 jari 565                 Vector geneSubset = getValidGenes(i, inputMatrix, nonMissingExpts); //getValidGenes() returns a Vector of genes that have valid values for all the non-missing expts
2 26 Feb 07 jari 566                 
2 26 Feb 07 jari 567                 //System.out.println(" Valid geneSubset.size() = " + geneSubset.size());
2 26 Feb 07 jari 568                 
2 26 Feb 07 jari 569                 /*
2 26 Feb 07 jari 570                 for (int j = 0; j < geneSubset.size(); j++) {
2 26 Feb 07 jari 571                     System.out.println(((Integer)geneSubset.get(j)).intValue());
2 26 Feb 07 jari 572                 }
2 26 Feb 07 jari 573                  */
2 26 Feb 07 jari 574                 //System.out.println("imputing KNN: current gene = " + i);
2 26 Feb 07 jari 575                 Vector kNearestGenes = getKNearestGenes(i, k, inputMatrix, geneSubset, nonMissingExpts);
2 26 Feb 07 jari 576                 
2 26 Feb 07 jari 577                 /*
2 26 Feb 07 jari 578                 System.out.println("k nearest genes of gene " + i + " : ");
2 26 Feb 07 jari 579                  
2 26 Feb 07 jari 580                 for (int j = 0; j < kNearestGenes.size(); j++) {
2 26 Feb 07 jari 581                     System.out.println("" + ((Integer)kNearestGenes.get(j)).intValue());
2 26 Feb 07 jari 582                 }
2 26 Feb 07 jari 583                  */
2 26 Feb 07 jari 584                 
2 26 Feb 07 jari 585                 //TESTED UPTO HERE -- 12/18/2002***********
2 26 Feb 07 jari 586                 //
2 26 Feb 07 jari 587                 /*
2 26 Feb 07 jari 588                 System.out.print("Gene " + i + " :\t");
2 26 Feb 07 jari 589                 for (int j = 0; j < numCols; j++) {
2 26 Feb 07 jari 590                     System.out.print("" +inputMatrix.A[i][j]);
2 26 Feb 07 jari 591                     System.out.print("\t");
2 26 Feb 07 jari 592                 }
2 26 Feb 07 jari 593                 System.out.println();
2 26 Feb 07 jari 594                 System.out.println("Matrix of k Nearest Genes");
2 26 Feb 07 jari 595                 printSubMatrix(kNearestGenes, inputMatrix);
2 26 Feb 07 jari 596                  */    //
2 26 Feb 07 jari 597                 for (int j = 0; j < numCols; j++) {
2 26 Feb 07 jari 598                     if (!Float.isNaN(inputMatrix.A[i][j])) {
2 26 Feb 07 jari 599                         resultMatrix.A[i][j] = inputMatrix.A[i][j];
2 26 Feb 07 jari 600                     } else {
2 26 Feb 07 jari 601                         
2 26 Feb 07 jari 602                         //System.out.println("just before entering getExptMean(): kNearestGenes.size() = " + kNearestGenes.size());
2 26 Feb 07 jari 603                         
2 26 Feb 07 jari 604                         //resultMatrix.A[i][j] = getExptMean(j, kNearestGenes, inputMatrix);
2 26 Feb 07 jari 605                         resultMatrix.A[i][j] = getExptWeightedMean(i, j, kNearestGenes, inputMatrix);
2 26 Feb 07 jari 606                     }
2 26 Feb 07 jari 607                 }
2 26 Feb 07 jari 608                 //DONE UPTO HERE
2 26 Feb 07 jari 609             }
2 26 Feb 07 jari 610             
2 26 Feb 07 jari 611             else {
2 26 Feb 07 jari 612                 for (int j = 0; j < numCols; j++) {
2 26 Feb 07 jari 613                     resultMatrix.A[i][j] = inputMatrix.A[i][j];
2 26 Feb 07 jari 614                 }
2 26 Feb 07 jari 615             }
2 26 Feb 07 jari 616         }
2 26 Feb 07 jari 617         
2 26 Feb 07 jari 618         return imputeRowAverageMatrix(resultMatrix);
2 26 Feb 07 jari 619     }   
2 26 Feb 07 jari 620     
2 26 Feb 07 jari 621     private FloatMatrix imputeRowAverageMatrix(FloatMatrix inputMatrix) throws AlgorithmException {
2 26 Feb 07 jari 622         int numRows = inputMatrix.getRowDimension();
2 26 Feb 07 jari 623         int numCols = inputMatrix.getColumnDimension();
2 26 Feb 07 jari 624         FloatMatrix resultMatrix = new FloatMatrix(numRows, numCols);
2 26 Feb 07 jari 625         
2 26 Feb 07 jari 626         AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, numGenes);
2 26 Feb 07 jari 627         fireValueChanged(event);
2 26 Feb 07 jari 628         event.setId(AlgorithmEvent.PROGRESS_VALUE);
2 26 Feb 07 jari 629         
2 26 Feb 07 jari 630         for (int i = 0; i < numRows; i++) {
2 26 Feb 07 jari 631             
2 26 Feb 07 jari 632             if (stop) {
2 26 Feb 07 jari 633                 throw new AbortException();
2 26 Feb 07 jari 634             }
2 26 Feb 07 jari 635             //event.setIntValue(i);
2 26 Feb 07 jari 636             //event.setDescription("Imputing missing values: Current gene = " + (i+ 1));
2 26 Feb 07 jari 637             //fireValueChanged(event);
2 26 Feb 07 jari 638             
2 26 Feb 07 jari 639             float[] currentRow = new float[numCols];
2 26 Feb 07 jari 640             float[] currentOrigRow = new float[numCols];
2 26 Feb 07 jari 641             for (int j = 0; j < numCols; j++) {
2 26 Feb 07 jari 642                 currentRow[j] = inputMatrix.A[i][j];
2 26 Feb 07 jari 643                 currentOrigRow[j] = inputMatrix.A[i][j];
2 26 Feb 07 jari 644             }
2 26 Feb 07 jari 645             for (int k = 0; k < numCols; k++) {
2 26 Feb 07 jari 646                 if (Float.isNaN(inputMatrix.A[i][k])) {
2 26 Feb 07 jari 647                     currentRow[k] = getMean(currentOrigRow);
2 26 Feb 07 jari 648                 }
2 26 Feb 07 jari 649             }
2 26 Feb 07 jari 650             
2 26 Feb 07 jari 651             for (int l = 0; l < numCols; l++) {
2 26 Feb 07 jari 652                 resultMatrix.A[i][l] = currentRow[l];
2 26 Feb 07 jari 653             }
2 26 Feb 07 jari 654         }
2 26 Feb 07 jari 655         
2 26 Feb 07 jari 656         return resultMatrix;
2 26 Feb 07 jari 657     }    
2 26 Feb 07 jari 658     
2 26 Feb 07 jari 659     private boolean isMissingValues(FloatMatrix mat, int row) {//returns true if the row of the matrix has any NaN values
2 26 Feb 07 jari 660         
2 26 Feb 07 jari 661         for (int i = 0; i < mat.getColumnDimension(); i++) {
2 26 Feb 07 jari 662             if (Float.isNaN(mat.A[row][i])) {
2 26 Feb 07 jari 663                 return true;
2 26 Feb 07 jari 664             }
2 26 Feb 07 jari 665         }
2 26 Feb 07 jari 666         
2 26 Feb 07 jari 667         return false;
2 26 Feb 07 jari 668     }   
2 26 Feb 07 jari 669     
2 26 Feb 07 jari 670     private Vector getValidGenes(int gene, FloatMatrix mat, Vector validExpts) { //returns the indices of those genes in "mat" that have valid values for all the validExpts
2 26 Feb 07 jari 671         Vector validGenes = new Vector();
2 26 Feb 07 jari 672         
2 26 Feb 07 jari 673         for (int i = 0; i < mat.getRowDimension(); i++) {
2 26 Feb 07 jari 674             if ((hasAllExpts(i, mat, validExpts)) && (gene != i)){//returns true if gene i in "mat" has valid values for all the validExpts
2 26 Feb 07 jari 675                 validGenes.add(new Integer(i));
2 26 Feb 07 jari 676             }
2 26 Feb 07 jari 677         }
2 26 Feb 07 jari 678         
2 26 Feb 07 jari 679         if (validGenes.size() < numNeighbors) { // if the number of valid genes is < k, other genes will be added to validGenes in increasing order of Euclidean distance until validGenes.size() = k
2 26 Feb 07 jari 680             int additionalGenesNeeded = numNeighbors - validGenes.size();
2 26 Feb 07 jari 681             Vector additionalGenes = getAdditionalGenes(gene, additionalGenesNeeded, validGenes, mat);
2 26 Feb 07 jari 682             for (int i = 0; i < additionalGenes.size(); i++) {
2 26 Feb 07 jari 683                 validGenes.add(additionalGenes.get(i));
2 26 Feb 07 jari 684             }
2 26 Feb 07 jari 685         }
2 26 Feb 07 jari 686         
2 26 Feb 07 jari 687         return validGenes;
2 26 Feb 07 jari 688     } 
2 26 Feb 07 jari 689     
2 26 Feb 07 jari 690     private Vector getAdditionalGenes(int currentGene, int numGenesNeeded, Vector alreadyPresentGenes, FloatMatrix mat) {
2 26 Feb 07 jari 691         Vector additionalGenes = new Vector();
2 26 Feb 07 jari 692         Vector allGenes = new Vector();
2 26 Feb 07 jari 693         Vector geneDistances = new Vector();
2 26 Feb 07 jari 694         
2 26 Feb 07 jari 695         for (int i = 0; i < mat.getRowDimension(); i++) {
2 26 Feb 07 jari 696             if (i != currentGene) {
2 26 Feb 07 jari 697                 float currentDistance = ExperimentUtil.geneEuclidianDistance(mat, null, i, currentGene, factor);
2 26 Feb 07 jari 698                 geneDistances.add(new Float(currentDistance));
2 26 Feb 07 jari 699                 allGenes.add(new Integer(i));
2 26 Feb 07 jari 700             }
2 26 Feb 07 jari 701         }
2 26 Feb 07 jari 702         
2 26 Feb 07 jari 703         float[] geneDistancesArray = new float[geneDistances.size()];
2 26 Feb 07 jari 704         for (int i = 0; i < geneDistances.size(); i++) {
2 26 Feb 07 jari 705             float currentDist = ((Float)geneDistances.get(i)).floatValue();
2 26 Feb 07 jari 706             geneDistancesArray[i] = currentDist;
2 26 Feb 07 jari 707         }
2 26 Feb 07 jari 708         
2 26 Feb 07 jari 709         QSort sortGeneDistances = new QSort(geneDistancesArray);
2 26 Feb 07 jari 710         float[] sortedDistances = sortGeneDistances.getSorted();
2 26 Feb 07 jari 711         int[] sortedDistanceIndices = sortGeneDistances.getOrigIndx();
2 26 Feb 07 jari 712         
2 26 Feb 07 jari 713         int counter = 0;
2 26 Feb 07 jari 714         
2 26 Feb 07 jari 715         for (int i = 0; i < sortedDistanceIndices.length; i++) {
2 26 Feb 07 jari 716             int currentIndex = sortedDistanceIndices[i];
2 26 Feb 07 jari 717             int currentNearestGene = ((Integer)allGenes.get(currentIndex)).intValue();
2 26 Feb 07 jari 718             if (belongsIn(alreadyPresentGenes, currentNearestGene)) {
2 26 Feb 07 jari 719                 continue;
2 26 Feb 07 jari 720             } else {
2 26 Feb 07 jari 721                 additionalGenes.add(new Integer(currentNearestGene));
2 26 Feb 07 jari 722                 counter++;
2 26 Feb 07 jari 723                 if (counter >= numGenesNeeded) {
2 26 Feb 07 jari 724                     break;
2 26 Feb 07 jari 725                 }
2 26 Feb 07 jari 726             }
2 26 Feb 07 jari 727         }
2 26 Feb 07 jari 728         
2 26 Feb 07 jari 729         return additionalGenes;
2 26 Feb 07 jari 730     }   
2 26 Feb 07 jari 731     
2 26 Feb 07 jari 732     Vector getKNearestGenes(int gene, int k, FloatMatrix mat, Vector geneSubset, Vector nonMissingExpts) {
2 26 Feb 07 jari 733         Vector allValidGenes = new Vector();
2 26 Feb 07 jari 734         Vector nearestGenes = new Vector();
2 26 Feb 07 jari 735         Vector geneDistances = new Vector();
2 26 Feb 07 jari 736         for (int i = 0; i < geneSubset.size(); i++) {
2 26 Feb 07 jari 737             int currentGene = ((Integer)geneSubset.get(i)).intValue();
2 26 Feb 07 jari 738             if (gene != currentGene) {
2 26 Feb 07 jari 739                 float currentDistance = ExperimentUtil.geneEuclidianDistance(mat, null, gene, currentGene, factor);
2 26 Feb 07 jari 740                 //System.out.println("Current distance = " + currentDistance);
2 26 Feb 07 jari 741                 geneDistances.add(new Float(currentDistance));
2 26 Feb 07 jari 742                 allValidGenes.add(new Integer(currentGene));
2 26 Feb 07 jari 743             }
2 26 Feb 07 jari 744         }
2 26 Feb 07 jari 745         
2 26 Feb 07 jari 746         float[] geneDistancesArray = new float[geneDistances.size()];
2 26 Feb 07 jari 747         for (int i = 0; i < geneDistances.size(); i++) {
2 26 Feb 07 jari 748             float currentDist = ((Float)geneDistances.get(i)).floatValue();
2 26 Feb 07 jari 749             geneDistancesArray[i] = currentDist;
2 26 Feb 07 jari 750         }
2 26 Feb 07 jari 751         
2 26 Feb 07 jari 752         QSort sortGeneDistances = new QSort(geneDistancesArray);
2 26 Feb 07 jari 753         float[] sortedDistances = sortGeneDistances.getSorted();
2 26 Feb 07 jari 754         int[] sortedDistanceIndices = sortGeneDistances.getOrigIndx();
2 26 Feb 07 jari 755         
2 26 Feb 07 jari 756         for (int i = 0; i < k; i++) {
2 26 Feb 07 jari 757             int currentGeneIndex = sortedDistanceIndices[i];
2 26 Feb 07 jari 758             int currentNearestGene = ((Integer)allValidGenes.get(currentGeneIndex)).intValue();
2 26 Feb 07 jari 759             nearestGenes.add(new Integer(currentNearestGene));
2 26 Feb 07 jari 760         }
2 26 Feb 07 jari 761         
2 26 Feb 07 jari 762         return nearestGenes;
2 26 Feb 07 jari 763     } 
2 26 Feb 07 jari 764     
2 26 Feb 07 jari 765     private float getExptWeightedMean(int gene, int expt, Vector geneVector, FloatMatrix mat) {
2 26 Feb 07 jari 766         float weightedMean = 0.0f;
2 26 Feb 07 jari 767         int validN = 0;
2 26 Feb 07 jari 768         float numerator = 0.0f;
2 26 Feb 07 jari 769         float recipNeighborDistances[] = new float[geneVector.size()];
2 26 Feb 07 jari 770         for (int i = 0; i < recipNeighborDistances.length; i++) {
2 26 Feb 07 jari 771             int currentGene = ((Integer)geneVector.get(i)).intValue();
2 26 Feb 07 jari 772             if (!Float.isNaN(mat.A[currentGene][expt])) {
2 26 Feb 07 jari 773                 float distance = ExperimentUtil.geneEuclidianDistance(mat, null, gene, currentGene, factor);
2 26 Feb 07 jari 774                 if (distance == 0.0f) {
2 26 Feb 07 jari 775                     distance = Float.MIN_VALUE;
2 26 Feb 07 jari 776                 }
2 26 Feb 07 jari 777                 recipNeighborDistances[i] = (float)(1.0f/distance);
2 26 Feb 07 jari 778                 numerator = numerator + (float)(recipNeighborDistances[i]*mat.A[currentGene][expt]);
2 26 Feb 07 jari 779                 validN++;
2 26 Feb 07 jari 780             } else {
2 26 Feb 07 jari 781                 recipNeighborDistances[i] = 0.0f;
2 26 Feb 07 jari 782             }
2 26 Feb 07 jari 783             
2 26 Feb 07 jari 784         }
2 26 Feb 07 jari 785         
2 26 Feb 07 jari 786         float denominator = 0.0f;
2 26 Feb 07 jari 787         for (int i = 0; i < recipNeighborDistances.length; i++) {
2 26 Feb 07 jari 788             denominator = denominator + recipNeighborDistances[i];
2 26 Feb 07 jari 789         }
2 26 Feb 07 jari 790         
2 26 Feb 07 jari 791         weightedMean = (float)(numerator/(float)denominator);
2 26 Feb 07 jari 792         return weightedMean;
2 26 Feb 07 jari 793     }   
2 26 Feb 07 jari 794     
2 26 Feb 07 jari 795     private float getMean(float[] row) {
2 26 Feb 07 jari 796         float mean = 0.0f;
2 26 Feb 07 jari 797         int validN = 0;
2 26 Feb 07 jari 798         
2 26 Feb 07 jari 799         for (int i = 0; i < row.length; i++) {
2 26 Feb 07 jari 800             if (!Float.isNaN(row[i])) {
2 26 Feb 07 jari 801                 mean = mean + row[i];
2 26 Feb 07 jari 802                 validN++;
2 26 Feb 07 jari 803             }
2 26 Feb 07 jari 804         }
2 26 Feb 07 jari 805         
2 26 Feb 07 jari 806         if (validN == 0) {
2 26 Feb 07 jari 807             validN = 1; // if the whole row is NaN, it will be set to zero;
2 26 Feb 07 jari 808         }
2 26 Feb 07 jari 809         
2 26 Feb 07 jari 810         mean = (float)(mean / validN);
2 26 Feb 07 jari 811         
2 26 Feb 07 jari 812         return mean;
2 26 Feb 07 jari 813     } 
2 26 Feb 07 jari 814     
2 26 Feb 07 jari 815     private boolean hasAllExpts(int gene, FloatMatrix mat, Vector validExpts) {//returns true if "gene" in "mat" has valid values for all the validExpts
2 26 Feb 07 jari 816         
2 26 Feb 07 jari 817         for (int i = 0; i < validExpts.size(); i++) {
2 26 Feb 07 jari 818             int expIndex = ((Integer)validExpts.get(i)).intValue();
2 26 Feb 07 jari 819             if (Float.isNaN(mat.A[gene][expIndex])) {
2 26 Feb 07 jari 820                 return false;
2 26 Feb 07 jari 821             }
2 26 Feb 07 jari 822         }
2 26 Feb 07 jari 823         
2 26 Feb 07 jari 824         return true;
2 26 Feb 07 jari 825     }
2 26 Feb 07 jari 826     
2 26 Feb 07 jari 827     private boolean belongsIn(Vector geneVector, int gene) {
2 26 Feb 07 jari 828         for (int i = 0; i < geneVector.size(); i++) {
2 26 Feb 07 jari 829             int currentGene = ((Integer)geneVector.get(i)).intValue();
2 26 Feb 07 jari 830             if (gene == currentGene) {
2 26 Feb 07 jari 831                 return true;
2 26 Feb 07 jari 832             }
2 26 Feb 07 jari 833         }
2 26 Feb 07 jari 834         
2 26 Feb 07 jari 835         return false;
2 26 Feb 07 jari 836     }    
2 26 Feb 07 jari 837     
2 26 Feb 07 jari 838 }

Generated on Thu Jul 20 09:25:34 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)