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

Code
Comments
Other
Rev Date Author Line
2 26 Feb 07 jari 1 /*
2 26 Feb 07 jari 2  * COA.java
2 26 Feb 07 jari 3  *
2 26 Feb 07 jari 4  * Created on September 15, 2004, 3:30 PM
2 26 Feb 07 jari 5  */
2 26 Feb 07 jari 6
2 26 Feb 07 jari 7 package org.tigr.microarray.mev.cluster.algorithm.impl;
2 26 Feb 07 jari 8
2 26 Feb 07 jari 9 import java.util.Vector;
2 26 Feb 07 jari 10
2 26 Feb 07 jari 11 import org.tigr.microarray.mev.cluster.algorithm.AbortException;
2 26 Feb 07 jari 12 import org.tigr.microarray.mev.cluster.algorithm.AbstractAlgorithm;
2 26 Feb 07 jari 13 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmData;
2 26 Feb 07 jari 14 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmEvent;
2 26 Feb 07 jari 15 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmException;
2 26 Feb 07 jari 16 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmParameters;
2 26 Feb 07 jari 17 import org.tigr.util.FloatMatrix;
2 26 Feb 07 jari 18 import org.tigr.util.QSort;
2 26 Feb 07 jari 19
2 26 Feb 07 jari 20 import Jama.Matrix;
2 26 Feb 07 jari 21 import Jama.SingularValueDecomposition;
2 26 Feb 07 jari 22
2 26 Feb 07 jari 23 /**
2 26 Feb 07 jari 24  *
2 26 Feb 07 jari 25  * @author  nbhagaba
2 26 Feb 07 jari 26  */
2 26 Feb 07 jari 27 public class COA extends AbstractAlgorithm {
2 26 Feb 07 jari 28     
2 26 Feb 07 jari 29     private boolean stop = false;    
2 26 Feb 07 jari 30     
2 26 Feb 07 jari 31     public FloatMatrix Org;     // Backup of original input matrix.
2 26 Feb 07 jari 32     private int numNeighbors;
2 26 Feb 07 jari 33     private float factor = 1.0f;
2 26 Feb 07 jari 34     public Matrix N;            // Input matrix as double values.
2 26 Feb 07 jari 35     public double[] R;          // Row sums.
2 26 Feb 07 jari 36     public double[] C;          // Column sums.
2 26 Feb 07 jari 37     public double N_Sum = 0;    // Sum of all values in N.
2 26 Feb 07 jari 38     public double[][] P;        // Matrix of N[i,j]/N_Sum values.
2 26 Feb 07 jari 39     public Matrix X;            // Matrix of (P[i,j]-(Ri*Cj))/sqrt(Ri*Cj) values.
2 26 Feb 07 jari 40     public Matrix Xt;           // Transpose of X.
2 26 Feb 07 jari 41     public Matrix DcRoot;           // Diagonal matrix of C.
2 26 Feb 07 jari 42     public Matrix Dr;           // Diagonal matrix of R.
2 26 Feb 07 jari 43     public Matrix B;            // Intermediate matrix of the algorithm.
2 26 Feb 07 jari 44     public Matrix geneUMatrix;  // Coordinates of genes.
2 26 Feb 07 jari 45     public Matrix exptUMatrix;  // Coordinates of experiments.
2 26 Feb 07 jari 46     public FloatMatrix gene;    // Coordinates of genes as FloatMatrix.
2 26 Feb 07 jari 47     public FloatMatrix expt;    // Coordinates of experiments as FloatMatrix.
2 26 Feb 07 jari 48     public double[] G_Sums;
2 26 Feb 07 jari 49     public Matrix G;
2 26 Feb 07 jari 50     
2 26 Feb 07 jari 51     private int numGenes, numExps;    
2 26 Feb 07 jari 52         
2 26 Feb 07 jari 53     
2 26 Feb 07 jari 54     /** Creates a new instance of COA */
2 26 Feb 07 jari 55     public COA() {
2 26 Feb 07 jari 56     }
2 26 Feb 07 jari 57     
2 26 Feb 07 jari 58     /** This method should interrupt the calculation.
2 26 Feb 07 jari 59      */
2 26 Feb 07 jari 60     public void abort() {
2 26 Feb 07 jari 61         stop = true;        
2 26 Feb 07 jari 62     }
2 26 Feb 07 jari 63     
2 26 Feb 07 jari 64     /** This method execute calculation and return result,
2 26 Feb 07 jari 65      * stored in <code>AlgorithmData</code> class.
2 26 Feb 07 jari 66      *
2 26 Feb 07 jari 67      * @param data the data to be calculated.
2 26 Feb 07 jari 68      */
2 26 Feb 07 jari 69     public AlgorithmData execute(AlgorithmData data) throws AlgorithmException {
2 26 Feb 07 jari 70   AlgorithmParameters map = data.getParams();
2 26 Feb 07 jari 71   numNeighbors = map.getInt("numNeighbors", 10);        
2 26 Feb 07 jari 72         
2 26 Feb 07 jari 73         Org = data.getMatrix("experiment");
2 26 Feb 07 jari 74   numGenes = this.Org.getRowDimension();
2 26 Feb 07 jari 75   numExps = this.Org.getColumnDimension(); 
2 26 Feb 07 jari 76         AlgorithmData result = Analysis(Org);
2 26 Feb 07 jari 77         return result;        
2 26 Feb 07 jari 78     }
2 26 Feb 07 jari 79     
2 26 Feb 07 jari 80     /* Turns A into a matrix with A values on the diagonal and 0 everywhere else. */
2 26 Feb 07 jari 81     public Matrix diagonal (double[] A) {
2 26 Feb 07 jari 82         System.out.println("Entered diagonal method");
2 26 Feb 07 jari 83         Matrix D = new Matrix(A.length,A.length);
2 26 Feb 07 jari 84         System.out.println("Created diagonal matrix");
2 26 Feb 07 jari 85         for (int i = 0; i < D.getRowDimension(); i++) {
2 26 Feb 07 jari 86             for (int j = 0; j < D.getColumnDimension(); j++) {
2 26 Feb 07 jari 87                 D.set(i, j, 0d);
2 26 Feb 07 jari 88             }
2 26 Feb 07 jari 89         }        
2 26 Feb 07 jari 90         for (int i=0; i<D.getRowDimension(); i++) {
2 26 Feb 07 jari 91             D.set(i,i,Math.sqrt(A[i]));
2 26 Feb 07 jari 92         }
2 26 Feb 07 jari 93         
2 26 Feb 07 jari 94         return D;
2 26 Feb 07 jari 95     }    
2 26 Feb 07 jari 96     
2 26 Feb 07 jari 97     public Matrix diagonalRoot(double[] A) {
2 26 Feb 07 jari 98         //System.out.println("Entered Diagonal Root method");
2 26 Feb 07 jari 99         Matrix D = new Matrix(A.length,A.length);
2 26 Feb 07 jari 100         //System.out.println("Created Diagonal Root Matrix");
2 26 Feb 07 jari 101         
2 26 Feb 07 jari 102         for (int i = 0; i < D.getRowDimension(); i++) {
2 26 Feb 07 jari 103             for (int j = 0; j < D.getColumnDimension(); j++) {
2 26 Feb 07 jari 104                 D.set(i, j, 0d);
2 26 Feb 07 jari 105             }
2 26 Feb 07 jari 106         }          
2 26 Feb 07 jari 107         
2 26 Feb 07 jari 108         for (int i=0; i<D.getRowDimension(); i++) {
2 26 Feb 07 jari 109             D.set(i,i,A[i]);
2 26 Feb 07 jari 110         }
2 26 Feb 07 jari 111         
2 26 Feb 07 jari 112         return D;        
2 26 Feb 07 jari 113     }
2 26 Feb 07 jari 114     
2 26 Feb 07 jari 115     /* Determines the average of all values in p */
2 26 Feb 07 jari 116     public double mean(double[] p) {
2 26 Feb 07 jari 117         double sum = 0;  // sum of all the elements
2 26 Feb 07 jari 118         
2 26 Feb 07 jari 119         for (int i=0; i<p.length; i++) {
2 26 Feb 07 jari 120             sum += p[i];
2 26 Feb 07 jari 121         }
2 26 Feb 07 jari 122         
2 26 Feb 07 jari 123         return (sum/p.length);
2 26 Feb 07 jari 124     }
2 26 Feb 07 jari 125     
2 26 Feb 07 jari 126     private FloatMatrix imputeKNearestMatrix(FloatMatrix inputMatrix, int k) throws AlgorithmException {
2 26 Feb 07 jari 127         int numRows = inputMatrix.getRowDimension();
2 26 Feb 07 jari 128         int numCols = inputMatrix.getColumnDimension();
2 26 Feb 07 jari 129         FloatMatrix resultMatrix = new FloatMatrix(numRows, numCols);
2 26 Feb 07 jari 130         
2 26 Feb 07 jari 131         AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, numGenes);
2 26 Feb 07 jari 132         event.setDescription("Imputing missing values");
2 26 Feb 07 jari 133         fireValueChanged(event);
2 26 Feb 07 jari 134         event.setId(AlgorithmEvent.PROGRESS_VALUE);
2 26 Feb 07 jari 135         
2 26 Feb 07 jari 136         for (int i = 0; i < numRows; i++) {
2 26 Feb 07 jari 137             if (stop) {
2 26 Feb 07 jari 138                 throw new AbortException();
2 26 Feb 07 jari 139             }
2 26 Feb 07 jari 140             //event.setIntValue(i);
2 26 Feb 07 jari 141             //event.setDescription("Imputing missing values: Current gene = " + (i+ 1));
2 26 Feb 07 jari 142             //fireValueChanged(event);
2 26 Feb 07 jari 143             
2 26 Feb 07 jari 144             if (isMissingValues(inputMatrix, i)) {
2 26 Feb 07 jari 145                 //System.out.println("gene " + i + " is missing values");
2 26 Feb 07 jari 146                 Vector nonMissingExpts = new Vector();
2 26 Feb 07 jari 147                 for (int j = 0; j < numCols; j++) {
2 26 Feb 07 jari 148                     if (!Float.isNaN(inputMatrix.A[i][j])) {
2 26 Feb 07 jari 149                         nonMissingExpts.add(new Integer(j));
2 26 Feb 07 jari 150                     }
2 26 Feb 07 jari 151                 }
2 26 Feb 07 jari 152                 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 153                 
2 26 Feb 07 jari 154                 //System.out.println(" Valid geneSubset.size() = " + geneSubset.size());
2 26 Feb 07 jari 155                 
2 26 Feb 07 jari 156                 /*
2 26 Feb 07 jari 157                 for (int j = 0; j < geneSubset.size(); j++) {
2 26 Feb 07 jari 158                     System.out.println(((Integer)geneSubset.get(j)).intValue());
2 26 Feb 07 jari 159                 }
2 26 Feb 07 jari 160                  */
2 26 Feb 07 jari 161                 //System.out.println("imputing KNN: current gene = " + i);
2 26 Feb 07 jari 162                 Vector kNearestGenes = getKNearestGenes(i, k, inputMatrix, geneSubset, nonMissingExpts);
2 26 Feb 07 jari 163                 
2 26 Feb 07 jari 164                 /*
2 26 Feb 07 jari 165                 System.out.println("k nearest genes of gene " + i + " : ");
2 26 Feb 07 jari 166                  
2 26 Feb 07 jari 167                 for (int j = 0; j < kNearestGenes.size(); j++) {
2 26 Feb 07 jari 168                     System.out.println("" + ((Integer)kNearestGenes.get(j)).intValue());
2 26 Feb 07 jari 169                 }
2 26 Feb 07 jari 170                  */
2 26 Feb 07 jari 171                 
2 26 Feb 07 jari 172                 //TESTED UPTO HERE -- 12/18/2002***********
2 26 Feb 07 jari 173                 //
2 26 Feb 07 jari 174                 /*
2 26 Feb 07 jari 175                 System.out.print("Gene " + i + " :\t");
2 26 Feb 07 jari 176                 for (int j = 0; j < numCols; j++) {
2 26 Feb 07 jari 177                     System.out.print("" +inputMatrix.A[i][j]);
2 26 Feb 07 jari 178                     System.out.print("\t");
2 26 Feb 07 jari 179                 }
2 26 Feb 07 jari 180                 System.out.println();
2 26 Feb 07 jari 181                 System.out.println("Matrix of k Nearest Genes");
2 26 Feb 07 jari 182                 printSubMatrix(kNearestGenes, inputMatrix);
2 26 Feb 07 jari 183                  */    //
2 26 Feb 07 jari 184                 for (int j = 0; j < numCols; j++) {
2 26 Feb 07 jari 185                     if (!Float.isNaN(inputMatrix.A[i][j])) {
2 26 Feb 07 jari 186                         resultMatrix.A[i][j] = inputMatrix.A[i][j];
2 26 Feb 07 jari 187                     } else {
2 26 Feb 07 jari 188                         
2 26 Feb 07 jari 189                         //System.out.println("just before entering getExptMean(): kNearestGenes.size() = " + kNearestGenes.size());
2 26 Feb 07 jari 190                         
2 26 Feb 07 jari 191                         //resultMatrix.A[i][j] = getExptMean(j, kNearestGenes, inputMatrix);
2 26 Feb 07 jari 192                         resultMatrix.A[i][j] = getExptWeightedMean(i, j, kNearestGenes, inputMatrix);
2 26 Feb 07 jari 193                     }
2 26 Feb 07 jari 194                 }
2 26 Feb 07 jari 195                 //DONE UPTO HERE
2 26 Feb 07 jari 196             }
2 26 Feb 07 jari 197             
2 26 Feb 07 jari 198             else {
2 26 Feb 07 jari 199                 for (int j = 0; j < numCols; j++) {
2 26 Feb 07 jari 200                     resultMatrix.A[i][j] = inputMatrix.A[i][j];
2 26 Feb 07 jari 201                 }
2 26 Feb 07 jari 202             }
2 26 Feb 07 jari 203         }
2 26 Feb 07 jari 204         
2 26 Feb 07 jari 205         return imputeRowAverageMatrix(resultMatrix);
2 26 Feb 07 jari 206     }  
2 26 Feb 07 jari 207     
2 26 Feb 07 jari 208     private FloatMatrix imputeRowAverageMatrix(FloatMatrix inputMatrix) throws AlgorithmException {
2 26 Feb 07 jari 209         int numRows = inputMatrix.getRowDimension();
2 26 Feb 07 jari 210         int numCols = inputMatrix.getColumnDimension();
2 26 Feb 07 jari 211         FloatMatrix resultMatrix = new FloatMatrix(numRows, numCols);
2 26 Feb 07 jari 212         
2 26 Feb 07 jari 213         AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, numGenes);
2 26 Feb 07 jari 214         fireValueChanged(event);
2 26 Feb 07 jari 215         event.setId(AlgorithmEvent.PROGRESS_VALUE);
2 26 Feb 07 jari 216         
2 26 Feb 07 jari 217         for (int i = 0; i < numRows; i++) {
2 26 Feb 07 jari 218             
2 26 Feb 07 jari 219             if (stop) {
2 26 Feb 07 jari 220                 throw new AbortException();
2 26 Feb 07 jari 221             }
2 26 Feb 07 jari 222             //event.setIntValue(i);
2 26 Feb 07 jari 223             //event.setDescription("Imputing missing values: Current gene = " + (i+ 1));
2 26 Feb 07 jari 224             //fireValueChanged(event);
2 26 Feb 07 jari 225             
2 26 Feb 07 jari 226             float[] currentRow = new float[numCols];
2 26 Feb 07 jari 227             float[] currentOrigRow = new float[numCols];
2 26 Feb 07 jari 228             for (int j = 0; j < numCols; j++) {
2 26 Feb 07 jari 229                 currentRow[j] = inputMatrix.A[i][j];
2 26 Feb 07 jari 230                 currentOrigRow[j] = inputMatrix.A[i][j];
2 26 Feb 07 jari 231             }
2 26 Feb 07 jari 232             for (int k = 0; k < numCols; k++) {
2 26 Feb 07 jari 233                 if (Float.isNaN(inputMatrix.A[i][k])) {
2 26 Feb 07 jari 234                     currentRow[k] = getMean(currentOrigRow);
2 26 Feb 07 jari 235                 }
2 26 Feb 07 jari 236             }
2 26 Feb 07 jari 237             
2 26 Feb 07 jari 238             for (int l = 0; l < numCols; l++) {
2 26 Feb 07 jari 239                 resultMatrix.A[i][l] = currentRow[l];
2 26 Feb 07 jari 240             }
2 26 Feb 07 jari 241         }
2 26 Feb 07 jari 242         
2 26 Feb 07 jari 243         return resultMatrix;
2 26 Feb 07 jari 244     }  
2 26 Feb 07 jari 245     
2 26 Feb 07 jari 246     private boolean isMissingValues(FloatMatrix mat, int row) {//returns true if the row of the matrix has any NaN values
2 26 Feb 07 jari 247         
2 26 Feb 07 jari 248         for (int i = 0; i < mat.getColumnDimension(); i++) {
2 26 Feb 07 jari 249             if (Float.isNaN(mat.A[row][i])) {
2 26 Feb 07 jari 250                 return true;
2 26 Feb 07 jari 251             }
2 26 Feb 07 jari 252         }
2 26 Feb 07 jari 253         
2 26 Feb 07 jari 254         return false;
2 26 Feb 07 jari 255     }   
2 26 Feb 07 jari 256     
2 26 Feb 07 jari 257     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 258         Vector validGenes = new Vector();
2 26 Feb 07 jari 259         
2 26 Feb 07 jari 260         for (int i = 0; i < mat.getRowDimension(); i++) {
2 26 Feb 07 jari 261             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 262                 validGenes.add(new Integer(i));
2 26 Feb 07 jari 263             }
2 26 Feb 07 jari 264         }
2 26 Feb 07 jari 265         
2 26 Feb 07 jari 266         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 267             int additionalGenesNeeded = numNeighbors - validGenes.size();
2 26 Feb 07 jari 268             Vector additionalGenes = getAdditionalGenes(gene, additionalGenesNeeded, validGenes, mat);
2 26 Feb 07 jari 269             for (int i = 0; i < additionalGenes.size(); i++) {
2 26 Feb 07 jari 270                 validGenes.add(additionalGenes.get(i));
2 26 Feb 07 jari 271             }
2 26 Feb 07 jari 272         }
2 26 Feb 07 jari 273         
2 26 Feb 07 jari 274         return validGenes;
2 26 Feb 07 jari 275     } 
2 26 Feb 07 jari 276     
2 26 Feb 07 jari 277     private Vector getAdditionalGenes(int currentGene, int numGenesNeeded, Vector alreadyPresentGenes, FloatMatrix mat) {
2 26 Feb 07 jari 278         Vector additionalGenes = new Vector();
2 26 Feb 07 jari 279         Vector allGenes = new Vector();
2 26 Feb 07 jari 280         Vector geneDistances = new Vector();
2 26 Feb 07 jari 281         
2 26 Feb 07 jari 282         for (int i = 0; i < mat.getRowDimension(); i++) {
2 26 Feb 07 jari 283             if (i != currentGene) {
2 26 Feb 07 jari 284                 float currentDistance = ExperimentUtil.geneEuclidianDistance(mat, null, i, currentGene, factor);
2 26 Feb 07 jari 285                 geneDistances.add(new Float(currentDistance));
2 26 Feb 07 jari 286                 allGenes.add(new Integer(i));
2 26 Feb 07 jari 287             }
2 26 Feb 07 jari 288         }
2 26 Feb 07 jari 289         
2 26 Feb 07 jari 290         float[] geneDistancesArray = new float[geneDistances.size()];
2 26 Feb 07 jari 291         for (int i = 0; i < geneDistances.size(); i++) {
2 26 Feb 07 jari 292             float currentDist = ((Float)geneDistances.get(i)).floatValue();
2 26 Feb 07 jari 293             geneDistancesArray[i] = currentDist;
2 26 Feb 07 jari 294         }
2 26 Feb 07 jari 295         
2 26 Feb 07 jari 296         QSort sortGeneDistances = new QSort(geneDistancesArray);
2 26 Feb 07 jari 297         float[] sortedDistances = sortGeneDistances.getSorted();
2 26 Feb 07 jari 298         int[] sortedDistanceIndices = sortGeneDistances.getOrigIndx();
2 26 Feb 07 jari 299         
2 26 Feb 07 jari 300         int counter = 0;
2 26 Feb 07 jari 301         
2 26 Feb 07 jari 302         for (int i = 0; i < sortedDistanceIndices.length; i++) {
2 26 Feb 07 jari 303             int currentIndex = sortedDistanceIndices[i];
2 26 Feb 07 jari 304             int currentNearestGene = ((Integer)allGenes.get(currentIndex)).intValue();
2 26 Feb 07 jari 305             if (belongsIn(alreadyPresentGenes, currentNearestGene)) {
2 26 Feb 07 jari 306                 continue;
2 26 Feb 07 jari 307             } else {
2 26 Feb 07 jari 308                 additionalGenes.add(new Integer(currentNearestGene));
2 26 Feb 07 jari 309                 counter++;
2 26 Feb 07 jari 310                 if (counter >= numGenesNeeded) {
2 26 Feb 07 jari 311                     break;
2 26 Feb 07 jari 312                 }
2 26 Feb 07 jari 313             }
2 26 Feb 07 jari 314         }
2 26 Feb 07 jari 315         
2 26 Feb 07 jari 316         return additionalGenes;
2 26 Feb 07 jari 317     }  
2 26 Feb 07 jari 318     
2 26 Feb 07 jari 319     // New.
2 26 Feb 07 jari 320     public void groups(int[][] selections, int numGroups) {       
2 26 Feb 07 jari 321         G = new Matrix(N.getRowDimension(),numGroups);
2 26 Feb 07 jari 322         G_Sums = new double[numGroups];
2 26 Feb 07 jari 323         
2 26 Feb 07 jari 324         for (int i=0; i<numGroups; i++) {
2 26 Feb 07 jari 325             // Find all values of G.
2 26 Feb 07 jari 326             double[][] sub = (N.getMatrix(0, N.getRowDimension()-1, selections[i])).getArrayCopy();
2 26 Feb 07 jari 327             for (int j=0; j<G.getRowDimension(); j++) {
2 26 Feb 07 jari 328                 G.set(j,i,mean(sub[j]));
2 26 Feb 07 jari 329             }
2 26 Feb 07 jari 330         
2 26 Feb 07 jari 331             // Find all values of G_Sums[i].
2 26 Feb 07 jari 332             for (int j=0; j<selections[i].length; j++) {
2 26 Feb 07 jari 333                 G_Sums[i] += C[selections[i][j]];
2 26 Feb 07 jari 334             }
2 26 Feb 07 jari 335         }
2 26 Feb 07 jari 336         
2 26 Feb 07 jari 337         System.out.println("G:");
2 26 Feb 07 jari 338         for (int i=0; i<G.getRowDimension(); i++) {
2 26 Feb 07 jari 339             for (int j=0; j<G.getColumnDimension(); j++) {
2 26 Feb 07 jari 340                 System.out.println("i: " + i + " , j: " + j + " , " + G.get(i,j));
2 26 Feb 07 jari 341             }
2 26 Feb 07 jari 342         }
2 26 Feb 07 jari 343         System.out.println();
2 26 Feb 07 jari 344         
2 26 Feb 07 jari 345         System.out.println("G_Sums:");
2 26 Feb 07 jari 346         for (int i=0; i<G_Sums.length; i++) {
2 26 Feb 07 jari 347             System.out.println("i: " + i + " , " + G_Sums[i]);
2 26 Feb 07 jari 348         }
2 26 Feb 07 jari 349         System.out.println();
2 26 Feb 07 jari 350     }
2 26 Feb 07 jari 351     
2 26 Feb 07 jari 352     
2 26 Feb 07 jari 353     public AlgorithmData Analysis(FloatMatrix input) throws AlgorithmException {
2 26 Feb 07 jari 354         
2 26 Feb 07 jari 355         input = imputeKNearestMatrix(input, numNeighbors); 
2 26 Feb 07 jari 356         
2 26 Feb 07 jari 357         // Prepare N.
2 26 Feb 07 jari 358         N = new Matrix(input.getRowDimension(),input.getColumnDimension());        
2 26 Feb 07 jari 359      
2 26 Feb 07 jari 360         boolean negFound = false;
2 26 Feb 07 jari 361         
2 26 Feb 07 jari 362         for (int i=0; i<input.getRowDimension(); i++) {
2 26 Feb 07 jari 363             if (negFound) break;
2 26 Feb 07 jari 364             for (int j=0; j<input.getColumnDimension(); j++) {
2 26 Feb 07 jari 365                 if (input.get(i,j) < 0f) {
2 26 Feb 07 jari 366                     negFound = true;
2 26 Feb 07 jari 367                     break;
2 26 Feb 07 jari 368                 }
2 26 Feb 07 jari 369             }
2 26 Feb 07 jari 370             //if (negFound) break;
2 26 Feb 07 jari 371         }
2 26 Feb 07 jari 372         
2 26 Feb 07 jari 373         if (negFound) {
2 26 Feb 07 jari 374             float min = Float.POSITIVE_INFINITY;
2 26 Feb 07 jari 375             for (int i=0; i<input.getRowDimension(); i++) {
2 26 Feb 07 jari 376                 for (int j=0; j<input.getColumnDimension(); j++) {
2 26 Feb 07 jari 377                     min = Math.min(min,input.get(i,j));
2 26 Feb 07 jari 378                 }
2 26 Feb 07 jari 379             }   
2 26 Feb 07 jari 380             
2 26 Feb 07 jari 381             //System.out.println("min = " + min);
2 26 Feb 07 jari 382             double additionFactor = (-1d)*min;
2 26 Feb 07 jari 383             //Make the matrix have all non-negative values.
2 26 Feb 07 jari 384             for (int i=0; i<input.getRowDimension(); i++) {
2 26 Feb 07 jari 385                 for (int j=0; j<input.getColumnDimension(); j++) {
2 26 Feb 07 jari 386                     N.set(i,j, (double)(input.get(i,j) + additionFactor));
2 26 Feb 07 jari 387                 }
2 26 Feb 07 jari 388             } 
2 26 Feb 07 jari 389             
2 26 Feb 07 jari 390         } else {
2 26 Feb 07 jari 391             for (int i=0; i<input.getRowDimension(); i++) {
2 26 Feb 07 jari 392                 for (int j=0; j<input.getColumnDimension(); j++) {
2 26 Feb 07 jari 393                     N.set(i,j,(double)input.get(i,j));
2 26 Feb 07 jari 394                 }
2 26 Feb 07 jari 395             }           
2 26 Feb 07 jari 396         }
2 26 Feb 07 jari 397         
2 26 Feb 07 jari 398         //System.out.println("finished making matrix positive");
2 26 Feb 07 jari 399         
2 26 Feb 07 jari 400         // Find the minimum value in the input matrix.
2 26 Feb 07 jari 401         /*
2 26 Feb 07 jari 402         float min = (float)(0);
2 26 Feb 07 jari 403         for (int i=0; i<input.getRowDimension(); i++) {
2 26 Feb 07 jari 404             for (int j=0; j<input.getColumnDimension(); j++) {
2 26 Feb 07 jari 405                 min = Math.min(min,input.get(i,j));
2 26 Feb 07 jari 406             }
2 26 Feb 07 jari 407         }
2 26 Feb 07 jari 408          */
2 26 Feb 07 jari 409         
2 26 Feb 07 jari 410         // Subtract the minimum value from all values in the input matrix
2 26 Feb 07 jari 411         // (Add, instead of subtract, if minimum value is negative) and place
2 26 Feb 07 jari 412         // the value, as a double, in N.
2 26 Feb 07 jari 413         /*
2 26 Feb 07 jari 414         for (int i=0; i<input.getRowDimension(); i++) {
2 26 Feb 07 jari 415             for (int j=0; j<input.getColumnDimension(); j++) {
2 26 Feb 07 jari 416                 N.set(i,j,(double)input.get(i,j)-min);
2 26 Feb 07 jari 417             }
2 26 Feb 07 jari 418         }
2 26 Feb 07 jari 419          */
2 26 Feb 07 jari 420         
2 26 Feb 07 jari 421       
2 26 Feb 07 jari 422         
2 26 Feb 07 jari 423         // If N has more columns than rows, transpose it.
2 26 Feb 07 jari 424         /*
2 26 Feb 07 jari 425         if (N.getColumnDimension() > N.getRowDimension()) {
2 26 Feb 07 jari 426             N = N.transpose();
2 26 Feb 07 jari 427         }
2 26 Feb 07 jari 428         */
2 26 Feb 07 jari 429         // Determine N_Sum
2 26 Feb 07 jari 430         for (int i=0; i<N.getRowDimension(); i++) {
2 26 Feb 07 jari 431             for (int j=0; j<N.getColumnDimension(); j++) {
2 26 Feb 07 jari 432                 N_Sum += N.get(i,j);
2 26 Feb 07 jari 433             }
2 26 Feb 07 jari 434         }
2 26 Feb 07 jari 435         //System.out.println("N_sum = " + N_Sum);
2 26 Feb 07 jari 436         // Find values of R
2 26 Feb 07 jari 437         R = new double[N.getRowDimension()];
2 26 Feb 07 jari 438         for (int i=0; i<N.getRowDimension(); i++) {
2 26 Feb 07 jari 439             double N_iSum = 0;
2 26 Feb 07 jari 440             
2 26 Feb 07 jari 441             for (int x=0; x<N.getColumnDimension(); x++) {
2 26 Feb 07 jari 442                 N_iSum += N.get(i,x);
2 26 Feb 07 jari 443             }
2 26 Feb 07 jari 444             
2 26 Feb 07 jari 445             R[i] = N_iSum/N_Sum;
2 26 Feb 07 jari 446         }
2 26 Feb 07 jari 447         
2 26 Feb 07 jari 448         // Find values of C
2 26 Feb 07 jari 449         C = new double[N.getColumnDimension()];
2 26 Feb 07 jari 450         for (int j=0; j<N.getColumnDimension(); j++) {
2 26 Feb 07 jari 451             double N_jSum = 0;
2 26 Feb 07 jari 452             
2 26 Feb 07 jari 453             for (int x=0; x<N.getRowDimension(); x++) {
2 26 Feb 07 jari 454                 N_jSum += N.get(x,j);
2 26 Feb 07 jari 455             }
2 26 Feb 07 jari 456             
2 26 Feb 07 jari 457             C[j] = N_jSum/N_Sum;
2 26 Feb 07 jari 458         }
2 26 Feb 07 jari 459         
2 26 Feb 07 jari 460         // Find values of P & X
2 26 Feb 07 jari 461         P = new double[N.getRowDimension()][N.getColumnDimension()];
2 26 Feb 07 jari 462         X = new Matrix(N.getRowDimension(),N.getColumnDimension());
2 26 Feb 07 jari 463         for (int i=0; i<N.getRowDimension(); i++) {
2 26 Feb 07 jari 464             for (int j=0; j<N.getColumnDimension(); j++) {
2 26 Feb 07 jari 465                 P[i][j] = N.get(i,j)/N_Sum;
2 26 Feb 07 jari 466                 double xValue = (P[i][j] - (R[i]*C[j])) / Math.sqrt(R[i]*C[j]);
2 26 Feb 07 jari 467                 X.set(i,j,xValue);
2 26 Feb 07 jari 468             }
2 26 Feb 07 jari 469         }
2 26 Feb 07 jari 470         
2 26 Feb 07 jari 471         SingularValueDecomposition svd = new SingularValueDecomposition(X);
2 26 Feb 07 jari 472         
2 26 Feb 07 jari 473         Matrix U = svd.getU();
2 26 Feb 07 jari 474         Matrix Lambda = svd.getS();
2 26 Feb 07 jari 475         Matrix V = svd.getV();
2 26 Feb 07 jari 476         
2 26 Feb 07 jari 477         geneUMatrix = U.times(Lambda);
2 26 Feb 07 jari 478         exptUMatrix = V.times(Lambda);
2 26 Feb 07 jari 479         
2 26 Feb 07 jari 480         for (int i = 0; i < geneUMatrix.getRowDimension(); i++) {
2 26 Feb 07 jari 481             for (int j = 0; j < geneUMatrix.getColumnDimension(); j++) {
2 26 Feb 07 jari 482                 double currVal = geneUMatrix.get(i,j);
2 26 Feb 07 jari 483                 geneUMatrix.set(i, j, currVal/Math.sqrt(R[i]));
2 26 Feb 07 jari 484             }
2 26 Feb 07 jari 485         }
2 26 Feb 07 jari 486         
2 26 Feb 07 jari 487         for (int i = 0; i < exptUMatrix.getRowDimension(); i++) {
2 26 Feb 07 jari 488             for (int j = 0; j < exptUMatrix.getColumnDimension(); j++) {
2 26 Feb 07 jari 489                 double currVal = exptUMatrix.get(i,j);
2 26 Feb 07 jari 490                 exptUMatrix.set(i, j, currVal/Math.sqrt(C[i]));
2 26 Feb 07 jari 491             }
2 26 Feb 07 jari 492         }
2 26 Feb 07 jari 493         
2 26 Feb 07 jari 494         //System.out.println("Finished computing R, C and X");
2 26 Feb 07 jari 495         
2 26 Feb 07 jari 496         // Determine Xt; Xt = transpose of X
2 26 Feb 07 jari 497         /*
2 26 Feb 07 jari 498         System.out.println("before Xt");
2 26 Feb 07 jari 499         
2 26 Feb 07 jari 500         Xt = X.transpose();
2 26 Feb 07 jari 501         
2 26 Feb 07 jari 502         System.out.println("Computed Xt");
2 26 Feb 07 jari 503         
2 26 Feb 07 jari 504         // Determine Dc; Dc = diagonal of C
2 26 Feb 07 jari 505         DcRoot = diagonalRoot(C);
2 26 Feb 07 jari 506         
2 26 Feb 07 jari 507         System.out.println("Computed DcRoot");
2 26 Feb 07 jari 508         // Determine Dr; Dr = diagonal of R
2 26 Feb 07 jari 509         Dr = diagonal(R);
2 26 Feb 07 jari 510         
2 26 Feb 07 jari 511         System.out.println("Computed Dr, before computing B");
2 26 Feb 07 jari 512         
2 26 Feb 07 jari 513         // Determine B; B = Dc*X*Dr*Xt*Dc
2 26 Feb 07 jari 514         B = (((DcRoot.times(Xt)).times(Dr)).times(X)).times(DcRoot);
2 26 Feb 07 jari 515         
2 26 Feb 07 jari 516         System.out.println("Finished computing B");
2 26 Feb 07 jari 517         
2 26 Feb 07 jari 518         // Get the Singular Value Decomposition of B
2 26 Feb 07 jari 519         SingularValueDecomposition mySVD = new SingularValueDecomposition(B);
2 26 Feb 07 jari 520         
2 26 Feb 07 jari 521         System.out.println("Finished SVD");
2 26 Feb 07 jari 522         
2 26 Feb 07 jari 523         Matrix T  = mySVD.getU();
2 26 Feb 07 jari 524         Matrix S = mySVD.getS();
2 26 Feb 07 jari 525         Matrix V = mySVD.getV();
2 26 Feb 07 jari 526         geneUMatrix = X.times(T);
2 26 Feb 07 jari 527         
2 26 Feb 07 jari 528         Matrix Q = mySVD.getU();
2 26 Feb 07 jari 529         Matrix D = mySVD.getS();
2 26 Feb 07 jari 530         //Matrix S= mySVD.getS();
2 26 Feb 07 jari 531         //Matrix V = mySVD.getV();
2 26 Feb 07 jari 532         for (int i=0;i<D.getRowDimension();i++) {
2 26 Feb 07 jari 533             D.set(i,i,1.0f/(float)Math.sqrt(D.get(i,i)));
2 26 Feb 07 jari 534         }
2 26 Feb 07 jari 535         T = X.times(Q.times(D));
2 26 Feb 07 jari 536         exptUMatrix = X.transpose().times(T);
2 26 Feb 07 jari 537         */
2 26 Feb 07 jari 538         
2 26 Feb 07 jari 539         // Fill gene will geneUMatrix's values as floats.
2 26 Feb 07 jari 540         gene = new FloatMatrix(geneUMatrix.getRowDimension(),geneUMatrix.getColumnDimension());
2 26 Feb 07 jari 541         for (int i=0; i<gene.getRowDimension(); i++) {
2 26 Feb 07 jari 542             for (int j=0; j<gene.getColumnDimension(); j++) {
2 26 Feb 07 jari 543                 gene.set(i,j,(float)geneUMatrix.get(i,j));
2 26 Feb 07 jari 544             }
2 26 Feb 07 jari 545         }
2 26 Feb 07 jari 546         
2 26 Feb 07 jari 547         // Fill expt will exptUMatrix's values as floats.
2 26 Feb 07 jari 548         expt = new FloatMatrix(exptUMatrix.getRowDimension(),exptUMatrix.getColumnDimension());
2 26 Feb 07 jari 549         for (int i=0; i<expt.getRowDimension(); i++) {
2 26 Feb 07 jari 550             for (int j=0; j<expt.getColumnDimension(); j++) {
2 26 Feb 07 jari 551                 expt.set(i,j,(float)exptUMatrix.get(i,j));
2 26 Feb 07 jari 552             }
2 26 Feb 07 jari 553         }
2 26 Feb 07 jari 554         
2 26 Feb 07 jari 555         FloatMatrix lambdaValues = new FloatMatrix(Lambda.getRowDimension(), 1); 
2 26 Feb 07 jari 556         
2 26 Feb 07 jari 557         for (int i = 0; i < Lambda.getRowDimension(); i++) {
2 26 Feb 07 jari 558             lambdaValues.set(i, 0, (float)(Lambda.get(i, i)));
2 26 Feb 07 jari 559         }
2 26 Feb 07 jari 560         
2 26 Feb 07 jari 561         AlgorithmData result = new AlgorithmData();
2 26 Feb 07 jari 562         result.addMatrix("gene",gene);
2 26 Feb 07 jari 563         result.addMatrix("expt",expt);
2 26 Feb 07 jari 564         result.addMatrix("lambdaValues", lambdaValues);
2 26 Feb 07 jari 565         
2 26 Feb 07 jari 566         //System.out.println("gene: rowDim = " + gene.getRowDimension() + ", colDim = " + gene.getColumnDimension());
2 26 Feb 07 jari 567         //System.out.println("expt: rowDim = " + expt.getRowDimension() + ", colDim = " + expt.getColumnDimension());
2 26 Feb 07 jari 568         
2 26 Feb 07 jari 569         return result;
2 26 Feb 07 jari 570     }   
2 26 Feb 07 jari 571     
2 26 Feb 07 jari 572     Vector getKNearestGenes(int gene, int k, FloatMatrix mat, Vector geneSubset, Vector nonMissingExpts) {
2 26 Feb 07 jari 573         Vector allValidGenes = new Vector();
2 26 Feb 07 jari 574         Vector nearestGenes = new Vector();
2 26 Feb 07 jari 575         Vector geneDistances = new Vector();
2 26 Feb 07 jari 576         for (int i = 0; i < geneSubset.size(); i++) {
2 26 Feb 07 jari 577             int currentGene = ((Integer)geneSubset.get(i)).intValue();
2 26 Feb 07 jari 578             if (gene != currentGene) {
2 26 Feb 07 jari 579                 float currentDistance = ExperimentUtil.geneEuclidianDistance(mat, null, gene, currentGene, factor);
2 26 Feb 07 jari 580                 //System.out.println("Current distance = " + currentDistance);
2 26 Feb 07 jari 581                 geneDistances.add(new Float(currentDistance));
2 26 Feb 07 jari 582                 allValidGenes.add(new Integer(currentGene));
2 26 Feb 07 jari 583             }
2 26 Feb 07 jari 584         }
2 26 Feb 07 jari 585         
2 26 Feb 07 jari 586         float[] geneDistancesArray = new float[geneDistances.size()];
2 26 Feb 07 jari 587         for (int i = 0; i < geneDistances.size(); i++) {
2 26 Feb 07 jari 588             float currentDist = ((Float)geneDistances.get(i)).floatValue();
2 26 Feb 07 jari 589             geneDistancesArray[i] = currentDist;
2 26 Feb 07 jari 590         }
2 26 Feb 07 jari 591         
2 26 Feb 07 jari 592         QSort sortGeneDistances = new QSort(geneDistancesArray);
2 26 Feb 07 jari 593         float[] sortedDistances = sortGeneDistances.getSorted();
2 26 Feb 07 jari 594         int[] sortedDistanceIndices = sortGeneDistances.getOrigIndx();
2 26 Feb 07 jari 595         
2 26 Feb 07 jari 596         for (int i = 0; i < k; i++) {
2 26 Feb 07 jari 597             int currentGeneIndex = sortedDistanceIndices[i];
2 26 Feb 07 jari 598             int currentNearestGene = ((Integer)allValidGenes.get(currentGeneIndex)).intValue();
2 26 Feb 07 jari 599             nearestGenes.add(new Integer(currentNearestGene));
2 26 Feb 07 jari 600         }
2 26 Feb 07 jari 601         
2 26 Feb 07 jari 602         return nearestGenes;
2 26 Feb 07 jari 603     } 
2 26 Feb 07 jari 604     
2 26 Feb 07 jari 605     private float getExptWeightedMean(int gene, int expt, Vector geneVector, FloatMatrix mat) {
2 26 Feb 07 jari 606         float weightedMean = 0.0f;
2 26 Feb 07 jari 607         int validN = 0;
2 26 Feb 07 jari 608         float numerator = 0.0f;
2 26 Feb 07 jari 609         float recipNeighborDistances[] = new float[geneVector.size()];
2 26 Feb 07 jari 610         for (int i = 0; i < recipNeighborDistances.length; i++) {
2 26 Feb 07 jari 611             int currentGene = ((Integer)geneVector.get(i)).intValue();
2 26 Feb 07 jari 612             if (!Float.isNaN(mat.A[currentGene][expt])) {
2 26 Feb 07 jari 613                 float distance = ExperimentUtil.geneEuclidianDistance(mat, null, gene, currentGene, factor);
2 26 Feb 07 jari 614                 if (distance == 0.0f) {
2 26 Feb 07 jari 615                     distance = Float.MIN_VALUE;
2 26 Feb 07 jari 616                 }
2 26 Feb 07 jari 617                 recipNeighborDistances[i] = (float)(1.0f/distance);
2 26 Feb 07 jari 618                 numerator = numerator + (float)(recipNeighborDistances[i]*mat.A[currentGene][expt]);
2 26 Feb 07 jari 619                 validN++;
2 26 Feb 07 jari 620             } else {
2 26 Feb 07 jari 621                 recipNeighborDistances[i] = 0.0f;
2 26 Feb 07 jari 622             }
2 26 Feb 07 jari 623             
2 26 Feb 07 jari 624         }
2 26 Feb 07 jari 625         
2 26 Feb 07 jari 626         float denominator = 0.0f;
2 26 Feb 07 jari 627         for (int i = 0; i < recipNeighborDistances.length; i++) {
2 26 Feb 07 jari 628             denominator = denominator + recipNeighborDistances[i];
2 26 Feb 07 jari 629         }
2 26 Feb 07 jari 630         
2 26 Feb 07 jari 631         weightedMean = (float)(numerator/(float)denominator);
2 26 Feb 07 jari 632         return weightedMean;
2 26 Feb 07 jari 633     }   
2 26 Feb 07 jari 634     
2 26 Feb 07 jari 635     private float getMean(float[] row) {
2 26 Feb 07 jari 636         float mean = 0.0f;
2 26 Feb 07 jari 637         int validN = 0;
2 26 Feb 07 jari 638         
2 26 Feb 07 jari 639         for (int i = 0; i < row.length; i++) {
2 26 Feb 07 jari 640             if (!Float.isNaN(row[i])) {
2 26 Feb 07 jari 641                 mean = mean + row[i];
2 26 Feb 07 jari 642                 validN++;
2 26 Feb 07 jari 643             }
2 26 Feb 07 jari 644         }
2 26 Feb 07 jari 645         
2 26 Feb 07 jari 646         if (validN == 0) {
2 26 Feb 07 jari 647             validN = 1; // if the whole row is NaN, it will be set to zero;
2 26 Feb 07 jari 648         }
2 26 Feb 07 jari 649         
2 26 Feb 07 jari 650         mean = (float)(mean / validN);
2 26 Feb 07 jari 651         
2 26 Feb 07 jari 652         return mean;
2 26 Feb 07 jari 653     } 
2 26 Feb 07 jari 654     
2 26 Feb 07 jari 655     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 656         
2 26 Feb 07 jari 657         for (int i = 0; i < validExpts.size(); i++) {
2 26 Feb 07 jari 658             int expIndex = ((Integer)validExpts.get(i)).intValue();
2 26 Feb 07 jari 659             if (Float.isNaN(mat.A[gene][expIndex])) {
2 26 Feb 07 jari 660                 return false;
2 26 Feb 07 jari 661             }
2 26 Feb 07 jari 662         }
2 26 Feb 07 jari 663         
2 26 Feb 07 jari 664         return true;
2 26 Feb 07 jari 665     }
2 26 Feb 07 jari 666     
2 26 Feb 07 jari 667     private boolean belongsIn(Vector geneVector, int gene) {
2 26 Feb 07 jari 668         for (int i = 0; i < geneVector.size(); i++) {
2 26 Feb 07 jari 669             int currentGene = ((Integer)geneVector.get(i)).intValue();
2 26 Feb 07 jari 670             if (gene == currentGene) {
2 26 Feb 07 jari 671                 return true;
2 26 Feb 07 jari 672             }
2 26 Feb 07 jari 673         }
2 26 Feb 07 jari 674         
2 26 Feb 07 jari 675         return false;
2 26 Feb 07 jari 676     }    
2 26 Feb 07 jari 677     
2 26 Feb 07 jari 678 }
2 26 Feb 07 jari 679

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