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

Code
Comments
Other
Rev Date Author Line
2 26 Feb 07 jari 1 /*
2 26 Feb 07 jari 2 Copyright @ 1999-2005, 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: DAM.java,v $
2 26 Feb 07 jari 7  * $Revision: 1.3 $
2 26 Feb 07 jari 8  * $Date: 2006/03/24 15:49:51 $
2 26 Feb 07 jari 9  * $Author: eleanorahowe $
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.text.DecimalFormat;
2 26 Feb 07 jari 15 import java.util.Vector;
2 26 Feb 07 jari 16
2 26 Feb 07 jari 17 import javax.swing.JOptionPane;
2 26 Feb 07 jari 18
2 26 Feb 07 jari 19 import org.tigr.microarray.mev.cluster.Cluster;
2 26 Feb 07 jari 20 import org.tigr.microarray.mev.cluster.Node;
2 26 Feb 07 jari 21 import org.tigr.microarray.mev.cluster.NodeList;
2 26 Feb 07 jari 22 import org.tigr.microarray.mev.cluster.algorithm.AbortException;
2 26 Feb 07 jari 23 import org.tigr.microarray.mev.cluster.algorithm.AbstractAlgorithm;
2 26 Feb 07 jari 24 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmData;
2 26 Feb 07 jari 25 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmEvent;
2 26 Feb 07 jari 26 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmException;
2 26 Feb 07 jari 27 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmParameters;
2 26 Feb 07 jari 28 import org.tigr.util.FloatMatrix;
2 26 Feb 07 jari 29
2 26 Feb 07 jari 30 import JSci.maths.statistics.TDistribution;
2 26 Feb 07 jari 31 import Jama.Matrix;
2 26 Feb 07 jari 32
2 26 Feb 07 jari 33 public class DAM extends AbstractAlgorithm {
2 26 Feb 07 jari 34     
2 26 Feb 07 jari 35     private boolean stop = false;
2 26 Feb 07 jari 36
2 26 Feb 07 jari 37     private Matrix expMatrix;
2 26 Feb 07 jari 38     private Matrix expMatrixTranspose;
2 26 Feb 07 jari 39     private Matrix responseMatrix;
2 26 Feb 07 jari 40     private Matrix testDataMatrix;
2 26 Feb 07 jari 41     private Matrix trainingMatrix;
2 26 Feb 07 jari 42
2 26 Feb 07 jari 43     private int numberOfGenes=0;
2 26 Feb 07 jari 44     private int numberOfSamples=0;
2 26 Feb 07 jari 45
2 26 Feb 07 jari 46     private int numberOfClasses=0;
2 26 Feb 07 jari 47     private int kValue=0;
2 26 Feb 07 jari 48     private int[] trainingIndices; // experiment indices array for classification 
2 26 Feb 07 jari 49     private int[] testIndices; // test data indices array for classification 
2 26 Feb 07 jari 50     private int[] classes;    // class number array for classification
2 26 Feb 07 jari 51
2 26 Feb 07 jari 52     private int whichAlgorithm = 0;
2 26 Feb 07 jari 53     private boolean isPDA = true;
2 26 Feb 07 jari 54
2 26 Feb 07 jari 55     private boolean preSelectGenes = true;
2 26 Feb 07 jari 56     private boolean performLOOCV = true;
2 26 Feb 07 jari 57
2 26 Feb 07 jari 58     private double alpha = 0.05; // default alpha is set to be 0.05
2 26 Feb 07 jari 59
2 26 Feb 07 jari 60     private int numberOfSelectedGenes=0;
2 26 Feb 07 jari 61
2 26 Feb 07 jari 62     private int[] geneRank;
2 26 Feb 07 jari 63     private int highestGeneRank;
2 26 Feb 07 jari 64     private Vector[] selectedGeneIndices;
2 26 Feb 07 jari 65
2 26 Feb 07 jari 66     private int[] usedGeneIndices;
2 26 Feb 07 jari 67     private int[] unusedGeneIndices;
2 26 Feb 07 jari 68
2 26 Feb 07 jari 69     final private int used=0;
2 26 Feb 07 jari 70     final private int unused=1;
2 26 Feb 07 jari 71
2 26 Feb 07 jari 72     private Vector[][] reducedGeneSetForA2;   // reduced gene set for A2 algorithm only
2 26 Feb 07 jari 73
2 26 Feb 07 jari 74     private Vector[] reducedGeneSet;
2 26 Feb 07 jari 75     private Vector[] clusters;
2 26 Feb 07 jari 76
2 26 Feb 07 jari 77     private Vector[] classified;
2 26 Feb 07 jari 78
2 26 Feb 07 jari 79     private Matrix classExpSumMatrix; // gene expression sum for each class
2 26 Feb 07 jari 80     private int[] classSampleArray; // number of samples in each class
2 26 Feb 07 jari 81  
2 26 Feb 07 jari 82     private Matrix geneComponentMatrix;
2 26 Feb 07 jari 83
2 26 Feb 07 jari 84     private Matrix [] beta;
2 26 Feb 07 jari 85     private Matrix [] A_Matrix;
2 26 Feb 07 jari 86     private Matrix [] C_Matrix;
2 26 Feb 07 jari 87
2 26 Feb 07 jari 88     private double [] cValues;
2 26 Feb 07 jari 89
2 26 Feb 07 jari 90     private boolean [] singularMatrix;
2 26 Feb 07 jari 91
2 26 Feb 07 jari 92     public AlgorithmData execute(AlgorithmData data) throws AlgorithmException {
2 26 Feb 07 jari 93
2 26 Feb 07 jari 94    expMatrix = getJamaMatrix(data.getMatrix("experiment"));
2 26 Feb 07 jari 95          expMatrixTranspose = expMatrix.transpose();
2 26 Feb 07 jari 96
2 26 Feb 07 jari 97    numberOfGenes = expMatrix.getRowDimension();
2 26 Feb 07 jari 98    numberOfSamples = expMatrix.getColumnDimension();
2 26 Feb 07 jari 99
2 26 Feb 07 jari 100    AlgorithmParameters map = data.getParams();
2 26 Feb 07 jari 101    int function = map.getInt("distance-function", COVARIANCE);
2 26 Feb 07 jari 102    float factor = map.getFloat("distance-factor", 1.0f);
2 26 Feb 07 jari 103    boolean absolute = map.getBoolean("distance-absolute", false);
2 26 Feb 07 jari 104    int mode = map.getInt("dam-mode", 0);
2 26 Feb 07 jari 105
2 26 Feb 07 jari 106          preSelectGenes = map.getBoolean("preSelectionGenes", false);
2 26 Feb 07 jari 107          whichAlgorithm = map.getInt("algorithmSelection", 0);
2 26 Feb 07 jari 108          isPDA = map.getBoolean("isPDA", true);
2 26 Feb 07 jari 109          preSelectGenes = map.getBoolean("preSelectGenes", true);
2 26 Feb 07 jari 110          numberOfClasses = map.getInt("numberOfClasses", 3);
2 26 Feb 07 jari 111          kValue = map.getInt("kValue", 3);
2 26 Feb 07 jari 112          alpha = (double) map.getFloat("alpha", 0.05f);
2 26 Feb 07 jari 113
2 26 Feb 07 jari 114  /*        System.out.println("");
2 26 Feb 07 jari 115          System.out.println("DAM.java: numberOfGenes = " + numberOfGenes);
2 26 Feb 07 jari 116          System.out.println("DAM.java: numberOfSamples = " + numberOfSamples);
2 26 Feb 07 jari 117    System.out.println("DAM.java: Algorithm = A" + whichAlgorithm);
2 26 Feb 07 jari 118    System.out.println("DAM.java: Perform PDA = " + isPDA);
2 26 Feb 07 jari 119    System.out.println("DAM.java: preSelectGenes = " + preSelectGenes);
2 26 Feb 07 jari 120    System.out.println("DAM.java: numberOfClasses = " + numberOfClasses);
2 26 Feb 07 jari 121    System.out.println("DAM.java: kValue = " + kValue);
2 26 Feb 07 jari 122    System.out.println("DAM.java: alpha = " + alpha);
2 26 Feb 07 jari 123          System.out.println("");
2 26 Feb 07 jari 124    */     
2 26 Feb 07 jari 125          singularMatrix = new boolean[numberOfSamples];
2 26 Feb 07 jari 126
2 26 Feb 07 jari 127          trainingIndices = data.getIntArray("trainingIndices");
2 26 Feb 07 jari 128          classes = data.getIntArray("classes");
2 26 Feb 07 jari 129          testIndices = data.getIntArray("testIndices");
2 26 Feb 07 jari 130
2 26 Feb 07 jari 131      //    System.out.println("DAM.java: trainingIndices: ");
2 26 Feb 07 jari 132 /*   for(int i=0; i<trainingIndices.length; i++) {
2 26 Feb 07 jari 133              System.out.print(trainingIndices[i] + ", ");
2 26 Feb 07 jari 134          }
2 26 Feb 07 jari 135          System.out.println("");
2 26 Feb 07 jari 136          System.out.println("");
2 26 Feb 07 jari 137
2 26 Feb 07 jari 138          System.out.println("DAM.java: classes: ");
2 26 Feb 07 jari 139    for(int i=0; i<classes.length; i++) {
2 26 Feb 07 jari 140              System.out.print(classes[i] + ", ");
2 26 Feb 07 jari 141          }
2 26 Feb 07 jari 142          System.out.println("");
2 26 Feb 07 jari 143          System.out.println("");
2 26 Feb 07 jari 144
2 26 Feb 07 jari 145          System.out.println("DAM.java: testIndices: ");
2 26 Feb 07 jari 146    for(int i=0; i<testIndices.length; i++) {
2 26 Feb 07 jari 147              System.out.print(testIndices[i] + ", ");
2 26 Feb 07 jari 148          }
2 26 Feb 07 jari 149
2 26 Feb 07 jari 150          System.out.println("");
2 26 Feb 07 jari 151          System.out.println("");
2 26 Feb 07 jari 152
2 26 Feb 07 jari 153
2 26 Feb 07 jari 154          System.out.print("DAM.java: expMatrix");
2 26 Feb 07 jari 155          for (int i=0; i<numberOfGenes; i++) {
2 26 Feb 07 jari 156              System.out.println("");
2 26 Feb 07 jari 157              for (int j=0; j<numberOfSamples; j++) {
2 26 Feb 07 jari 158                  System.out.print(expMatrix.get(i,j) + ",   ");
2 26 Feb 07 jari 159              }
2 26 Feb 07 jari 160              System.out.println("");
2 26 Feb 07 jari 161              System.out.println("");
2 26 Feb 07 jari 162          }
2 26 Feb 07 jari 163          System.out.println("");
2 26 Feb 07 jari 164 */
2 26 Feb 07 jari 165
2 26 Feb 07 jari 166
2 26 Feb 07 jari 167          responseMatrix = new Matrix(numberOfSamples, numberOfClasses);
2 26 Feb 07 jari 168
2 26 Feb 07 jari 169          for (int row=0; row<numberOfSamples; row++) {
2 26 Feb 07 jari 170              for (int column=0; column<numberOfClasses; column++) {
2 26 Feb 07 jari 171            responseMatrix.set(row, column, 0);
2 26 Feb 07 jari 172              }
2 26 Feb 07 jari 173          }
2 26 Feb 07 jari 174
2 26 Feb 07 jari 175    for (int i=0; i<trainingIndices.length; i++) {
2 26 Feb 07 jari 176        responseMatrix.set(trainingIndices[i], classes[i]-1, 1.0);
2 26 Feb 07 jari 177    }
2 26 Feb 07 jari 178
2 26 Feb 07 jari 179 /*
2 26 Feb 07 jari 180          System.out.println("");
2 26 Feb 07 jari 181          System.out.println("DAM.java: responseMatrix");
2 26 Feb 07 jari 182          for(int i=0; i<numberOfSamples; i++) {
2 26 Feb 07 jari 183              for(int j=0; j<numberOfClasses; j++) {
2 26 Feb 07 jari 184                  System.out.print(responseMatrix.get(i,j) + ",  ");
2 26 Feb 07 jari 185              }
2 26 Feb 07 jari 186              System.out.println("");
2 26 Feb 07 jari 187          }
2 26 Feb 07 jari 188          System.out.println("");
2 26 Feb 07 jari 189 */
2 26 Feb 07 jari 190
2 26 Feb 07 jari 191          if (trainingIndices.length > 0) {
2 26 Feb 07 jari 192              trainingMatrix = new Matrix(numberOfGenes, trainingIndices.length);
2 26 Feb 07 jari 193              trainingMatrix = expMatrix.getMatrix(0, numberOfGenes-1, trainingIndices);
2 26 Feb 07 jari 194
2 26 Feb 07 jari 195 /*
2 26 Feb 07 jari 196        System.out.println("");
2 26 Feb 07 jari 197        System.out.println("DAM.java: trainingMatrix "+ trainingMatrix.getRowDimension() + "x" + trainingMatrix.getColumnDimension());
2 26 Feb 07 jari 198
2 26 Feb 07 jari 199        for (int i=0; i<trainingMatrix.getRowDimension(); i++) {
2 26 Feb 07 jari 200      for (int j=0; j<trainingMatrix.getColumnDimension(); j++) {
2 26 Feb 07 jari 201          System.out.print(trainingMatrix.get(i,j) + ", ");
2 26 Feb 07 jari 202      }
2 26 Feb 07 jari 203      System.out.println("");
2 26 Feb 07 jari 204      System.out.println("");
2 26 Feb 07 jari 205        }
2 26 Feb 07 jari 206        System.out.println("");
2 26 Feb 07 jari 207 */
2 26 Feb 07 jari 208          }
2 26 Feb 07 jari 209
2 26 Feb 07 jari 210          if (testIndices.length > 0) {
2 26 Feb 07 jari 211              testDataMatrix = new Matrix(numberOfGenes, testIndices.length);
2 26 Feb 07 jari 212              testDataMatrix = expMatrix.getMatrix(0, numberOfGenes-1, testIndices);
2 26 Feb 07 jari 213
2 26 Feb 07 jari 214 /*
2 26 Feb 07 jari 215        System.out.println("");
2 26 Feb 07 jari 216        System.out.println("DAM.java: testMatrix " + testDataMatrix.getRowDimension() + "x" + testDataMatrix.getColumnDimension());
2 26 Feb 07 jari 217        System.out.println("");
2 26 Feb 07 jari 218 */
2 26 Feb 07 jari 219          }
2 26 Feb 07 jari 220
2 26 Feb 07 jari 221
2 26 Feb 07 jari 222          classExpSumMatrix = new Matrix(numberOfGenes,numberOfClasses);
2 26 Feb 07 jari 223          classExpSumMatrix = expMatrix.times(responseMatrix);
2 26 Feb 07 jari 224
2 26 Feb 07 jari 225 /*
2 26 Feb 07 jari 226          System.out.print("DAM.java: classExpSumMatrix ");
2 26 Feb 07 jari 227          for (int i=0; i<25; i++) {
2 26 Feb 07 jari 228              System.out.println("");
2 26 Feb 07 jari 229              for (int j=0; j<numberOfClasses; j++) {
2 26 Feb 07 jari 230                  System.out.print(classExpSumMatrix.get(i,j) + ",  ");
2 26 Feb 07 jari 231              }
2 26 Feb 07 jari 232          }
2 26 Feb 07 jari 233          System.out.println("");
2 26 Feb 07 jari 234          System.out.println("");
2 26 Feb 07 jari 235 */
2 26 Feb 07 jari 236
2 26 Feb 07 jari 237          classSampleArray = new int[numberOfClasses]; 
2 26 Feb 07 jari 238
2 26 Feb 07 jari 239 //         System.out.println("DAM.java: classSampleArray: ");
2 26 Feb 07 jari 240
2 26 Feb 07 jari 241          for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 242             for (int sample=0; sample<numberOfSamples; sample++) {
2 26 Feb 07 jari 243                 classSampleArray[classId] += responseMatrix.get(sample, classId); 
2 26 Feb 07 jari 244             }
2 26 Feb 07 jari 245             
2 26 Feb 07 jari 246  //           System.out.print(classSampleArray[classId] + ",  ");
2 26 Feb 07 jari 247          }
2 26 Feb 07 jari 248 //         System.out.println("");
2 26 Feb 07 jari 249 //         System.out.println("");
2 26 Feb 07 jari 250
2 26 Feb 07 jari 251          geneRank = new int[numberOfGenes];
2 26 Feb 07 jari 252
2 26 Feb 07 jari 253          reducedGeneSet = new Vector[2];
2 26 Feb 07 jari 254          reducedGeneSet[used] = new Vector();
2 26 Feb 07 jari 255          reducedGeneSet[unused] = new Vector();
2 26 Feb 07 jari 256
2 26 Feb 07 jari 257    reducedGeneSetForA2 = new Vector[numberOfSamples][2];
2 26 Feb 07 jari 258
2 26 Feb 07 jari 259          for (int leaveOutSample=0; leaveOutSample<numberOfSamples; leaveOutSample++) {
2 26 Feb 07 jari 260      reducedGeneSetForA2[leaveOutSample][used] = new Vector();
2 26 Feb 07 jari 261      reducedGeneSetForA2[leaveOutSample][unused] = new Vector();
2 26 Feb 07 jari 262          }
2 26 Feb 07 jari 263
2 26 Feb 07 jari 264          Matrix probFunction = new Matrix(numberOfSamples, numberOfClasses);
2 26 Feb 07 jari 265
2 26 Feb 07 jari 266          Matrix [] beta = new Matrix[numberOfClasses];
2 26 Feb 07 jari 267
2 26 Feb 07 jari 268    Vector[] classifiers = new Vector[numberOfClasses + 1];
2 26 Feb 07 jari 269    classified = new Vector[numberOfClasses + 1];
2 26 Feb 07 jari 270    Vector[] classifierPlusClassified = new Vector[numberOfClasses + 1];
2 26 Feb 07 jari 271
2 26 Feb 07 jari 272    for (int i = 0; i < numberOfClasses+1; i++) {
2 26 Feb 07 jari 273        classifiers[i] = new Vector();
2 26 Feb 07 jari 274        classified[i] = new Vector();
2 26 Feb 07 jari 275        classifierPlusClassified[i] = new Vector();
2 26 Feb 07 jari 276    }
2 26 Feb 07 jari 277
2 26 Feb 07 jari 278
2 26 Feb 07 jari 279          //jcb, add in initial classification to set the proper response matrix
2 26 Feb 07 jari 280          //to verify.  Note that responseMatrix is set to initial result.         
2 26 Feb 07 jari 281      //    probFunction = InitialClassification(expMatrix);         
2 26 Feb 07 jari 282      //    if(whichAlgorithm != 3) { //if the algorithm is not just a classification try 
2 26 Feb 07 jari 283                                     // to adjust response matrix to initial classification result
2 26 Feb 07 jari 284      //        responseMatrix = extractResult(probFunction);             
2 26 Feb 07 jari 285       //   }
2 26 Feb 07 jari 286          
2 26 Feb 07 jari 287          
2 26 Feb 07 jari 288          switch (whichAlgorithm) {
2 26 Feb 07 jari 289             case 0: 
2 26 Feb 07 jari 290                     probFunction = A0Algorithm(expMatrix);
2 26 Feb 07 jari 291                     break;
2 26 Feb 07 jari 292             case 1: 
2 26 Feb 07 jari 293                     probFunction = A1Algorithm(expMatrix);
2 26 Feb 07 jari 294                     break;
2 26 Feb 07 jari 295             case 2: 
2 26 Feb 07 jari 296                     probFunction = A2Algorithm(expMatrix);
2 26 Feb 07 jari 297                     break;
2 26 Feb 07 jari 298             case 3: 
2 26 Feb 07 jari 299                     probFunction = InitialClassification(expMatrix);
2 26 Feb 07 jari 300                     break;
2 26 Feb 07 jari 301             default: 
2 26 Feb 07 jari 302                     A0Algorithm(expMatrix);
2 26 Feb 07 jari 303                     break;
2 26 Feb 07 jari 304          }
2 26 Feb 07 jari 305
2 26 Feb 07 jari 306          if (probFunction == null) {
2 26 Feb 07 jari 307             throw new AbortException();
2 26 Feb 07 jari 308          }
2 26 Feb 07 jari 309
2 26 Feb 07 jari 310    for (int i = 0; i < trainingIndices.length; i++) {
2 26 Feb 07 jari 311        classifiers[classes[i]].add(new Integer(trainingIndices[i])); 
2 26 Feb 07 jari 312
2 26 Feb 07 jari 313 //             System.out.println("DAM.java - classifiers[" + classes[i] + "] add: " + trainingIndices[i]);
2 26 Feb 07 jari 314    }
2 26 Feb 07 jari 315
2 26 Feb 07 jari 316
2 26 Feb 07 jari 317          if (whichAlgorithm == 3) {  // for InitialClassification only
2 26 Feb 07 jari 318        for (int i = 0; i < classifiers.length; i++) {
2 26 Feb 07 jari 319      for (int j=0; j<classifiers[i].size(); j++) {
2 26 Feb 07 jari 320          classifierPlusClassified[i].add(classifiers[i].get(j));;
2 26 Feb 07 jari 321      }              
2 26 Feb 07 jari 322        }
2 26 Feb 07 jari 323        for (int i = 0; i < classified.length; i++) {
2 26 Feb 07 jari 324      for (int j=0; j<classified[i].size(); j++) {
2 26 Feb 07 jari 325          classifierPlusClassified[i].add(classified[i].get(j));;
2 26 Feb 07 jari 326      }
2 26 Feb 07 jari 327        }
2 26 Feb 07 jari 328          }
2 26 Feb 07 jari 329
2 26 Feb 07 jari 330   /*       System.out.println(" ");
2 26 Feb 07 jari 331    for (int i = 0; i < classifiers.length; i++) {
2 26 Feb 07 jari 332              System.out.println("DAM.java - classifiers[" + i + "].size() = " + classifiers[i].size());
2 26 Feb 07 jari 333    }
2 26 Feb 07 jari 334  
2 26 Feb 07 jari 335          System.out.println(" ");
2 26 Feb 07 jari 336    for (int i = 0; i < classified.length; i++) {
2 26 Feb 07 jari 337              System.out.println("DAM.java - classified[" + i + "].size() = " + classified[i].size());
2 26 Feb 07 jari 338    }
2 26 Feb 07 jari 339
2 26 Feb 07 jari 340          System.out.println(" ");
2 26 Feb 07 jari 341    for (int i = 0; i < classifierPlusClassified.length; i++) {
2 26 Feb 07 jari 342              System.out.println("DAM.java - classifierPlusClassified[" + i + "].size() = " + classifierPlusClassified[i].size());
2 26 Feb 07 jari 343    }
2 26 Feb 07 jari 344 */
2 26 Feb 07 jari 345
2 26 Feb 07 jari 346          clusters = new Vector[numberOfClasses*3]; 
2 26 Feb 07 jari 347          for (int i = 1; i <= numberOfClasses; i++) {
2 26 Feb 07 jari 348              clusters[i-1] = classifiers[i];
2 26 Feb 07 jari 349              clusters[i-1 + numberOfClasses] = classified[i];
2 26 Feb 07 jari 350              clusters[i-1 + numberOfClasses*2] = classifierPlusClassified[i];
2 26 Feb 07 jari 351          }
2 26 Feb 07 jari 352
2 26 Feb 07 jari 353  /*        for (int i = 1; i <= numberOfClasses; i++) {
2 26 Feb 07 jari 354              System.out.println("DAM.java - clusters 1 size " + clusters[i-1].size());
2 26 Feb 07 jari 355              System.out.println("DAM.java - clusters 2 size " + clusters[i-1 + numberOfClasses].size());
2 26 Feb 07 jari 356              System.out.println("DAM.java - clusters 3 size " + clusters[i-1 + numberOfClasses*2].size());
2 26 Feb 07 jari 357          }
2 26 Feb 07 jari 358  */
2 26 Feb 07 jari 359           
2 26 Feb 07 jari 360          Matrix means = getMeans(clusters);
2 26 Feb 07 jari 361          Matrix variances = getVariances(clusters, means);
2 26 Feb 07 jari 362
2 26 Feb 07 jari 363          Cluster result_cluster = new Cluster();
2 26 Feb 07 jari 364          NodeList nodeList = result_cluster.getNodeList();
2 26 Feb 07 jari 365          int[] features;
2 26 Feb 07 jari 366          for (int i=0; i<clusters.length; i++) {
2 26 Feb 07 jari 367              if (stop) {
2 26 Feb 07 jari 368                  throw new AbortException();
2 26 Feb 07 jari 369              }
2 26 Feb 07 jari 370              features = convert2int(clusters[i]);
2 26 Feb 07 jari 371
2 26 Feb 07 jari 372              Node node = new Node(features);
2 26 Feb 07 jari 373              nodeList.addNode(node);
2 26 Feb 07 jari 374          }
2 26 Feb 07 jari 375
2 26 Feb 07 jari 376          if (whichAlgorithm == 2) {
2 26 Feb 07 jari 377        usedGeneIndices = new int[reducedGeneSetForA2[0][used].size()];
2 26 Feb 07 jari 378        unusedGeneIndices = new int[reducedGeneSetForA2[0][unused].size()];
2 26 Feb 07 jari 379
2 26 Feb 07 jari 380        for (int i = 0; i < usedGeneIndices.length; i++) {
2 26 Feb 07 jari 381     usedGeneIndices[i] = ((Integer)(reducedGeneSetForA2[0][used].get(i))).intValue();
2 26 Feb 07 jari 382        }
2 26 Feb 07 jari 383
2 26 Feb 07 jari 384        for (int i = 0; i < unusedGeneIndices.length; i++) {
2 26 Feb 07 jari 385     unusedGeneIndices[i] = ((Integer)(reducedGeneSetForA2[0][unused].get(i))).intValue();
2 26 Feb 07 jari 386        }
2 26 Feb 07 jari 387          } else {
2 26 Feb 07 jari 388
2 26 Feb 07 jari 389        usedGeneIndices = new int[reducedGeneSet[used].size()];
2 26 Feb 07 jari 390        unusedGeneIndices = new int[reducedGeneSet[unused].size()];
2 26 Feb 07 jari 391
2 26 Feb 07 jari 392        for (int i = 0; i < usedGeneIndices.length; i++) {
2 26 Feb 07 jari 393     usedGeneIndices[i] = ((Integer)(reducedGeneSet[used].get(i))).intValue();
2 26 Feb 07 jari 394        }
2 26 Feb 07 jari 395
2 26 Feb 07 jari 396        for (int i = 0; i < unusedGeneIndices.length; i++) {
2 26 Feb 07 jari 397     unusedGeneIndices[i] = ((Integer)(reducedGeneSet[unused].get(i))).intValue();
2 26 Feb 07 jari 398        }
2 26 Feb 07 jari 399          }
2 26 Feb 07 jari 400
2 26 Feb 07 jari 401          int rowDim = means.getRowDimension();
2 26 Feb 07 jari 402
2 26 Feb 07 jari 403          Cluster gene_cluster = new Cluster();
2 26 Feb 07 jari 404          nodeList = gene_cluster.getNodeList();
2 26 Feb 07 jari 405          for (int i=0; i<reducedGeneSet.length; i++) {
2 26 Feb 07 jari 406              if (stop) {
2 26 Feb 07 jari 407                  throw new AbortException();
2 26 Feb 07 jari 408              }
2 26 Feb 07 jari 409              features = convert2int(reducedGeneSet[i]);
2 26 Feb 07 jari 410
2 26 Feb 07 jari 411              Node node = new Node(features);
2 26 Feb 07 jari 412              nodeList.addNode(node);
2 26 Feb 07 jari 413          }
2 26 Feb 07 jari 414
2 26 Feb 07 jari 415          Matrix means_used = getMeansForGenes(reducedGeneSet);
2 26 Feb 07 jari 416          Matrix variances_used = getVariancesForGenes(reducedGeneSet, means_used);
2 26 Feb 07 jari 417
2 26 Feb 07 jari 418          Matrix means_unused = getMeansForGenes(reducedGeneSet);
2 26 Feb 07 jari 419          Matrix variances_unused = getVariancesForGenes(reducedGeneSet, means_used);
2 26 Feb 07 jari 420
2 26 Feb 07 jari 421
2 26 Feb 07 jari 422          // construct the result
2 26 Feb 07 jari 423          AlgorithmData result = new AlgorithmData();
2 26 Feb 07 jari 424
2 26 Feb 07 jari 425          result.addParam("numberOfGenes", String.valueOf(numberOfGenes));
2 26 Feb 07 jari 426          result.addMatrix("probFunction", getFloatMatrix(probFunction));
2 26 Feb 07 jari 427
2 26 Feb 07 jari 428 /*
2 26 Feb 07 jari 429          System.out.println("DAM.java -- geneComponentMatrix ");
2 26 Feb 07 jari 430          printMatrix(geneComponentMatrix);
2 26 Feb 07 jari 431 */
2 26 Feb 07 jari 432
2 26 Feb 07 jari 433          FloatMatrix tempMatrix = getFloatMatrix(geneComponentMatrix).transpose();  
2 26 Feb 07 jari 434          FloatMatrix matrix3D;  
2 26 Feb 07 jari 435
2 26 Feb 07 jari 436          if(tempMatrix.getColumnDimension() == 2) {
2 26 Feb 07 jari 437              matrix3D = new FloatMatrix(tempMatrix.getRowDimension(), 3);
2 26 Feb 07 jari 438              matrix3D.setMatrix(0, tempMatrix.getRowDimension()-1, 0, tempMatrix.getColumnDimension()-1, tempMatrix);
2 26 Feb 07 jari 439              for (int i=0; i< tempMatrix.getRowDimension(); i++) {
2 26 Feb 07 jari 440                  matrix3D.set(i, 2, 0);
2 26 Feb 07 jari 441              }
2 26 Feb 07 jari 442          } else {
2 26 Feb 07 jari 443              matrix3D = tempMatrix;
2 26 Feb 07 jari 444          }
2 26 Feb 07 jari 445
2 26 Feb 07 jari 446 //         System.out.println("DAM.java -- matrix3D ");
2 26 Feb 07 jari 447 //         printMatrix(matrix3D);
2 26 Feb 07 jari 448  
2 26 Feb 07 jari 449          result.addMatrix("matrix3D", matrix3D);
2 26 Feb 07 jari 450
2 26 Feb 07 jari 451          result.addCluster("cluster", result_cluster);
2 26 Feb 07 jari 452          result.addCluster("geneCluster", gene_cluster);
2 26 Feb 07 jari 453
2 26 Feb 07 jari 454          result.addMatrix("clusters_means", getFloatMatrix(means));
2 26 Feb 07 jari 455          result.addMatrix("clusters_variances", getFloatMatrix(variances));
2 26 Feb 07 jari 456
2 26 Feb 07 jari 457          result.addMatrix("clusters_means_used", getFloatMatrix(means_used));
2 26 Feb 07 jari 458          result.addMatrix("clusters_variances_used", getFloatMatrix(variances_used));
2 26 Feb 07 jari 459
2 26 Feb 07 jari 460          result.addMatrix("clusters_means_unused", getFloatMatrix(means_unused));
2 26 Feb 07 jari 461          result.addMatrix("clusters_variances_unused", getFloatMatrix(variances_unused));
2 26 Feb 07 jari 462
2 26 Feb 07 jari 463 /*
2 26 Feb 07 jari 464          if (whichAlgorithm == 2) {
2 26 Feb 07 jari 465        for (int leaveOutSample=0; leaveOutSample<numberOfSamples; leaveOutSample++) {
2 26 Feb 07 jari 466
2 26 Feb 07 jari 467               System.out.println("DAM.java: Used Gene Index for sample "  + leaveOutSample); 
2 26 Feb 07 jari 468
2 26 Feb 07 jari 469            for(int i=0; i< reducedGeneSetForA2[leaveOutSample][used].size(); i++) {
2 26 Feb 07 jari 470          System.out.print(((Integer)(reducedGeneSetForA2[leaveOutSample][used].get(i))).intValue() + ", ");
2 26 Feb 07 jari 471                  }
2 26 Feb 07 jari 472      System.out.println(" ");
2 26 Feb 07 jari 473      System.out.println(" ");
2 26 Feb 07 jari 474
2 26 Feb 07 jari 475              System.out.println("DAM.java: Unsed Gene Index for sample "  + leaveOutSample); 
2 26 Feb 07 jari 476
2 26 Feb 07 jari 477            for(int i=0; i< reducedGeneSetForA2[leaveOutSample][unused].size(); i++) {
2 26 Feb 07 jari 478          System.out.print(((Integer)(reducedGeneSetForA2[leaveOutSample][unused].get(i))).intValue() + ", ");
2 26 Feb 07 jari 479                  }
2 26 Feb 07 jari 480      System.out.println(" ");
2 26 Feb 07 jari 481      System.out.println(" ");
2 26 Feb 07 jari 482        }
2 26 Feb 07 jari 483        System.out.println(" ");
2 26 Feb 07 jari 484
2 26 Feb 07 jari 485          } else {
2 26 Feb 07 jari 486        System.out.println("DAM.java: usedGeneIndices size: " + usedGeneIndices.length); 
2 26 Feb 07 jari 487        for(int i=0; i< usedGeneIndices.length; i++) {
2 26 Feb 07 jari 488      System.out.print(usedGeneIndices[i] + ", "); 
2 26 Feb 07 jari 489        }
2 26 Feb 07 jari 490        System.out.println(" ");
2 26 Feb 07 jari 491        System.out.println(" ");
2 26 Feb 07 jari 492        System.out.println("DAM.java: unusedGeneIndices size: " + unusedGeneIndices.length); 
2 26 Feb 07 jari 493        for(int i=0; i< unusedGeneIndices.length; i++) {
2 26 Feb 07 jari 494      System.out.print(unusedGeneIndices[i] + ", "); 
2 26 Feb 07 jari 495        }
2 26 Feb 07 jari 496        System.out.println(" ");
2 26 Feb 07 jari 497        System.out.println(" ");
2 26 Feb 07 jari 498          }
2 26 Feb 07 jari 499
2 26 Feb 07 jari 500 */
2 26 Feb 07 jari 501
2 26 Feb 07 jari 502          result.addIntArray("usedGeneIndices", usedGeneIndices);
2 26 Feb 07 jari 503          result.addIntArray("unusedGeneIndices", unusedGeneIndices);
2 26 Feb 07 jari 504
2 26 Feb 07 jari 505    return result;
2 26 Feb 07 jari 506     }
2 26 Feb 07 jari 507
2 26 Feb 07 jari 508
2 26 Feb 07 jari 509     
2 26 Feb 07 jari 510     //jcb, get a new responseMatrix from the initial result
2 26 Feb 07 jari 511     
2 26 Feb 07 jari 512     private Matrix extractResult(Matrix pMatrix) {
2 26 Feb 07 jari 513     
2 26 Feb 07 jari 514          Matrix newMatrix = new Matrix(numberOfSamples, numberOfClasses);
2 26 Feb 07 jari 515          int rows = pMatrix.getRowDimension();
2 26 Feb 07 jari 516          int cols = pMatrix.getColumnDimension();
2 26 Feb 07 jari 517          
2 26 Feb 07 jari 518         //set to zero
2 26 Feb 07 jari 519                  for (int row=0; row<rows; row++) {
2 26 Feb 07 jari 520              for (int column=0; column<cols; column++) {
2 26 Feb 07 jari 521            newMatrix.set(row, column, 0);
2 26 Feb 07 jari 522              }
2 26 Feb 07 jari 523          }
2 26 Feb 07 jari 524         
2 26 Feb 07 jari 525        double [][] array = pMatrix.getArray();
2 26 Feb 07 jari 526        int classID = 0; 
2 26 Feb 07 jari 527        
2 26 Feb 07 jari 528         for(int row = 0; row < rows; row++) {
2 26 Feb 07 jari 529             classID = getClassID(array[row]);
2 26 Feb 07 jari 530             newMatrix.set(row, classID, 1.0);
2 26 Feb 07 jari 531         }
2 26 Feb 07 jari 532         
2 26 Feb 07 jari 533         return newMatrix;
2 26 Feb 07 jari 534     }
2 26 Feb 07 jari 535     
2 26 Feb 07 jari 536     //jcb, returns classID for a given row in the prob matrix
2 26 Feb 07 jari 537     
2 26 Feb 07 jari 538     private int getClassID(double [] probs) {
2 26 Feb 07 jari 539         int hasMax = 0;
2 26 Feb 07 jari 540         double max = 0.0, currMax = 0.0;
2 26 Feb 07 jari 541         for(int i = 0; i < probs.length; i++) {
2 26 Feb 07 jari 542             currMax = Math.max(max, probs[i]);
2 26 Feb 07 jari 543             if(currMax > max) {
2 26 Feb 07 jari 544                 hasMax = i;
2 26 Feb 07 jari 545                 max = currMax;
2 26 Feb 07 jari 546             }
2 26 Feb 07 jari 547         }
2 26 Feb 07 jari 548         return hasMax;
2 26 Feb 07 jari 549     }
2 26 Feb 07 jari 550
2 26 Feb 07 jari 551     /**
2 26 Feb 07 jari 552      *  geneSelection(): 
2 26 Feb 07 jari 553      *  This method compares all (G+1, 2) pairwise (absolute) mean differences ( G+1 = numberOfClasses )
2 26 Feb 07 jari 554      *  |mean(Xk)-mean(Xk')| (k!=k') to a critical score
2 26 Feb 07 jari 555      *     t*sqrt(MSE(1/nk + 1/nk'))
2 26 Feb 07 jari 556      *  where MSE(mean square error) is the estimate of vairability from ANOVA model
2 26 Feb 07 jari 557      *  with one factor and G+1 groups, t is t(alpha/2, N-(G+1)) value of t-distribution.
2 26 Feb 07 jari 558      *  Each gene (j=1, 2, ... p) is ranked according to the number of times the pairwise 
2 26 Feb 07 jari 559      *  absolute mean difference exceeded the critical score.
2 26 Feb 07 jari 560      *
2 26 Feb 07 jari 561      *  Input: none
2 26 Feb 07 jari 562      *  Output: An array of rank vectors that contains the gene indices of that rank 
2 26 Feb 07 jari 563      */
2 26 Feb 07 jari 564     public Vector[] geneSelection(Matrix expMatrix) {
2 26 Feb 07 jari 565
2 26 Feb 07 jari 566          
2 26 Feb 07 jari 567 //   System.out.println(" ");
2 26 Feb 07 jari 568 //   System.out.println("**************** Begin geneSelection() *************** ");
2 26 Feb 07 jari 569 //   System.out.println(" ");
2 26 Feb 07 jari 570
2 26 Feb 07 jari 571          Vector[] geneIndices;
2 26 Feb 07 jari 572
2 26 Feb 07 jari 573    double tValue=0;
2 26 Feb 07 jari 574    double meanSquareError=0;
2 26 Feb 07 jari 575
2 26 Feb 07 jari 576          int numOfGenes = expMatrix.getRowDimension();
2 26 Feb 07 jari 577          int numOfSamples = expMatrix.getColumnDimension();
2 26 Feb 07 jari 578
2 26 Feb 07 jari 579    double [][] criticalScores = new double[numberOfClasses][numberOfClasses];
2 26 Feb 07 jari 580          double [][] meanDifferences = new double[numberOfClasses][numberOfClasses];
2 26 Feb 07 jari 581
2 26 Feb 07 jari 582    // Obtain the t value for the given alpha value and numberOfClasses
2 26 Feb 07 jari 583
2 26 Feb 07 jari 584          AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, numberOfGenes, "Gene Screening\n");
2 26 Feb 07 jari 585          fireValueChanged(event);
2 26 Feb 07 jari 586          event.setId(AlgorithmEvent.PROGRESS_VALUE);            
2 26 Feb 07 jari 587
2 26 Feb 07 jari 588  /*        System.out.println("");
2 26 Feb 07 jari 589          System.out.println("geneSelection() - numOfGenes = " + numOfGenes);
2 26 Feb 07 jari 590          System.out.println("geneSelection() - numOfSamples = " + numOfSamples);
2 26 Feb 07 jari 591          System.out.println("geneSelection() - numberOfClasses = " + numberOfClasses);
2 26 Feb 07 jari 592          System.out.println("geneSelection() - alpha = " + alpha);
2 26 Feb 07 jari 593 */
2 26 Feb 07 jari 594          TDistribution testT = new TDistribution(numOfSamples-(numberOfClasses+1));
2 26 Feb 07 jari 595    tValue = -testT.inverse(alpha/2.0d);
2 26 Feb 07 jari 596
2 26 Feb 07 jari 597    double tValue1 = getTValue(numOfSamples-(numberOfClasses+1), alpha);
2 26 Feb 07 jari 598
2 26 Feb 07 jari 599      //    System.out.println("geneSelection() - tValue = " + tValue);
2 26 Feb 07 jari 600      //    System.out.println("geneSelection() - tValue1 = " + tValue1);
2 26 Feb 07 jari 601
2 26 Feb 07 jari 602          if (tValue < - 1.0) {
2 26 Feb 07 jari 603              return null;
2 26 Feb 07 jari 604          }
2 26 Feb 07 jari 605
2 26 Feb 07 jari 606     //     System.out.println("");
2 26 Feb 07 jari 607     //     System.out.println("geneSelection() - tValue = " + tValue);
2 26 Feb 07 jari 608
2 26 Feb 07 jari 609    // Obtain the Mean Square Error Estimate
2 26 Feb 07 jari 610    meanSquareError = getSampleVariance(expMatrix, getSampleMeans(expMatrix));
2 26 Feb 07 jari 611      //    System.out.println("geneSelection() - meanSquareError = " + meanSquareError);
2 26 Feb 07 jari 612  
2 26 Feb 07 jari 613          // For each gene
2 26 Feb 07 jari 614          for (int gene=0; gene < numOfGenes; gene++) {
2 26 Feb 07 jari 615
2 26 Feb 07 jari 616        for(int class1=0; class1<numberOfClasses; class1++) {
2 26 Feb 07 jari 617
2 26 Feb 07 jari 618      for(int class2=class1+1; class2<numberOfClasses; class2++) {
2 26 Feb 07 jari 619
2 26 Feb 07 jari 620              // Calculate the critical score for each pair of classes among all classes
2 26 Feb 07 jari 621        criticalScores[class1][class2] = tValue * 
2 26 Feb 07 jari 622            Math.sqrt(meanSquareError*(1.0/classSampleArray[class1] + 1.0/classSampleArray[class2]));
2 26 Feb 07 jari 623
2 26 Feb 07 jari 624              // Calculate the pairwise absolute mean differneces for all classes 
2 26 Feb 07 jari 625        meanDifferences[class1][class2] = Math.abs(
2 26 Feb 07 jari 626                           classExpSumMatrix.get(gene, class1)/classSampleArray[class1]-
2 26 Feb 07 jari 627                           classExpSumMatrix.get(gene, class2)/classSampleArray[class2]); 
2 26 Feb 07 jari 628
2 26 Feb 07 jari 629                    // rank the genes
2 26 Feb 07 jari 630                    if (meanDifferences[class1][class2] > criticalScores[class1][class2]) {
2 26 Feb 07 jari 631                        geneRank[gene] ++ ;
2 26 Feb 07 jari 632                    }
2 26 Feb 07 jari 633                 }
2 26 Feb 07 jari 634        }
2 26 Feb 07 jari 635          }
2 26 Feb 07 jari 636
2 26 Feb 07 jari 637          // Find highest rank
2 26 Feb 07 jari 638          highestGeneRank = 0;
2 26 Feb 07 jari 639
2 26 Feb 07 jari 640          for (int gene=0; gene<numOfGenes; gene++) {
2 26 Feb 07 jari 641             if (geneRank[gene] > highestGeneRank) {
2 26 Feb 07 jari 642                highestGeneRank = geneRank[gene];
2 26 Feb 07 jari 643             }
2 26 Feb 07 jari 644          }
2 26 Feb 07 jari 645
2 26 Feb 07 jari 646     //     System.out.println("");
2 26 Feb 07 jari 647     //     System.out.println("geneSelection() - highestGeneRank = " + highestGeneRank);
2 26 Feb 07 jari 648
2 26 Feb 07 jari 649          // Calculate the rank vector for storing gene indices
2 26 Feb 07 jari 650          geneIndices = new Vector[highestGeneRank+1]; 
2 26 Feb 07 jari 651
2 26 Feb 07 jari 652
2 26 Feb 07 jari 653          for (int rank=0; rank<highestGeneRank+1; rank++) {
2 26 Feb 07 jari 654             geneIndices[rank] = new Vector(); 
2 26 Feb 07 jari 655             for (int gene=0; gene<numOfGenes; gene++) {
2 26 Feb 07 jari 656
2 26 Feb 07 jari 657                 if (geneRank[gene] == rank) {
2 26 Feb 07 jari 658                     geneIndices[rank].add(new Integer(gene));
2 26 Feb 07 jari 659                     
2 26 Feb 07 jari 660                     if (rank > 0) {
2 26 Feb 07 jari 661                         geneIndices[0].remove(new Integer(gene));
2 26 Feb 07 jari 662                     }
2 26 Feb 07 jari 663                 }
2 26 Feb 07 jari 664
2 26 Feb 07 jari 665             }
2 26 Feb 07 jari 666          }
2 26 Feb 07 jari 667
2 26 Feb 07 jari 668      /*    for (int rank=0; rank<highestGeneRank+1; rank++) {
2 26 Feb 07 jari 669              System.out.println("geneSelection() - rank = " + rank + ", geneIndices[rank].size()=" + geneIndices[rank].size());
2 26 Feb 07 jari 670          }
2 26 Feb 07 jari 671          System.out.println("");
2 26 Feb 07 jari 672
2 26 Feb 07 jari 673    System.out.println(" ");
2 26 Feb 07 jari 674    System.out.println("**************** End geneSelection() *************** ");
2 26 Feb 07 jari 675    System.out.println(" ");
2 26 Feb 07 jari 676 */
2 26 Feb 07 jari 677          return geneIndices;
2 26 Feb 07 jari 678
2 26 Feb 07 jari 679     } // end of geneSelection
2 26 Feb 07 jari 680
2 26 Feb 07 jari 681
2 26 Feb 07 jari 682
2 26 Feb 07 jari 683     /**
2 26 Feb 07 jari 684      * MLE Algorithm:
2 26 Feb 07 jari 685      *     qualitative response variable: y = k where k = 0, 1, 2, ... G
2 26 Feb 07 jari 686      *     gene expression levels from one experiment: x = (x1, x2, ... xp) 
2 26 Feb 07 jari 687      *     probability: P(y = k | x) = exp(gk(x))/(1+exp(g0(x)) + exp(g1(x)) + ... exp(gK(x))) 
2 26 Feb 07 jari 688      *     where gk(x) = betak0*x0 + betak1*x1 +...+ betakp*xp
2 26 Feb 07 jari 689      *     beta = (beta1, .... betak) with betak= (betak0, betak1,..., betakp);
2 26 Feb 07 jari 690      *     bata can be estimated from MLE when numberOfSamples > numberOfClasses * numberOfGenes.
2 26 Feb 07 jari 691      * 
2 26 Feb 07 jari 692      *     beta_Matrix: beta Matrix of size [numberOfClasses*(numberOfGenes+1), 1];
2 26 Feb 07 jari 693      *     Z_Matrix: Indicator Matrix of size [numberOfSamples, (numberOfClasses+1)]
2 26 Feb 07 jari 694      *     P_Matrix: Probability Matrix of size [numberOfGenes, (numberOfClasses+1)]
2 26 Feb 07 jari 695      *     I_Matrix: Information Matrix of size 
2 26 Feb 07 jari 696      *               [numberOfClasses*(numberOfGenes+1), numberOfClasses*(numberOfGenes+1)]
2 26 Feb 07 jari 697      *     S_Matrix: Score Matrix of size [numberOfClasses*(numberOfGenes+1), 1]
2 26 Feb 07 jari 698      *     W_Matrices[numberOfClasses*numberOfClasses]: 
2 26 Feb 07 jari 699      *           An array of diagonal matrix of size [numberOfSamples, numberOfSamples]
2 26 Feb 07 jari 700      *
2 26 Feb 07 jari 701      *     S(beta) = (delta(1)(x(1) - X'w(1))) + ... + (delta(N)(x(N) - X'w(N))
2 26 Feb 07 jari 702      *     Information Matrix:
2 26 Feb 07 jari 703      *     I(beta) = delta(1)[X'w(1)X - (X'w(1)(w(1)'X)] + ... +
2 26 Feb 07 jari 704      *     delta(1)[X'w(N)X - (X'w(N)(w(N)'X)] 
2 26 Feb 07 jari 705      * 
2 26 Feb 07 jari 706      *     NewTon-Raphson algorithm:
2 26 Feb 07 jari 707      *       beta(s+1) = beta(s) + inverse(I(beta))*(beta(s)*S(beta(s)))
2 26 Feb 07 jari 708      * 
2 26 Feb 07 jari 709      *     If the NewTon-Raphson algorithm converges, the vector of coefficients is 
2 26 Feb 07 jari 710      *     the maximum partial likelihood estimate of beta.
2 26 Feb 07 jari 711      */                     
2 26 Feb 07 jari 712     public Matrix [] mleAlgorithm(Matrix compMatrix, Matrix compResponseMatrix) throws AlgorithmException{
2 26 Feb 07 jari 713
2 26 Feb 07 jari 714 /*  System.out.println(" ");
2 26 Feb 07 jari 715    System.out.println("******************** Begin mleAlgorithm() ********************");
2 26 Feb 07 jari 716   System.out.println(" ");
2 26 Feb 07 jari 717
2 26 Feb 07 jari 718 */
2 26 Feb 07 jari 719         int numOfGenes = compMatrix.getRowDimension();
2 26 Feb 07 jari 720         int numOfSamples = compMatrix.getColumnDimension();
2 26 Feb 07 jari 721 /*
2 26 Feb 07 jari 722         System.out.println(" ");
2 26 Feb 07 jari 723         System.out.println("mleAlgorithm() - numOfGenes = " + numOfGenes);
2 26 Feb 07 jari 724         System.out.println("mleAlgorithm() - numOfSamples = " + numOfSamples);
2 26 Feb 07 jari 725 */
2 26 Feb 07 jari 726
2 26 Feb 07 jari 727         AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, numberOfGenes, "MLE Algorithm \n");
2 26 Feb 07 jari 728         fireValueChanged(event);
2 26 Feb 07 jari 729         event.setId(AlgorithmEvent.PROGRESS_VALUE);            
2 26 Feb 07 jari 730
2 26 Feb 07 jari 731         Matrix [] X = new Matrix[numOfSamples]; // an array of sample matrices
2 26 Feb 07 jari 732         for (int i=0; i<numOfSamples; i++) {
2 26 Feb 07 jari 733             X[i] = new Matrix(numOfGenes+1, 1);
2 26 Feb 07 jari 734             X[i].set(0, 0, 1);
2 26 Feb 07 jari 735             X[i].setMatrix(1, numOfGenes, 0, 0, compMatrix.getMatrix(0, numOfGenes-1, i, i));
2 26 Feb 07 jari 736         }
2 26 Feb 07 jari 737
2 26 Feb 07 jari 738         // X_Matrix: obtained from compMatrix by adding 1's in the first row
2 26 Feb 07 jari 739         Matrix X_Matrix = new Matrix(numOfGenes+1, numOfSamples);
2 26 Feb 07 jari 740         for (int i=0; i<numOfSamples; i++) {
2 26 Feb 07 jari 741             X_Matrix.set(0, i, 1);
2 26 Feb 07 jari 742         }
2 26 Feb 07 jari 743         X_Matrix.setMatrix(1, numOfGenes, 0, numOfSamples-1, compMatrix);
2 26 Feb 07 jari 744
2 26 Feb 07 jari 745 /*
2 26 Feb 07 jari 746   System.out.println(" ");
2 26 Feb 07 jari 747   System.out.println("mleAlogrithm() - X_Matrix : ");
2 26 Feb 07 jari 748   for (int i=0; i<X_Matrix.getRowDimension(); i++) {
2 26 Feb 07 jari 749       for (int j=0; j<X_Matrix.getColumnDimension(); j++) {
2 26 Feb 07 jari 750           System.out.print(X_Matrix.get(i,j) + ", ");
2 26 Feb 07 jari 751       }
2 26 Feb 07 jari 752       System.out.println(" ");
2 26 Feb 07 jari 753       System.out.println(" ");
2 26 Feb 07 jari 754   }
2 26 Feb 07 jari 755 */
2 26 Feb 07 jari 756
2 26 Feb 07 jari 757         Matrix beta_Matrix = new Matrix(numberOfClasses*(numOfGenes+1), 1);
2 26 Feb 07 jari 758         Matrix new_beta_Matrix = new Matrix(numberOfClasses*(numOfGenes+1), 1);
2 26 Feb 07 jari 759
2 26 Feb 07 jari 760         Matrix [] oldBeta = new Matrix[numberOfClasses];
2 26 Feb 07 jari 761         Matrix [] newBeta = new Matrix[numberOfClasses];
2 26 Feb 07 jari 762         for (int k=0; k< numberOfClasses; k++) {
2 26 Feb 07 jari 763             oldBeta[k] = new Matrix(numOfGenes+1, 1);
2 26 Feb 07 jari 764             newBeta[k] = new Matrix(numOfGenes+1, 1);
2 26 Feb 07 jari 765         }
2 26 Feb 07 jari 766
2 26 Feb 07 jari 767         // Probability Matrix (numOfSamples x numberOfClasses )
2 26 Feb 07 jari 768         Matrix P_Matrix = new Matrix(numOfSamples, numberOfClasses); 
2 26 Feb 07 jari 769
2 26 Feb 07 jari 770         // Score Matrix (numOfSamples*(numberOfGenes+1) x 1 )
2 26 Feb 07 jari 771         Matrix S_Matrix = new Matrix(numberOfClasses*(numOfGenes+1), 1); 
2 26 Feb 07 jari 772
2 26 Feb 07 jari 773         // A Matrix of Diagonal Matrices
2 26 Feb 07 jari 774         Matrix [][] W_Matrices = new Matrix[numberOfClasses][numberOfClasses]; 
2 26 Feb 07 jari 775
2 26 Feb 07 jari 776         for (int i=0; i< numberOfClasses; i++) {
2 26 Feb 07 jari 777             for (int j=0; j< numberOfClasses; j++) {
2 26 Feb 07 jari 778                 W_Matrices[i][j] = new Matrix(numOfSamples,numOfSamples);
2 26 Feb 07 jari 779             }
2 26 Feb 07 jari 780         }
2 26 Feb 07 jari 781
2 26 Feb 07 jari 782         // Information Matrix (numberOfClasses*(numOfGenes+1), numberOfClasses*(numOfGenes+1));
2 26 Feb 07 jari 783         Matrix I_Matrix = new Matrix(numberOfClasses*(numOfGenes+1), numberOfClasses*(numOfGenes+1));
2 26 Feb 07 jari 784
2 26 Feb 07 jari 785         Matrix [][] I = new Matrix[numberOfClasses][numberOfClasses];
2 26 Feb 07 jari 786         for (int i=0; i<numberOfClasses; i++) {
2 26 Feb 07 jari 787             for (int j=0; j<numberOfClasses; j++) {
2 26 Feb 07 jari 788                 I[i][j] = new Matrix (numOfGenes+1, numOfGenes+1);
2 26 Feb 07 jari 789             }
2 26 Feb 07 jari 790         }
2 26 Feb 07 jari 791
2 26 Feb 07 jari 792
2 26 Feb 07 jari 793         // Initialize oldBeta and newBeta variables 
2 26 Feb 07 jari 794         for (int k=0; k<numberOfClasses; k++) {
2 26 Feb 07 jari 795             for(int row=0; row<numOfGenes+1; row++) {
2 26 Feb 07 jari 796                 oldBeta[k].set(row, 0, 0.1f);
2 26 Feb 07 jari 797                 newBeta[k].set(row, 0, 0.1f);
2 26 Feb 07 jari 798             }
2 26 Feb 07 jari 799         }
2 26 Feb 07 jari 800
2 26 Feb 07 jari 801         // Initialize beta_Matrix
2 26 Feb 07 jari 802         for (int i=0; i<numberOfClasses; i++) {
2 26 Feb 07 jari 803             for(int j=0; j<numOfGenes+1; j++) {
2 26 Feb 07 jari 804                 beta_Matrix.set(i*j, 0, oldBeta[i].get(j,0));
2 26 Feb 07 jari 805                 new_beta_Matrix.set(i*j, 0, oldBeta[i].get(j,0));
2 26 Feb 07 jari 806             }
2 26 Feb 07 jari 807         }
2 26 Feb 07 jari 808
2 26 Feb 07 jari 809         // Obtain the Z_matrix: Indicator Matrix ( numOfSamples x numberOfClasses )
2 26 Feb 07 jari 810         Matrix Z_Matrix = compResponseMatrix.getMatrix(0, numOfSamples-1, 0, numberOfClasses-1); 
2 26 Feb 07 jari 811
2 26 Feb 07 jari 812 /*
2 26 Feb 07 jari 813   System.out.println(" ");
2 26 Feb 07 jari 814   System.out.println("mleAlogrithm() - responseMatrix : ");
2 26 Feb 07 jari 815   for (int i=0; i<responseMatrix.getRowDimension(); i++) {
2 26 Feb 07 jari 816       for (int j=0; j<responseMatrix.getColumnDimension(); j++) {
2 26 Feb 07 jari 817           System.out.print(responseMatrix.get(i,j) + ", ");
2 26 Feb 07 jari 818       }
2 26 Feb 07 jari 819       System.out.println(" ");
2 26 Feb 07 jari 820       System.out.println(" ");
2 26 Feb 07 jari 821   }
2 26 Feb 07 jari 822
2 26 Feb 07 jari 823   System.out.println(" ");
2 26 Feb 07 jari 824   System.out.println("mleAlogrithm() - Z_Matrix : ");
2 26 Feb 07 jari 825   for (int i=0; i<Z_Matrix.getRowDimension(); i++) {
2 26 Feb 07 jari 826       for (int j=0; j<Z_Matrix.getColumnDimension(); j++) {
2 26 Feb 07 jari 827           System.out.print(Z_Matrix.get(i,j) + ", ");
2 26 Feb 07 jari 828       }
2 26 Feb 07 jari 829       System.out.println(" ");
2 26 Feb 07 jari 830       System.out.println(" ");
2 26 Feb 07 jari 831   }
2 26 Feb 07 jari 832 */
2 26 Feb 07 jari 833
2 26 Feb 07 jari 834         double g=0.0, c=0.0, localSum=0.0;
2 26 Feb 07 jari 835         
2 26 Feb 07 jari 836         // Main iteration for obtain beta estimator
2 26 Feb 07 jari 837         for (int iteration = 0; iteration < 40; iteration++) {
2 26 Feb 07 jari 838
2 26 Feb 07 jari 839             oldBeta = newBeta;
2 26 Feb 07 jari 840
2 26 Feb 07 jari 841 /*
2 26 Feb 07 jari 842       System.out.println(" ");
2 26 Feb 07 jari 843       System.out.println("mleAlogrithm() - oldBeta[i] : ");
2 26 Feb 07 jari 844       for (int i=0; i<numberOfClasses; i++) {
2 26 Feb 07 jari 845     for (int j=0; j<numOfGenes+1; j++) {
2 26 Feb 07 jari 846         System.out.print(oldBeta[i].get(j,0) + ", ");
2 26 Feb 07 jari 847     }
2 26 Feb 07 jari 848     System.out.println(" ");
2 26 Feb 07 jari 849     System.out.println(" ");
2 26 Feb 07 jari 850       }
2 26 Feb 07 jari 851 */
2 26 Feb 07 jari 852
2 26 Feb 07 jari 853             // Calculate P_Matrix 
2 26 Feb 07 jari 854       for (int i=0; i<numOfSamples; i++) {
2 26 Feb 07 jari 855     for (int j=0; j<numberOfClasses; j++) {
2 26 Feb 07 jari 856         localSum += Math.exp(((X[i].transpose()).times(oldBeta[j])).get(0,0));
2 26 Feb 07 jari 857     }
2 26 Feb 07 jari 858
2 26 Feb 07 jari 859     c = Math.log(1 + localSum);
2 26 Feb 07 jari 860
2 26 Feb 07 jari 861     for (int k=0; k<numberOfClasses; k++) {
2 26 Feb 07 jari 862
2 26 Feb 07 jari 863         g = ((X[i].transpose()).times(oldBeta[k])).get(0,0);
2 26 Feb 07 jari 864
2 26 Feb 07 jari 865         P_Matrix.set(i, k, Math.exp(g-c)); 
2 26 Feb 07 jari 866     }
2 26 Feb 07 jari 867
2 26 Feb 07 jari 868                 localSum = 0.0;
2 26 Feb 07 jari 869       }
2 26 Feb 07 jari 870
2 26 Feb 07 jari 871 /*
2 26 Feb 07 jari 872       System.out.println(" ");
2 26 Feb 07 jari 873       System.out.println("mleAlogrithm() - P_Matrix : ");
2 26 Feb 07 jari 874       System.out.println("mleAlogrithm() - P_Matrix : Row Dimension = " + P_Matrix.getRowDimension());
2 26 Feb 07 jari 875       System.out.println("mleAlogrithm() - P_Matrix : Column Dimension = " + P_Matrix.getColumnDimension());
2 26 Feb 07 jari 876       for (int i=0; i<P_Matrix.getRowDimension(); i++) {
2 26 Feb 07 jari 877     for (int j=0; j<P_Matrix.getColumnDimension(); j++) {
2 26 Feb 07 jari 878         System.out.print(P_Matrix.get(i,j) + ", ");
2 26 Feb 07 jari 879     }
2 26 Feb 07 jari 880     System.out.println(" ");
2 26 Feb 07 jari 881     System.out.println(" ");
2 26 Feb 07 jari 882       }
2 26 Feb 07 jari 883 */
2 26 Feb 07 jari 884
2 26 Feb 07 jari 885             // Calculate local sum matrix for sub score matrix 
2 26 Feb 07 jari 886             Matrix [] sum = new Matrix[numberOfClasses];
2 26 Feb 07 jari 887             for (int k=0; k<numberOfClasses; k++) {
2 26 Feb 07 jari 888                 sum[k] = new Matrix(numOfGenes+1, 1);
2 26 Feb 07 jari 889
2 26 Feb 07 jari 890                 for (int i=0; i<numOfSamples; i++) {
2 26 Feb 07 jari 891                     sum[k] = sum[k].plus(X[i].times(Z_Matrix.get(i,k) - P_Matrix.get(i,k)));
2 26 Feb 07 jari 892                 }
2 26 Feb 07 jari 893             }
2 26 Feb 07 jari 894
2 26 Feb 07 jari 895             // Calculate S_Matrix
2 26 Feb 07 jari 896             for (int k=0; k<numberOfClasses; k++) {
2 26 Feb 07 jari 897                 S_Matrix.setMatrix(k*(numOfGenes+1), (k+1)*(numOfGenes+1)-1, 0, 0, sum[k]);
2 26 Feb 07 jari 898             }
2 26 Feb 07 jari 899
2 26 Feb 07 jari 900             // Calculate W_Matrix
2 26 Feb 07 jari 901             for (int i=0; i< numberOfClasses; i++) {
2 26 Feb 07 jari 902                 for (int j=0; j< numberOfClasses; j++) {
2 26 Feb 07 jari 903                     if (i==j) {
2 26 Feb 07 jari 904                        for (int k=0; k<numOfSamples; k++) {
2 26 Feb 07 jari 905                            W_Matrices[i][j].set(k, k, P_Matrix.get(k,i)*(1-P_Matrix.get(k,i)));
2 26 Feb 07 jari 906                        }
2 26 Feb 07 jari 907                     } else {
2 26 Feb 07 jari 908                        for (int k=0; k<numOfSamples; k++) {
2 26 Feb 07 jari 909                            W_Matrices[i][j].set(k, k, P_Matrix.get(k,i)*P_Matrix.get(k,j));
2 26 Feb 07 jari 910                        }
2 26 Feb 07 jari 911                     }
2 26 Feb 07 jari 912                 }
2 26 Feb 07 jari 913             }
2 26 Feb 07 jari 914
2 26 Feb 07 jari 915             // Calculate sub I matrices
2 26 Feb 07 jari 916             for (int k1=0; k1<numberOfClasses; k1++) {
2 26 Feb 07 jari 917          for (int k2=0; k2<numberOfClasses; k2++) {
2 26 Feb 07 jari 918                   if (k1==k2) {
2 26 Feb 07 jari 919                       I[k1][k2] = X_Matrix.times(W_Matrices[k1][k2]).times(X_Matrix.transpose());
2 26 Feb 07 jari 920                   } else {
2 26 Feb 07 jari 921                       I[k1][k2] = (X_Matrix.times(W_Matrices[k1][k2]).times(X_Matrix.transpose())).uminus();
2 26 Feb 07 jari 922                   }
2 26 Feb 07 jari 923          }
2 26 Feb 07 jari 924             }
2 26 Feb 07 jari 925              
2 26 Feb 07 jari 926 /*
2 26 Feb 07 jari 927       System.out.println("mleAlogrithm() - I[0][0] : ");
2 26 Feb 07 jari 928       for (int i=0; i<I[0][0].getRowDimension(); i++) {
2 26 Feb 07 jari 929     for (int j=0; j<I[0][0].getColumnDimension(); j++) {
2 26 Feb 07 jari 930          System.out.print(I[0][0].get(i,j) + ", ");
2 26 Feb 07 jari 931     }
2 26 Feb 07 jari 932     System.out.println(" ");
2 26 Feb 07 jari 933       }
2 26 Feb 07 jari 934         System.out.println(" ");
2 26 Feb 07 jari 935 */
2 26 Feb 07 jari 936
2 26 Feb 07 jari 937             // Calculate I_Matrix
2 26 Feb 07 jari 938             for (int row=0; row<numberOfClasses; row++) {
2 26 Feb 07 jari 939     for (int column=0; column<numberOfClasses; column++) {
2 26 Feb 07 jari 940                     I_Matrix.setMatrix(row*(numOfGenes+1), (row+1)*(numOfGenes+1)-1, column*(numOfGenes+1), (column+1)*(numOfGenes+1)-1, I[row][column]);
2 26 Feb 07 jari 941     }
2 26 Feb 07 jari 942             }
2 26 Feb 07 jari 943
2 26 Feb 07 jari 944 /*
2 26 Feb 07 jari 945             System.out.println("mleAlgorithm() - I_Matrix : ");
2 26 Feb 07 jari 946             for (int i=0; i<I_Matrix.getRowDimension(); i++) {
2 26 Feb 07 jari 947                 for (int j=0; j<I_Matrix.getColumnDimension(); j++) {
2 26 Feb 07 jari 948                      System.out.print(I_Matrix.get(i,j) + ", ");
2 26 Feb 07 jari 949                 }
2 26 Feb 07 jari 950                 System.out.println(" ");
2 26 Feb 07 jari 951             }
2 26 Feb 07 jari 952             System.out.println(" ");
2 26 Feb 07 jari 953 */
2 26 Feb 07 jari 954
2 26 Feb 07 jari 955
2 26 Feb 07 jari 956             // Calculate newBeta
2 26 Feb 07 jari 957             for (int i=0; i<numberOfClasses; i++) {
2 26 Feb 07 jari 958                 try {
2 26 Feb 07 jari 959                     Matrix inverseMatrix = I_Matrix.inverse();
2 26 Feb 07 jari 960                     newBeta[i] = oldBeta[i].plus(((I_Matrix.inverse()).times(S_Matrix)).getMatrix(i*(numOfGenes+1), (i+1)*(numOfGenes+1)-1, 0, 0));
2 26 Feb 07 jari 961                 } catch (RuntimeException ex) {
2 26 Feb 07 jari 962                     JOptionPane.showMessageDialog(null, "MLE Algorithm: Information Matrix is Singular", "Alert", JOptionPane.WARNING_MESSAGE);
2 26 Feb 07 jari 963 //                    throw new AbortException();
2 26 Feb 07 jari 964                     throw new AlgorithmException("MLE Algorithm: Singular Matrix");
2 26 Feb 07 jari 965                 }
2 26 Feb 07 jari 966             }
2 26 Feb 07 jari 967
2 26 Feb 07 jari 968
2 26 Feb 07 jari 969 /*
2 26 Feb 07 jari 970       System.out.println(" ");
2 26 Feb 07 jari 971       System.out.println("mleAlogrithm() - iteration = " + iteration);
2 26 Feb 07 jari 972       System.out.println(" ");
2 26 Feb 07 jari 973       System.out.println("mleAlogrithm() - newBeta[i] : ");
2 26 Feb 07 jari 974       for (int i=0; i<numberOfClasses; i++) {
2 26 Feb 07 jari 975     for (int j=0; j<numOfGenes+1; j++) {
2 26 Feb 07 jari 976         System.out.print(newBeta[i].get(j,0) + ", ");
2 26 Feb 07 jari 977     }
2 26 Feb 07 jari 978     System.out.println(" ");
2 26 Feb 07 jari 979     System.out.println(" ");
2 26 Feb 07 jari 980       }
2 26 Feb 07 jari 981 */
2 26 Feb 07 jari 982
2 26 Feb 07 jari 983         }
2 26 Feb 07 jari 984 /*
2 26 Feb 07 jari 985   System.out.println(" ");
2 26 Feb 07 jari 986    System.out.println("******************** End mleAlgorithm() ********************");
2 26 Feb 07 jari 987   System.out.println(" ");
2 26 Feb 07 jari 988 */
2 26 Feb 07 jari 989         return newBeta;
2 26 Feb 07 jari 990
2 26 Feb 07 jari 991     } // end of mleAlgorithm()
2 26 Feb 07 jari 992
2 26 Feb 07 jari 993
2 26 Feb 07 jari 994
2 26 Feb 07 jari 995     /**
2 26 Feb 07 jari 996      * PDA Algorithm:
2 26 Feb 07 jari 997      * qualitative response variable: y = k where k = 0, 1, 2, ... G
2 26 Feb 07 jari 998      * gene expression levels from one experiment: x = (x1, x2, ... xp) 
2 26 Feb 07 jari 999      * probability: P(y = k | x) = exp(gk(x))/(1+exp(g0(x)) + exp(g1(x)) + ... exp(gK(x))) 
2 26 Feb 07 jari 1000      * where gk(x) = betak0*x0 + betak1*x1 +...+ betakp*xp
2 26 Feb 07 jari 1001      * beta = (beta1, .... betak) with betak= (betak0, betak1,..., betakp);
2 26 Feb 07 jari 1002      * bata can be estimated from MLE when 
2 26 Feb 07 jari 1003      * numberOfSamples > numberOfClasses * numberOfGenes.
2 26 Feb 07 jari 1004      * 
2 26 Feb 07 jari 1005      * D matrix: N x N Risk Indicator Matrix
2 26 Feb 07 jari 1006      * u vector: N x 1
2 26 Feb 07 jari 1007      * a vector: N x 1
2 26 Feb 07 jari 1008      * W matrix: N x N matrix 
2 26 Feb 07 jari 1009      * W(i)(j) = (1/a(i))u(j)*D(i)(j)
2 26 Feb 07 jari 1010      * Wdiag[i] array of NxN matrices 
2 26 Feb 07 jari 1011      * delta vector: N x 1, Censoring Indicator
2 26 Feb 07 jari 1012      * delta(i) = I(t(i)= min(y(i), z(i)));
2 26 Feb 07 jari 1013      * 
2 26 Feb 07 jari 1014      * Log Partial likelihood: 
2 26 Feb 07 jari 1015      * L(beta) = beta' * X' * delta - delta'* a
2 26 Feb 07 jari 1016      * Score Vector:
2 26 Feb 07 jari 1017      * S(beta) = (delta(1)(x(1) - X'W(1))) + ... + (delta(N)(x(N) - X'W(N))
2 26 Feb 07 jari 1018      * Information Matrix:
2 26 Feb 07 jari 1019      * I(beta) = delta(1)[X'W(1)X - (X'W(1)(W(1)'X)] + ... +
2 26 Feb 07 jari 1020      * delta(1)[X'W(N)X - (X'W(N)(W(N)'X)] 
2 26 Feb 07 jari 1021      *
2 26 Feb 07 jari 1022      * NewTon-Raphson algorithm:
2 26 Feb 07 jari 1023      *    beta(s+1) = beta(s) + inverse(I(beta))*(beta(s)*S(beta(s)))
2 26 Feb 07 jari 1024      *
2 26 Feb 07 jari 1025      * If the NewTon-Raphson algorithm converges, the vector of coefficients is 
2 26 Feb 07 jari 1026      * the maximum partial likelihood estimate of beta.
2 26 Feb 07 jari 1027      */                 
2 26 Feb 07 jari 1028     public double[] pdaAlgorithm(Matrix testMatrix) {
2 26 Feb 07 jari 1029 /*
2 26 Feb 07 jari 1030   System.out.println(" ");
2 26 Feb 07 jari 1031   System.out.println("******************** Begin pdaAlgorithm() ********************");
2 26 Feb 07 jari 1032   System.out.println(" ");
2 26 Feb 07 jari 1033 */
2 26 Feb 07 jari 1034         int numOfGenes = testMatrix.getRowDimension();
2 26 Feb 07 jari 1035         int numOfTestSamples = testMatrix.getColumnDimension();
2 26 Feb 07 jari 1036       
2 26 Feb 07 jari 1037         Matrix X_Matrix = new Matrix(numOfGenes+1, numOfTestSamples);
2 26 Feb 07 jari 1038         
2 26 Feb 07 jari 1039         for (int i=0; i<numOfTestSamples; i++) {
2 26 Feb 07 jari 1040             X_Matrix.set(0, i, 1);
2 26 Feb 07 jari 1041         }
2 26 Feb 07 jari 1042         X_Matrix.setMatrix(1, numOfGenes, 0, numOfTestSamples-1, testMatrix);
2 26 Feb 07 jari 1043
2 26 Feb 07 jari 1044   double [] probFunction = new double [numberOfClasses];
2 26 Feb 07 jari 1045
2 26 Feb 07 jari 1046 /*
2 26 Feb 07 jari 1047   System.out.println("pdaAlgorithm() -- testMatrix: ");
2 26 Feb 07 jari 1048   printMatrix(testMatrix);
2 26 Feb 07 jari 1049
2 26 Feb 07 jari 1050   System.out.println(" ");
2 26 Feb 07 jari 1051   System.out.println("pdaAlgorithm() -- X_Matrix");
2 26 Feb 07 jari 1052   for (int row = 0; row < X_Matrix.getRowDimension(); row++) {
2 26 Feb 07 jari 1053        for (int column = 0; column < testMatrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 1054      System.out.print(X_Matrix.get(row, column) + ", ");
2 26 Feb 07 jari 1055        }
2 26 Feb 07 jari 1056        System.out.println(" "); 
2 26 Feb 07 jari 1057   }
2 26 Feb 07 jari 1058
2 26 Feb 07 jari 1059   System.out.println(" ");
2 26 Feb 07 jari 1060         for(int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1061       System.out.println("pdaAlgorithm() -- beta[" + classId + "]");
2 26 Feb 07 jari 1062       for (int row = 0; row < beta[classId].getRowDimension(); row++) {
2 26 Feb 07 jari 1063      for (int column = 0; column < beta[classId].getColumnDimension(); column++) {
2 26 Feb 07 jari 1064          System.out.print(beta[classId].get(row, column) + ", ");
2 26 Feb 07 jari 1065      }
2 26 Feb 07 jari 1066      System.out.println(" "); 
2 26 Feb 07 jari 1067        }
2 26 Feb 07 jari 1068           System.out.println(" "); 
2 26 Feb 07 jari 1069         }
2 26 Feb 07 jari 1070 */
2 26 Feb 07 jari 1071
2 26 Feb 07 jari 1072         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1073             probFunction[classId] = ((X_Matrix.transpose()).times(beta[classId])).get(0,0);
2 26 Feb 07 jari 1074         }
2 26 Feb 07 jari 1075
2 26 Feb 07 jari 1076 /*
2 26 Feb 07 jari 1077         System.out.println(" ");
2 26 Feb 07 jari 1078         System.out.println("pdaAlgorithm() -- *******************************");
2 26 Feb 07 jari 1079         System.out.println("pdaAlgorithm() -- probFunction: ");
2 26 Feb 07 jari 1080         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1081             System.out.print(probFunction[classId] + ", ");
2 26 Feb 07 jari 1082         }
2 26 Feb 07 jari 1083         System.out.println(" ");
2 26 Feb 07 jari 1084         System.out.println("pdaAlgorithm() -- *******************************");
2 26 Feb 07 jari 1085 */
2 26 Feb 07 jari 1086 /*
2 26 Feb 07 jari 1087   System.out.println(" ");
2 26 Feb 07 jari 1088   System.out.println("******************** End pdaAlgorithm() ********************");
2 26 Feb 07 jari 1089   System.out.println(" ");
2 26 Feb 07 jari 1090 */
2 26 Feb 07 jari 1091         return probFunction;
2 26 Feb 07 jari 1092
2 26 Feb 07 jari 1093     } // end of pdaAlgorithm
2 26 Feb 07 jari 1094
2 26 Feb 07 jari 1095
2 26 Feb 07 jari 1096
2 26 Feb 07 jari 1097     /**
2 26 Feb 07 jari 1098      * Calculate Parameters for QDA Algorithm:
2 26 Feb 07 jari 1099      *   A(i) = -0.5 * covarianceMatrix(i).inverse()    
2 26 Feb 07 jari 1100      *   C(i) = covarianceMatrix(i).inverse() * meanVector(i) and
2 26 Feb 07 jari 1101      *   c(i) = log(P(Y=i)) - 0.5log(covarianceMatrix(i).det()) 
2 26 Feb 07 jari 1102      *   - 0.5 * meanVector'(i) * covarianceMatrix(i).inverse() * meanVector(i);
2 26 Feb 07 jari 1103      */
2 26 Feb 07 jari 1104     public int calculateQDAParameters(Matrix compMatrix, Matrix responseMatrix) throws AlgorithmException {
2 26 Feb 07 jari 1105
2 26 Feb 07 jari 1106         int numOfGenes = compMatrix.getRowDimension();
2 26 Feb 07 jari 1107         int numOfSamples = compMatrix.getColumnDimension();
2 26 Feb 07 jari 1108
2 26 Feb 07 jari 1109 /*
2 26 Feb 07 jari 1110         System.out.println(" ");
2 26 Feb 07 jari 1111         System.out.println("********** Begin calculateQDAParameters() **********");
2 26 Feb 07 jari 1112         System.out.println(" ");
2 26 Feb 07 jari 1113 */
2 26 Feb 07 jari 1114
2 26 Feb 07 jari 1115   Matrix [] subCompMatrix = new Matrix[numberOfClasses];
2 26 Feb 07 jari 1116
2 26 Feb 07 jari 1117         int [] numOfSamplesInClass = new int[numberOfClasses];
2 26 Feb 07 jari 1118
2 26 Feb 07 jari 1119 /*
2 26 Feb 07 jari 1120         System.out.println("calculateQDAParameters() - compMatrix: " + compMatrix.getRowDimension() + " X " + compMatrix.getColumnDimension() ); 
2 26 Feb 07 jari 1121         printMatrix(compMatrix);
2 26 Feb 07 jari 1122 */
2 26 Feb 07 jari 1123
2 26 Feb 07 jari 1124         for(int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1125
2 26 Feb 07 jari 1126             numOfSamplesInClass[classId] = 0;
2 26 Feb 07 jari 1127             for (int i=0; i<numOfSamples; i++) {
2 26 Feb 07 jari 1128                   numOfSamplesInClass[classId] += responseMatrix.get(i, classId);
2 26 Feb 07 jari 1129             }
2 26 Feb 07 jari 1130
2 26 Feb 07 jari 1131 //            System.out.println("number of samples in class " + classId + " is " + numOfSamplesInClass[classId]);
2 26 Feb 07 jari 1132
2 26 Feb 07 jari 1133             int k = numOfSamplesInClass[classId];
2 26 Feb 07 jari 1134             int [] classIndices = new int[k];
2 26 Feb 07 jari 1135
2 26 Feb 07 jari 1136             int j=0;
2 26 Feb 07 jari 1137             for (int i = 0; i<numOfSamples; i++ ) {
2 26 Feb 07 jari 1138     if (responseMatrix.get(i, classId) == 1) {
2 26 Feb 07 jari 1139                     classIndices[j] = i;
2 26 Feb 07 jari 1140                     j++;
2 26 Feb 07 jari 1141                 }
2 26 Feb 07 jari 1142             }
2 26 Feb 07 jari 1143
2 26 Feb 07 jari 1144             subCompMatrix[classId] = new Matrix(numOfGenes, numOfSamplesInClass[classId]);
2 26 Feb 07 jari 1145       subCompMatrix[classId] = compMatrix.getMatrix(0, numOfGenes-1, classIndices);
2 26 Feb 07 jari 1146         }
2 26 Feb 07 jari 1147
2 26 Feb 07 jari 1148
2 26 Feb 07 jari 1149 /*
2 26 Feb 07 jari 1150         for (int classId=0; classId <numberOfClasses; classId++) { 
2 26 Feb 07 jari 1151             System.out.println("calculateQDAParameters() -- subCompMatrix for class " + classId);
2 26 Feb 07 jari 1152             System.out.println("calculateQDAParameters() - Dimension: " + subCompMatrix[classId].getRowDimension() + " X " + subCompMatrix[classId].getColumnDimension() ); 
2 26 Feb 07 jari 1153             printMatrix(subCompMatrix[classId]);
2 26 Feb 07 jari 1154         }
2 26 Feb 07 jari 1155 */
2 26 Feb 07 jari 1156
2 26 Feb 07 jari 1157         // means for each class
2 26 Feb 07 jari 1158   Matrix [] means = new Matrix[numberOfClasses];
2 26 Feb 07 jari 1159         for (int classId=0; classId < numberOfClasses; classId++) {
2 26 Feb 07 jari 1160              means[classId] = new Matrix(1, numOfGenes);
2 26 Feb 07 jari 1161         }
2 26 Feb 07 jari 1162
2 26 Feb 07 jari 1163         // covariance matrix for each class
2 26 Feb 07 jari 1164   Matrix [] covarianceMatrix = new Matrix[numberOfClasses];
2 26 Feb 07 jari 1165
2 26 Feb 07 jari 1166         for(int classId=0; classId < numberOfClasses; classId++) {
2 26 Feb 07 jari 1167             covarianceMatrix[classId] = new Matrix(numOfGenes, numOfGenes);
2 26 Feb 07 jari 1168         }
2 26 Feb 07 jari 1169
2 26 Feb 07 jari 1170         // inverse of covariance matrix for each class
2 26 Feb 07 jari 1171   Matrix [] covInverseMatrix = new Matrix[numberOfClasses];
2 26 Feb 07 jari 1172         for(int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1173             covInverseMatrix[classId] = new Matrix(numOfGenes, numOfGenes);
2 26 Feb 07 jari 1174         }
2 26 Feb 07 jari 1175
2 26 Feb 07 jari 1176   double [] probability = new double[numberOfClasses]; // probability of a sample in a class
2 26 Feb 07 jari 1177
2 26 Feb 07 jari 1178   double[] sumVector = new double[numOfSamples];
2 26 Feb 07 jari 1179
2 26 Feb 07 jari 1180         A_Matrix = new Matrix[numberOfClasses];
2 26 Feb 07 jari 1181
2 26 Feb 07 jari 1182         C_Matrix = new Matrix[numberOfClasses];
2 26 Feb 07 jari 1183         for (int classId=0; classId < numberOfClasses; classId++) {
2 26 Feb 07 jari 1184             C_Matrix[classId] = new Matrix(numOfGenes,1);
2 26 Feb 07 jari 1185         }
2 26 Feb 07 jari 1186
2 26 Feb 07 jari 1187         cValues = new double[numberOfClasses];
2 26 Feb 07 jari 1188  
2 26 Feb 07 jari 1189         // Calculate probability vector
2 26 Feb 07 jari 1190         int responseSum=0;
2 26 Feb 07 jari 1191         for (int classId=0; classId < numberOfClasses; classId++) {
2 26 Feb 07 jari 1192             responseSum=0;
2 26 Feb 07 jari 1193             for (int sample=0; sample < numOfSamples; sample++) {
2 26 Feb 07 jari 1194                 responseSum += responseMatrix.get(sample, classId);
2 26 Feb 07 jari 1195             }
2 26 Feb 07 jari 1196
2 26 Feb 07 jari 1197             probability[classId] = (double)responseSum/(double)numOfSamples; 
2 26 Feb 07 jari 1198
2 26 Feb 07 jari 1199         }
2 26 Feb 07 jari 1200
2 26 Feb 07 jari 1201         // Calculate gene mean for all the samples for each class 
2 26 Feb 07 jari 1202         double geneSumInOneClass=0.0;
2 26 Feb 07 jari 1203         for (int classId=0; classId < numberOfClasses; classId++) {
2 26 Feb 07 jari 1204             for (int gene=0; gene < numOfGenes; gene++) {
2 26 Feb 07 jari 1205                 geneSumInOneClass=0.0;
2 26 Feb 07 jari 1206                 for (int sample=0; sample < subCompMatrix[classId].getColumnDimension(); sample++) {
2 26 Feb 07 jari 1207                     geneSumInOneClass += subCompMatrix[classId].get(gene, sample);
2 26 Feb 07 jari 1208                 }
2 26 Feb 07 jari 1209                 means[classId].set(0, gene, geneSumInOneClass/subCompMatrix[classId].getColumnDimension()); 
2 26 Feb 07 jari 1210             }
2 26 Feb 07 jari 1211         }
2 26 Feb 07 jari 1212
2 26 Feb 07 jari 1213         // Calculate covariance matrix for the given gene expression matrix that's in one class
2 26 Feb 07 jari 1214
2 26 Feb 07 jari 1215         double [][]tempDouble = new double [numOfGenes][numOfGenes];
2 26 Feb 07 jari 1216         Matrix tempMatrix = new Matrix(numOfGenes, numOfGenes);
2 26 Feb 07 jari 1217
2 26 Feb 07 jari 1218         for (int classId=0; classId < numberOfClasses; classId++) {
2 26 Feb 07 jari 1219
2 26 Feb 07 jari 1220             // calculate covariance matrix for this class
2 26 Feb 07 jari 1221             covarianceMatrix[classId] = getCovarianceMatrix(subCompMatrix[classId], means[classId]);
2 26 Feb 07 jari 1222
2 26 Feb 07 jari 1223 /*
2 26 Feb 07 jari 1224       System.out.println("calculateQDAParameters() -- covariance matrix for class " + classId);
2 26 Feb 07 jari 1225             printMatrix(covarianceMatrix[classId]);
2 26 Feb 07 jari 1226 */
2 26 Feb 07 jari 1227
2 26 Feb 07 jari 1228       double [][] tempArray = covarianceMatrix[classId].getArray();
2 26 Feb 07 jari 1229
2 26 Feb 07 jari 1230       for (int i=0; i<numOfGenes; i++) {
2 26 Feb 07 jari 1231          for (int j=0; j<numOfGenes; j++) {
2 26 Feb 07 jari 1232        tempMatrix.set(i, j, (double) (((int)(tempArray[i][j] * 1E6))/1E6) ); 
2 26 Feb 07 jari 1233          }
2 26 Feb 07 jari 1234       }
2 26 Feb 07 jari 1235
2 26 Feb 07 jari 1236 /*
2 26 Feb 07 jari 1237       System.out.println("calculateQDAParameters() -- temp matrix: ");
2 26 Feb 07 jari 1238       System.out.println("calculateQDAParameters() -- det of covariance matrix = " + tempMatrix.det());
2 26 Feb 07 jari 1239             printMatrix(tempMatrix);
2 26 Feb 07 jari 1240 */
2 26 Feb 07 jari 1241
2 26 Feb 07 jari 1242             if (Math.abs(tempMatrix.det()) < 1E-10) {
2 26 Feb 07 jari 1243
2 26 Feb 07 jari 1244  //               System.out.println("calculateQDAParameters() -- Covariance Matrix is singular for class Id # " + classId);
2 26 Feb 07 jari 1245  //           System.out.println("calculateQDAParameters() -- det of covariance matrix for class " + 
2 26 Feb 07 jari 1246  //                   (covarianceMatrix[classId]).det());
2 26 Feb 07 jari 1247  //               System.out.println(" ");
2 26 Feb 07 jari 1248
2 26 Feb 07 jari 1249                 JOptionPane.showMessageDialog(null, "QDA Algorithm: Covariance Matrix is Singular", "Alert", JOptionPane.WARNING_MESSAGE);
2 26 Feb 07 jari 1250
2 26 Feb 07 jari 1251                 throw new AlgorithmException("QDA Algorithm: Singular Matrix");
2 26 Feb 07 jari 1252             }
2 26 Feb 07 jari 1253
2 26 Feb 07 jari 1254             // calculate inverse of the covariance matrix for this class
2 26 Feb 07 jari 1255 //            covInverseMatrix[classId] = covarianceMatrix[classId].inverse();
2 26 Feb 07 jari 1256             covInverseMatrix[classId] = tempMatrix.inverse();
2 26 Feb 07 jari 1257
2 26 Feb 07 jari 1258 /*
2 26 Feb 07 jari 1259       System.out.println("calculateQDAParameters() -- Inverse of covariance matrix for class " + classId);
2 26 Feb 07 jari 1260        printMatrix(covInverseMatrix[classId]);
2 26 Feb 07 jari 1261 */
2 26 Feb 07 jari 1262
2 26 Feb 07 jari 1263         }
2 26 Feb 07 jari 1264
2 26 Feb 07 jari 1265
2 26 Feb 07 jari 1266 /*
2 26 Feb 07 jari 1267         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1268       System.out.println("calculateQDAParameters() -- mean for class " + classId);
2 26 Feb 07 jari 1269       printMatrix(means[classId]);
2 26 Feb 07 jari 1270         }
2 26 Feb 07 jari 1271 */
2 26 Feb 07 jari 1272
2 26 Feb 07 jari 1273
2 26 Feb 07 jari 1274         // Calculate A matrix for group i
2 26 Feb 07 jari 1275         // A(i) = -0.5 * covarianceMatrix(i).inverse()    
2 26 Feb 07 jari 1276         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1277             A_Matrix[classId] = new Matrix(numOfGenes, numOfGenes);
2 26 Feb 07 jari 1278             A_Matrix[classId] = covInverseMatrix[classId].times((-0.5));
2 26 Feb 07 jari 1279         }
2 26 Feb 07 jari 1280
2 26 Feb 07 jari 1281 /*
2 26 Feb 07 jari 1282         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1283             System.out.println("calculateQDAParameters() -- A_Matrix for class  " + classId);
2 26 Feb 07 jari 1284             printMatrix(A_Matrix[classId]);
2 26 Feb 07 jari 1285         }
2 26 Feb 07 jari 1286 */
2 26 Feb 07 jari 1287
2 26 Feb 07 jari 1288         // Calculate C vector for each class i
2 26 Feb 07 jari 1289         // C(i) = covarianceMatrix(i).inverse() * meanVector(i) and
2 26 Feb 07 jari 1290         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1291             C_Matrix[classId] = covInverseMatrix[classId].times((means[classId]).transpose());
2 26 Feb 07 jari 1292         }
2 26 Feb 07 jari 1293
2 26 Feb 07 jari 1294
2 26 Feb 07 jari 1295 /*
2 26 Feb 07 jari 1296         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1297             System.out.println("calculateQDAParameters() -- C_Matrix for class  " + classId);
2 26 Feb 07 jari 1298       printMatrix(C_Matrix[classId]);
2 26 Feb 07 jari 1299         }
2 26 Feb 07 jari 1300 */
2 26 Feb 07 jari 1301
2 26 Feb 07 jari 1302
2 26 Feb 07 jari 1303         // Calculate c value for group i
2 26 Feb 07 jari 1304         // c(i) = log(P(Y=i)) - 0.5log(covarianceMatrix(i).det()) 
2 26 Feb 07 jari 1305         //        - 0.5 * meanVector'(i) * covarianceMatrix(i).inverse() * meanVector(i);
2 26 Feb 07 jari 1306         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1307     cValues[classId] = Math.log(probability[classId]);
2 26 Feb 07 jari 1308
2 26 Feb 07 jari 1309 //    cValues[classId] -= Math.log(covarianceMatrix[classId].det()) * 0.5;
2 26 Feb 07 jari 1310
2 26 Feb 07 jari 1311     cValues[classId] -= Math.log(tempMatrix.det()) * 0.5;
2 26 Feb 07 jari 1312     cValues[classId] -= (means[classId].times(C_Matrix[classId])).get(0,0) * 0.5;
2 26 Feb 07 jari 1313         }
2 26 Feb 07 jari 1314
2 26 Feb 07 jari 1315
2 26 Feb 07 jari 1316    //     System.out.println(" ");
2 26 Feb 07 jari 1317    //     for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1318    //         System.out.println("calculateQDAParameters() -- cValues for class #" + classId + " is " + cValues[classId]);
2 26 Feb 07 jari 1319    //         System.out.println(" ");
2 26 Feb 07 jari 1320     //    }
2 26 Feb 07 jari 1321    //     System.out.println(" ");
2 26 Feb 07 jari 1322
2 26 Feb 07 jari 1323    //     System.out.println(" ");
2 26 Feb 07 jari 1324    //     System.out.println("********** End calculateQDAParameters() **********");
2 26 Feb 07 jari 1325    //     System.out.println(" ");
2 26 Feb 07 jari 1326
2 26 Feb 07 jari 1327         return 1;
2 26 Feb 07 jari 1328
2 26 Feb 07 jari 1329     } // end of calculateQDAParameters
2 26 Feb 07 jari 1330
2 26 Feb 07 jari 1331
2 26 Feb 07 jari 1332
2 26 Feb 07 jari 1333     /**
2 26 Feb 07 jari 1334      * QDA Algorithm:
2 26 Feb 07 jari 1335      *   The probability of Y=k given X is P(Y=k|X) = exp(Q(k,X))/Sum(exp(Q(i,X))),
2 26 Feb 07 jari 1336      *   where Q(i,X)=X'A(i)X + C'X + c(i) with
2 26 Feb 07 jari 1337      *   A(i) = -0.5 * covarianceMatrix(i).inverse()    
2 26 Feb 07 jari 1338      *   C(i) = covarianceMatrix(i).inverse() * meanVector(i) and
2 26 Feb 07 jari 1339      *   c(i) = log(P(Y=i)) - 0.5log(covarianceMatrix(i).det()) 
2 26 Feb 07 jari 1340      *   - 0.5 * meanVector'(i) * covarianceMatrix(i).inverse() * meanVector(i);
2 26 Feb 07 jari 1341      */
2 26 Feb 07 jari 1342     public double[] qdaAlgorithm(Matrix testMatrix) {
2 26 Feb 07 jari 1343
2 26 Feb 07 jari 1344 //  System.out.println(" ");
2 26 Feb 07 jari 1345 //  System.out.println("************* Begin qdaAlgorithm() **************");
2 26 Feb 07 jari 1346 //  System.out.println(" ");
2 26 Feb 07 jari 1347
2 26 Feb 07 jari 1348         // Calculate Q value for class i
2 26 Feb 07 jari 1349         // Q(i,X)=X'A(i)X + C'X + c(i) 
2 26 Feb 07 jari 1350
2 26 Feb 07 jari 1351   Matrix Q_Matrix = new Matrix(numberOfClasses, 1);
2 26 Feb 07 jari 1352         Matrix temp = new Matrix(1,1);
2 26 Feb 07 jari 1353
2 26 Feb 07 jari 1354 //  System.out.println("qdaAlgorithm() -- testMatrix: ");
2 26 Feb 07 jari 1355 //  printMatrix(testMatrix);
2 26 Feb 07 jari 1356
2 26 Feb 07 jari 1357         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1358
2 26 Feb 07 jari 1359             temp = (testMatrix.transpose()).times(A_Matrix[classId]).times(testMatrix);
2 26 Feb 07 jari 1360             temp = temp.plus((C_Matrix[classId].transpose()).times(testMatrix));
2 26 Feb 07 jari 1361
2 26 Feb 07 jari 1362             Q_Matrix.set(classId, 0, temp.get(0,0) + cValues[classId]);
2 26 Feb 07 jari 1363         }
2 26 Feb 07 jari 1364
2 26 Feb 07 jari 1365   //      System.out.println("qdaAlgorithm() -- Q_Matrix: ");
2 26 Feb 07 jari 1366   //      printMatrix(Q_Matrix);
2 26 Feb 07 jari 1367  
2 26 Feb 07 jari 1368         double sum = 0.0d;
2 26 Feb 07 jari 1369  
2 26 Feb 07 jari 1370         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1371             sum += Math.exp(Q_Matrix.get(classId, 0));
2 26 Feb 07 jari 1372         }
2 26 Feb 07 jari 1373
2 26 Feb 07 jari 1374         // Calculate probability P value for Y=k given X
2 26 Feb 07 jari 1375         // The probability of Y=k given X is P(Y=k|X)=exp(Q(k,X))/Sum(exp(Q(i,X)))
2 26 Feb 07 jari 1376
2 26 Feb 07 jari 1377   double [] probFunction = new double [numberOfClasses];
2 26 Feb 07 jari 1378
2 26 Feb 07 jari 1379         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1380 //            probFunction[classId] = Q_Matrix.get(classId,0);
2 26 Feb 07 jari 1381             probFunction[classId] = Math.exp(Q_Matrix.get(classId,0))/sum;
2 26 Feb 07 jari 1382         }
2 26 Feb 07 jari 1383
2 26 Feb 07 jari 1384  //       System.out.println("qdaAlgorithm() -- probFunction: ");
2 26 Feb 07 jari 1385  //       for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1386  //           System.out.print(probFunction[classId] + ", ");
2 26 Feb 07 jari 1387  //       }
2 26 Feb 07 jari 1388  //       System.out.println(" ");
2 26 Feb 07 jari 1389
2 26 Feb 07 jari 1390 //  System.out.println(" ");
2 26 Feb 07 jari 1391 //  System.out.println("************* End qdaAlgorithm() **************");
2 26 Feb 07 jari 1392 //  System.out.println(" ");
2 26 Feb 07 jari 1393
2 26 Feb 07 jari 1394         return probFunction;
2 26 Feb 07 jari 1395
2 26 Feb 07 jari 1396     } // end of qdaAlgorithm
2 26 Feb 07 jari 1397
2 26 Feb 07 jari 1398
2 26 Feb 07 jari 1399
2 26 Feb 07 jari 1400
2 26 Feb 07 jari 1401     /**
2 26 Feb 07 jari 1402      * MPLS Algorithm:
2 26 Feb 07 jari 1403      *     Set convergence criterion e = 1E-12;
2 26 Feb 07 jari 1404      *     X1 = X 
2 26 Feb 07 jari 1405      *     Y1 = Y 
2 26 Feb 07 jari 1406      * 
2 26 Feb 07 jari 1407      *     for (int k=1; k<Kp; k++) {
2 26 Feb 07 jari 1408      *       u = first column of Yk
2 26 Feb 07 jari 1409      *       delta = 1.0;
2 26 Feb 07 jari 1410      *       while (delta > e) {
2 26 Feb 07 jari 1411      *         w = Xk'u/u'u
2 26 Feb 07 jari 1412      *         t = Xkw
2 26 Feb 07 jari 1413      *         c = Yk't/t't
2 26 Feb 07 jari 1414      *         scale c to unit length
2 26 Feb 07 jari 1415      *         u = Ykc
2 26 Feb 07 jari 1416      *         delta=(w-w_prep)'(w-w_prep), w_prep is previous value of w
2 26 Feb 07 jari 1417      *       }
2 26 Feb 07 jari 1418      * 
2 26 Feb 07 jari 1419      *       ck = c
2 26 Feb 07 jari 1420      *       pk = X't/(t't)
2 26 Feb 07 jari 1421      *       tk = tcp, cp=(pk'pk)^0.5
2 26 Feb 07 jari 1422      *       wk=wcp
2 26 Feb 07 jari 1423      *       bk=u't/(t't)
2 26 Feb 07 jari 1424      *       Xk+1 = Xk - tkpk'
2 26 Feb 07 jari 1425      *       Yk+1 = Yk - bktkck'
2 26 Feb 07 jari 1426      * 
2 26 Feb 07 jari 1427      *       compute test components based on training information from 5
2 26 Feb 07 jari 1428      *       set X*1=X* and compute t*1=X* and compute t*1=X*1W1. Subsequent test components are computed as T*k=X*kWk
2 26 Feb 07 jari 1429      *       where X*k=X*k-1 - t*k-1PK-1, for k=1,2,...kp-1.
2 26 Feb 07 jari 1430      *     }
2 26 Feb 07 jari 1431      *     End of MPLS Alogorithm
2 26 Feb 07 jari 1432      *
2 26 Feb 07 jari 1433      *  Input: subTrainingMatrix cantainning pre-selected genes obtained from geneSelection() algorithm 
2 26 Feb 07 jari 1434      *  Output: geneCompMatrix tMatrix of size Kp x N (where N = numberOfSamples)
2 26 Feb 07 jari 1435      */
2 26 Feb 07 jari 1436     public Matrix mplsAlgorithm(Matrix trainingMatrix, Matrix trainingResponseMatrix, Matrix testMatrix) {
2 26 Feb 07 jari 1437
2 26 Feb 07 jari 1438 //  System.out.println(" ");
2 26 Feb 07 jari 1439 //  System.out.println("************** Begin mplsAlgorithm() **************");
2 26 Feb 07 jari 1440 //  System.out.println(" ");
2 26 Feb 07 jari 1441
2 26 Feb 07 jari 1442   Matrix [] matrix_X = new Matrix[kValue+1]; // training matrix for each k iteration
2 26 Feb 07 jari 1443   Matrix [] matrix_Y = new Matrix[kValue+1]; // response matrix for each k iteration
2 26 Feb 07 jari 1444         Matrix [] matrix_T = new Matrix[kValue]; // normalized test data matrix for each k iteration
2 26 Feb 07 jari 1445
2 26 Feb 07 jari 1446         double tempMean, tempVariance;
2 26 Feb 07 jari 1447
2 26 Feb 07 jari 1448         Matrix training_means_X, training_means_Y;
2 26 Feb 07 jari 1449         Matrix test_means_X;
2 26 Feb 07 jari 1450         Matrix training_variances_X, training_variances_Y;
2 26 Feb 07 jari 1451         Matrix test_variances_X;
2 26 Feb 07 jari 1452         
2 26 Feb 07 jari 1453         int numOfGenes = trainingMatrix.getRowDimension();
2 26 Feb 07 jari 1454         int numOfTrainingSamples = trainingMatrix.getColumnDimension();
2 26 Feb 07 jari 1455         int numOfTestSamples = testMatrix.getColumnDimension();
2 26 Feb 07 jari 1456
2 26 Feb 07 jari 1457 //        System.out.println("mplsAlgorithm() - numOfGenes = " + numOfGenes);
2 26 Feb 07 jari 1458 //        System.out.println("mplsAlgorithm() - numOfTrainingSamples = " + numOfTrainingSamples);
2 26 Feb 07 jari 1459 //        System.out.println("mplsAlgorithm() - numOfTestSamples = " + numOfTestSamples);
2 26 Feb 07 jari 1460 //        System.out.println(" ");
2 26 Feb 07 jari 1461
2 26 Feb 07 jari 1462         AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, kValue, "Dimension Reduction (MPLS)\n");
2 26 Feb 07 jari 1463         fireValueChanged(event);
2 26 Feb 07 jari 1464         event.setId(AlgorithmEvent.PROGRESS_VALUE);            
2 26 Feb 07 jari 1465
2 26 Feb 07 jari 1466         for (int k=0; k < kValue+1; k++) {
2 26 Feb 07 jari 1467            matrix_X[k] = new Matrix(numOfGenes, numOfTrainingSamples); 
2 26 Feb 07 jari 1468      matrix_Y[k] = new Matrix(numOfTrainingSamples, numberOfClasses);
2 26 Feb 07 jari 1469      if (k < kValue) {
2 26 Feb 07 jari 1470               matrix_T[k] = new Matrix(numOfGenes, numOfTestSamples);
2 26 Feb 07 jari 1471            }
2 26 Feb 07 jari 1472         }
2 26 Feb 07 jari 1473  
2 26 Feb 07 jari 1474         
2 26 Feb 07 jari 1475 //        System.out.println("mplsAlgorithm() - Call getSampleMeans for trainingMatrix: ");
2 26 Feb 07 jari 1476         training_means_X = getSampleMeans(trainingMatrix);
2 26 Feb 07 jari 1477         training_variances_X = getSampleVariances(trainingMatrix, training_means_X);
2 26 Feb 07 jari 1478
2 26 Feb 07 jari 1479         // Normalize input training data set (matrix of numOfGenes x numOfTrainingSamples)
2 26 Feb 07 jari 1480         for (int column=0; column<numOfTrainingSamples; column++) {
2 26 Feb 07 jari 1481
2 26 Feb 07 jari 1482            tempMean = training_means_X.get(0, column);
2 26 Feb 07 jari 1483            tempVariance = training_variances_X.get(0, column);
2 26 Feb 07 jari 1484
2 26 Feb 07 jari 1485            for (int row=0; row<numOfGenes; row++) {
2 26 Feb 07 jari 1486               matrix_X[0].set(row, column, (trainingMatrix.get(row,column)-tempMean)/tempVariance); 
2 26 Feb 07 jari 1487            }
2 26 Feb 07 jari 1488         }
2 26 Feb 07 jari 1489
2 26 Feb 07 jari 1490 //        System.out.println("mplsAlgorithm() - Call getSampleMeans for trainingResponseMatrix: ");
2 26 Feb 07 jari 1491         training_means_Y = getSampleMeans(trainingResponseMatrix);
2 26 Feb 07 jari 1492         training_variances_Y = getSampleVariances(trainingResponseMatrix, training_means_Y);
2 26 Feb 07 jari 1493
2 26 Feb 07 jari 1494         // Normalize input response data set
2 26 Feb 07 jari 1495         for (int column=0; column<numberOfClasses; column++) {
2 26 Feb 07 jari 1496
2 26 Feb 07 jari 1497            tempMean = training_means_Y.get(0, column);
2 26 Feb 07 jari 1498            tempVariance = training_variances_Y.get(0, column);
2 26 Feb 07 jari 1499
2 26 Feb 07 jari 1500            for (int row=0; row<numOfTrainingSamples; row++) {
2 26 Feb 07 jari 1501               matrix_Y[0].set(row, column, (trainingResponseMatrix.get(row,column)-tempMean)/tempVariance); 
2 26 Feb 07 jari 1502            }
2 26 Feb 07 jari 1503         }
2 26 Feb 07 jari 1504
2 26 Feb 07 jari 1505         test_means_X = getSampleMeans(testMatrix);
2 26 Feb 07 jari 1506         test_variances_X = getSampleVariances(testMatrix, test_means_X);
2 26 Feb 07 jari 1507
2 26 Feb 07 jari 1508         // Normalize input test data set
2 26 Feb 07 jari 1509         for (int column=0; column<numOfTestSamples; column++) {
2 26 Feb 07 jari 1510
2 26 Feb 07 jari 1511            tempMean = test_means_X.get(0, column);
2 26 Feb 07 jari 1512            tempVariance = test_variances_X.get(0, column);
2 26 Feb 07 jari 1513
2 26 Feb 07 jari 1514            for (int row=0; row<numOfGenes; row++) {
2 26 Feb 07 jari 1515               matrix_T[0].set(row, column, (testMatrix.get(row, column)-tempMean)/tempVariance); 
2 26 Feb 07 jari 1516            }
2 26 Feb 07 jari 1517         }
2 26 Feb 07 jari 1518
2 26 Feb 07 jari 1519         double e = 1E-12;
2 26 Feb 07 jari 1520   Matrix  u = new Matrix(numOfTrainingSamples, 1);  // first column of matrix_Y
2 26 Feb 07 jari 1521         double uu = 0.0; // u'u
2 26 Feb 07 jari 1522
2 26 Feb 07 jari 1523   Matrix  w = new Matrix(numOfGenes, 1);
2 26 Feb 07 jari 1524         Matrix [] wk = new Matrix [kValue];
2 26 Feb 07 jari 1525         for (int i=0; i<kValue; i++) {
2 26 Feb 07 jari 1526           wk[i] = new Matrix(numOfGenes, 1);
2 26 Feb 07 jari 1527         }
2 26 Feb 07 jari 1528
2 26 Feb 07 jari 1529         Matrix w_prep = new Matrix(numOfGenes, 1);
2 26 Feb 07 jari 1530         for (int i=0; i<numOfGenes; i++) {
2 26 Feb 07 jari 1531            w_prep.set(i, 0, 0);
2 26 Feb 07 jari 1532         }
2 26 Feb 07 jari 1533
2 26 Feb 07 jari 1534   Matrix  t = new Matrix(numOfTrainingSamples, 1);
2 26 Feb 07 jari 1535         double tt = 0.0; // t't
2 26 Feb 07 jari 1536         double cc = 0.0; // c'c
2 26 Feb 07 jari 1537         double ww = 0.0; // w'w
2 26 Feb 07 jari 1538         double pp = 0.0; // p'p
2 26 Feb 07 jari 1539
2 26 Feb 07 jari 1540   Matrix  c = new Matrix(numberOfClasses, 1);
2 26 Feb 07 jari 1541
2 26 Feb 07 jari 1542         Matrix [] tk = new Matrix[kValue];
2 26 Feb 07 jari 1543         Matrix [] ck = new Matrix[kValue];
2 26 Feb 07 jari 1544         Matrix [] pk = new Matrix[kValue];
2 26 Feb 07 jari 1545   double [] bk = new double[kValue];
2 26 Feb 07 jari 1546         for (int i=0; i<kValue; i++) {
2 26 Feb 07 jari 1547           tk[i] = new Matrix(numOfTrainingSamples, 1);
2 26 Feb 07 jari 1548           ck[i] = new Matrix(numberOfClasses, 1);
2 26 Feb 07 jari 1549           pk[i] = new Matrix(numOfGenes, 1);
2 26 Feb 07 jari 1550           bk[i] = 0;
2 26 Feb 07 jari 1551         }
2 26 Feb 07 jari 1552
2 26 Feb 07 jari 1553         double cp = 1.0;
2 26 Feb 07 jari 1554        
2 26 Feb 07 jari 1555
2 26 Feb 07 jari 1556         for (int k=0; k < kValue; k++) {
2 26 Feb 07 jari 1557
2 26 Feb 07 jari 1558            // set u = first column of Yk
2 26 Feb 07 jari 1559            for (int row=0; row<numOfTrainingSamples; row++) {
2 26 Feb 07 jari 1560         u.set(row, 0, matrix_Y[k].get(row, 0));
2 26 Feb 07 jari 1561            }
2 26 Feb 07 jari 1562
2 26 Feb 07 jari 1563            // get uu = u'u
2 26 Feb 07 jari 1564            uu = ((u.transpose()).times(u)).get(0,0);
2 26 Feb 07 jari 1565
2 26 Feb 07 jari 1566            int loop = 0;
2 26 Feb 07 jari 1567            double delta = 1.0;
2 26 Feb 07 jari 1568
2 26 Feb 07 jari 1569            while ((delta>e) && loop < 1000) {
2 26 Feb 07 jari 1570  
2 26 Feb 07 jari 1571               loop++;
2 26 Feb 07 jari 1572
2 26 Feb 07 jari 1573               // get w = X[k]'u/u'u (note: here X[k] is genes x samples, so no need to transpose)
2 26 Feb 07 jari 1574               w = matrix_X[k].times(u); // (numOfGenes x 1)
2 26 Feb 07 jari 1575
2 26 Feb 07 jari 1576               for (int row = 0; row < numOfGenes; row++) {
2 26 Feb 07 jari 1577                 w.set(row, 0, w.get(row, 0)/uu);
2 26 Feb 07 jari 1578               }
2 26 Feb 07 jari 1579
2 26 Feb 07 jari 1580               // scale w to unit length
2 26 Feb 07 jari 1581               ww = Math.sqrt(w.transpose().times(w).get(0,0));
2 26 Feb 07 jari 1582               for (int row = 0; row < numOfGenes; row++) {
2 26 Feb 07 jari 1583                 w.set(row, 0, (w.get(row, 0)/ww));
2 26 Feb 07 jari 1584               }
2 26 Feb 07 jari 1585
2 26 Feb 07 jari 1586               // t = X[k]w (numOfSampes x 1)
2 26 Feb 07 jari 1587               t = (matrix_X[k].transpose()).times(w);
2 26 Feb 07 jari 1588
2 26 Feb 07 jari 1589               // c = Yk't/t't (numberOfClasses x 1)
2 26 Feb 07 jari 1590               // scale c to unit length
2 26 Feb 07 jari 1591               tt = ((t.transpose()).times(t)).get(0,0);
2 26 Feb 07 jari 1592               c = (matrix_Y[k].transpose()).times(t);
2 26 Feb 07 jari 1593
2 26 Feb 07 jari 1594               for (int row=0; row<numberOfClasses; row++) {
2 26 Feb 07 jari 1595                   c.set(row, 0, c.get(row, 0)/tt);
2 26 Feb 07 jari 1596               }
2 26 Feb 07 jari 1597
2 26 Feb 07 jari 1598               // scale c to unit length
2 26 Feb 07 jari 1599               cc = Math.sqrt(c.transpose().times(c).get(0,0));
2 26 Feb 07 jari 1600               for (int row = 0; row < numberOfClasses; row++) {
2 26 Feb 07 jari 1601                 c.set(row, 0, (c.get(row, 0)/cc));
2 26 Feb 07 jari 1602               }
2 26 Feb 07 jari 1603
2 26 Feb 07 jari 1604
2 26 Feb 07 jari 1605               // u = Y[k]c
2 26 Feb 07 jari 1606               u = matrix_Y[k].times(c);
2 26 Feb 07 jari 1607
2 26 Feb 07 jari 1608
2 26 Feb 07 jari 1609 /*
2 26 Feb 07 jari 1610           System.out.println(" ");
2 26 Feb 07 jari 1611               if ((loop==6) || (loop == 7) || (loop==8) || (loop==9) || (loop==10) ) {
2 26 Feb 07 jari 1612
2 26 Feb 07 jari 1613       System.out.println(" ");
2 26 Feb 07 jari 1614       System.out.print("mplsAlgorithm() - w-w_prep: ");
2 26 Feb 07 jari 1615       for (int i=0; i<numOfGenes; i++ ) {
2 26 Feb 07 jari 1616           System.out.print((w.minus(w_prep)).get(i,0) + ", ");
2 26 Feb 07 jari 1617       }
2 26 Feb 07 jari 1618       System.out.println(" ");
2 26 Feb 07 jari 1619               System.out.println(" ");
2 26 Feb 07 jari 1620               }
2 26 Feb 07 jari 1621 */
2 26 Feb 07 jari 1622
2 26 Feb 07 jari 1623               // delta=(w-w_prep)'(w-w_prep), w_prep is previous value of w
2 26 Feb 07 jari 1624               if (loop > 1) {
2 26 Feb 07 jari 1625                   delta = (((w.minus(w_prep)).transpose()).times(w.minus(w_prep))).get(0,0);
2 26 Feb 07 jari 1626               }
2 26 Feb 07 jari 1627
2 26 Feb 07 jari 1628  //             System.out.println("mplsAlgorithm() - loop = " + loop + "  delta = " + delta);
2 26 Feb 07 jari 1629
2 26 Feb 07 jari 1630         for(int i=0; i<numOfGenes; i++) {
2 26 Feb 07 jari 1631       w_prep.set(i, 0, w.get(i, 0));
2 26 Feb 07 jari 1632         }
2 26 Feb 07 jari 1633
2 26 Feb 07 jari 1634           }
2 26 Feb 07 jari 1635
2 26 Feb 07 jari 1636           // ck[k] = c
2 26 Feb 07 jari 1637              ck[k] = c;
2 26 Feb 07 jari 1638
2 26 Feb 07 jari 1639           // pk[k] = Xk't/(t't)
2 26 Feb 07 jari 1640           pk[k] = matrix_X[k].times(t);
2 26 Feb 07 jari 1641           for (int row = 0; row < numOfGenes; row++) {
2 26 Feb 07 jari 1642                 pk[k].set(row, 0,  pk[k].get(row, 0)/tt);
2 26 Feb 07 jari 1643           }
2 26 Feb 07 jari 1644
2 26 Feb 07 jari 1645           // cp=(pk[k]'pk[k])^0.5, tk[k] = t*cp
2 26 Feb 07 jari 1646           cp = Math.sqrt(pk[k].transpose().times(pk[k]).get(0,0));
2 26 Feb 07 jari 1647
2 26 Feb 07 jari 1648     // scale pk to unit length
2 26 Feb 07 jari 1649     for (int row = 0; row < numOfGenes; row++) {
2 26 Feb 07 jari 1650       pk[k].set(row, 0, (pk[k].get(row, 0)/cp));
2 26 Feb 07 jari 1651     }
2 26 Feb 07 jari 1652
2 26 Feb 07 jari 1653 //          tk[k] = t.times(cp);
2 26 Feb 07 jari 1654           tk[k] = t;
2 26 Feb 07 jari 1655
2 26 Feb 07 jari 1656           // wk[k]=w*cp
2 26 Feb 07 jari 1657 //          wk[k] = w.times(cp); 
2 26 Feb 07 jari 1658           wk[k] = w; 
2 26 Feb 07 jari 1659
2 26 Feb 07 jari 1660     // bk[k]=u't/(t't)
2 26 Feb 07 jari 1661           bk[k] = (((u.transpose()).times(t)).get(0,0))/tt; 
2 26 Feb 07 jari 1662
2 26 Feb 07 jari 1663           if (k < kValue) {
2 26 Feb 07 jari 1664         // X[k+1] = X[k] - tk[k]*pk[k]'
2 26 Feb 07 jari 1665         matrix_X[k+1] = matrix_X[k].minus((tk[k].times(pk[k].transpose())).transpose());
2 26 Feb 07 jari 1666
2 26 Feb 07 jari 1667         // Y[k+1] = Y[k] - bk[k]*tk[k]*ck'[k]
2 26 Feb 07 jari 1668         matrix_Y[k+1] = matrix_Y[k].minus((tk[k].times(ck[k].transpose())).times(bk[k]));
2 26 Feb 07 jari 1669           }
2 26 Feb 07 jari 1670
2 26 Feb 07 jari 1671           event.setIntValue(k);
2 26 Feb 07 jari 1672           event.setDescription("Calculating Component # " + (k + 1) + "\n");
2 26 Feb 07 jari 1673           fireValueChanged(event); 
2 26 Feb 07 jari 1674        }
2 26 Feb 07 jari 1675
2 26 Feb 07 jari 1676
2 26 Feb 07 jari 1677        // compute test component matrix 
2 26 Feb 07 jari 1678        Matrix [] tMatrix = new Matrix[kValue];
2 26 Feb 07 jari 1679        for (int i=0; i<kValue; i++) {
2 26 Feb 07 jari 1680            tMatrix[i] = new Matrix(numOfTestSamples, 1);
2 26 Feb 07 jari 1681        } 
2 26 Feb 07 jari 1682
2 26 Feb 07 jari 1683        // compute test components based on training information.
2 26 Feb 07 jari 1684        // set X*[1]=X* and compute t*[1]=X* and compute t*1=X*1W1. 
2 26 Feb 07 jari 1685        // Subsequent test components are computed as T*k=X*kWk
2 26 Feb 07 jari 1686        // where X*[k]=X*[k-1] - t*[k-1]P[k-1], for k=1,2,...kp-1.
2 26 Feb 07 jari 1687
2 26 Feb 07 jari 1688        tMatrix[0] = (matrix_T[0].transpose()).times(wk[0]);
2 26 Feb 07 jari 1689
2 26 Feb 07 jari 1690        for (int k=1; k < kValue; k++) {
2 26 Feb 07 jari 1691
2 26 Feb 07 jari 1692            matrix_T[k] = matrix_T[k-1].minus((tMatrix[k-1].times(pk[k-1].transpose())).transpose());
2 26 Feb 07 jari 1693
2 26 Feb 07 jari 1694            tMatrix[k] = (matrix_T[k].transpose()).times(wk[k]); 
2 26 Feb 07 jari 1695
2 26 Feb 07 jari 1696        }
2 26 Feb 07 jari 1697
2 26 Feb 07 jari 1698        Matrix geneCompMatrix = new Matrix(kValue, numOfTestSamples);
2 26 Feb 07 jari 1699
2 26 Feb 07 jari 1700
2 26 Feb 07 jari 1701        for(int gene=0; gene<kValue; gene++) {
2 26 Feb 07 jari 1702            for (int sample=0; sample<numOfTestSamples; sample++) {
2 26 Feb 07 jari 1703                geneCompMatrix.set(gene, sample, tMatrix[gene].get(sample, 0));
2 26 Feb 07 jari 1704            }
2 26 Feb 07 jari 1705        }
2 26 Feb 07 jari 1706
2 26 Feb 07 jari 1707 //  System.out.println(" ");
2 26 Feb 07 jari 1708 //  System.out.println("************** End mplsAlgorithm() **************");
2 26 Feb 07 jari 1709 //  System.out.println(" ");
2 26 Feb 07 jari 1710
2 26 Feb 07 jari 1711        return geneCompMatrix;
2 26 Feb 07 jari 1712
2 26 Feb 07 jari 1713     } // end of mplsAlgorithm
2 26 Feb 07 jari 1714
2 26 Feb 07 jari 1715
2 26 Feb 07 jari 1716
2 26 Feb 07 jari 1717     public Matrix InitialClassification(Matrix expMatrix) throws AlgorithmException {
2 26 Feb 07 jari 1718
2 26 Feb 07 jari 1719 //   System.out.println(" ");
2 26 Feb 07 jari 1720 //     System.out.println("******************** Begin InitialClassification() ********************");
2 26 Feb 07 jari 1721 //   System.out.println(" ");
2 26 Feb 07 jari 1722  
2 26 Feb 07 jari 1723    Matrix selectedTrainingMatrix, trainingRespMatrix, testMatrix;
2 26 Feb 07 jari 1724
2 26 Feb 07 jari 1725    Matrix geneCompMatrix, trainingCompMatrix;
2 26 Feb 07 jari 1726
2 26 Feb 07 jari 1727          Matrix selectedExpMatrix;
2 26 Feb 07 jari 1728
2 26 Feb 07 jari 1729          Matrix compMatrix, compResponseMatrix; 
2 26 Feb 07 jari 1730  
2 26 Feb 07 jari 1731          double[][] probFunction = new double[numberOfSamples][numberOfClasses];
2 26 Feb 07 jari 1732
2 26 Feb 07 jari 1733          int numOfSelectedGenes=1;
2 26 Feb 07 jari 1734          double [] selectedGenesArray = new double[numberOfGenes];
2 26 Feb 07 jari 1735
2 26 Feb 07 jari 1736          int numOfGenes;
2 26 Feb 07 jari 1737          int numOfTrainingSamples = trainingMatrix.getColumnDimension();
2 26 Feb 07 jari 1738          int numOfTestSamples;
2 26 Feb 07 jari 1739          
2 26 Feb 07 jari 1740          AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, numberOfGenes, "Initial Classification\n");
2 26 Feb 07 jari 1741          fireValueChanged(event);
2 26 Feb 07 jari 1742          event.setId(AlgorithmEvent.PROGRESS_VALUE);            
2 26 Feb 07 jari 1743
2 26 Feb 07 jari 1744  
2 26 Feb 07 jari 1745          //  1. (InitialClassification) Select Genes - Select a set, geneSet, of p* genes giving an expression 
2 26 Feb 07 jari 1746          //     matrix expMatrix of size p* x N
2 26 Feb 07 jari 1747
2 26 Feb 07 jari 1748          // Perform gene selection from the training matrix if required
2 26 Feb 07 jari 1749          if (preSelectGenes && (expMatrix.getRowDimension() > kValue)) {
2 26 Feb 07 jari 1750              
2 26 Feb 07 jari 1751              selectedGeneIndices = geneSelection(trainingMatrix);
2 26 Feb 07 jari 1752
2 26 Feb 07 jari 1753              if (selectedGeneIndices == null) return null;
2 26 Feb 07 jari 1754              
2 26 Feb 07 jari 1755              if(highestGeneRank > 1) {
2 26 Feb 07 jari 1756                  numOfSelectedGenes = selectedGeneIndices[highestGeneRank-1].size() + 
2 26 Feb 07 jari 1757                                       selectedGeneIndices[highestGeneRank].size();
2 26 Feb 07 jari 1758              } else if (highestGeneRank == 1){
2 26 Feb 07 jari 1759                  numOfSelectedGenes = selectedGeneIndices[highestGeneRank].size();
2 26 Feb 07 jari 1760              } else if (highestGeneRank == 0) {
2 26 Feb 07 jari 1761                  numOfSelectedGenes = selectedGeneIndices[0].size();
2 26 Feb 07 jari 1762              }
2 26 Feb 07 jari 1763
2 26 Feb 07 jari 1764 //             usedGeneIndices = new int[selectedGeneIndices[highestGeneRank].size()];
2 26 Feb 07 jari 1765 //             unusedGeneIndices = new int[numberOfGenes-selectedGeneIndices[highestGeneRank].size()];
2 26 Feb 07 jari 1766              
2 26 Feb 07 jari 1767 //             System.out.println("***********************************");
2 26 Feb 07 jari 1768 //             System.out.println("selected gene indices size for highest rank= "+selectedGeneIndices[highestGeneRank].size());
2 26 Feb 07 jari 1769 //             System.out.println("highest rank = "+highestGeneRank);
2 26 Feb 07 jari 1770              
2 26 Feb 07 jari 1771              for (int i = 0; i < selectedGeneIndices[highestGeneRank].size(); i++) {
2 26 Feb 07 jari 1772 //                  usedGeneIndices[i] = ((Integer)(selectedGeneIndices[highestGeneRank]).get(i)).intValue();
2 26 Feb 07 jari 1773                  reducedGeneSet[used].add(selectedGeneIndices[highestGeneRank].get(i));  
2 26 Feb 07 jari 1774        }
2 26 Feb 07 jari 1775
2 26 Feb 07 jari 1776              if (highestGeneRank > 1) {
2 26 Feb 07 jari 1777      for (int i = 0; i < selectedGeneIndices[highestGeneRank-1].size(); i++) {
2 26 Feb 07 jari 1778          reducedGeneSet[used].add(selectedGeneIndices[highestGeneRank-1].get(i));  
2 26 Feb 07 jari 1779      }
2 26 Feb 07 jari 1780              }
2 26 Feb 07 jari 1781             
2 26 Feb 07 jari 1782  
2 26 Feb 07 jari 1783              int j=0;
2 26 Feb 07 jari 1784          for (int i = 0; i < numberOfGenes; i++) {
2 26 Feb 07 jari 1785
2 26 Feb 07 jari 1786      if (!isFoundInVector(i, reducedGeneSet[used])) {
2 26 Feb 07 jari 1787                        reducedGeneSet[unused].add(new Integer(i));
2 26 Feb 07 jari 1788                  }
2 26 Feb 07 jari 1789        }
2 26 Feb 07 jari 1790
2 26 Feb 07 jari 1791 /*
2 26 Feb 07 jari 1792              System.out.println("DAM.java: usedGeneIndices size: " + usedGeneIndices.length); 
2 26 Feb 07 jari 1793              for(int i=0; i< usedGeneIndices.length; i++) {
2 26 Feb 07 jari 1794                  System.out.print(usedGeneIndices[i] + ", "); 
2 26 Feb 07 jari 1795              }
2 26 Feb 07 jari 1796              System.out.println(" ");
2 26 Feb 07 jari 1797              System.out.println(" ");
2 26 Feb 07 jari 1798              System.out.println("DAM.java: unusedGeneIndices size: " + unusedGeneIndices.length); 
2 26 Feb 07 jari 1799              for(int i=0; i< unusedGeneIndices.length; i++) {
2 26 Feb 07 jari 1800                  System.out.print(unusedGeneIndices[i] + ", "); 
2 26 Feb 07 jari 1801              }
2 26 Feb 07 jari 1802              System.out.println(" ");
2 26 Feb 07 jari 1803              System.out.println(" ");
2 26 Feb 07 jari 1804 */
2 26 Feb 07 jari 1805
2 26 Feb 07 jari 1806        for (int gene=0; gene<numberOfGenes; gene++) {
2 26 Feb 07 jari 1807      if (isFoundInVector(gene, selectedGeneIndices[highestGeneRank]) || 
2 26 Feb 07 jari 1808                      isFoundInVector(gene, selectedGeneIndices[highestGeneRank-1])  ) {
2 26 Feb 07 jari 1809          selectedGenesArray[gene] = 1;
2 26 Feb 07 jari 1810      } else {
2 26 Feb 07 jari 1811          selectedGenesArray[gene] = 0;
2 26 Feb 07 jari 1812      }
2 26 Feb 07 jari 1813        }
2 26 Feb 07 jari 1814
2 26 Feb 07 jari 1815              //  Construct selectedTrainingMatrix for MPLS algorithm
2 26 Feb 07 jari 1816              selectedTrainingMatrix = new Matrix(numOfSelectedGenes, numOfTrainingSamples);
2 26 Feb 07 jari 1817
2 26 Feb 07 jari 1818              selectedExpMatrix = new Matrix(numOfSelectedGenes, numberOfSamples);
2 26 Feb 07 jari 1819
2 26 Feb 07 jari 1820        int gene=0;
2 26 Feb 07 jari 1821        for (int selected=0; selected<numOfSelectedGenes; selected++) {
2 26 Feb 07 jari 1822
2 26 Feb 07 jari 1823      if (selectedGenesArray[gene] == 1) {
2 26 Feb 07 jari 1824          for (int sample=0; sample<numOfTrainingSamples; sample++) {
2 26 Feb 07 jari 1825        selectedTrainingMatrix.set(selected, sample, trainingMatrix.get(gene, sample));
2 26 Feb 07 jari 1826          }
2 26 Feb 07 jari 1827
2 26 Feb 07 jari 1828                      for (int sample=0; sample<numberOfSamples; sample++) {
2 26 Feb 07 jari 1829        selectedExpMatrix.set(selected, sample, expMatrix.get(gene, sample));
2 26 Feb 07 jari 1830          } 
2 26 Feb 07 jari 1831      }
2 26 Feb 07 jari 1832      gene++;
2 26 Feb 07 jari 1833        }
2 26 Feb 07 jari 1834
2 26 Feb 07 jari 1835 //       System.out.println(" "); 
2 26 Feb 07 jari 1836 //       System.out.println("InitialClassification() - numOfSelectedGenes = " + numOfSelectedGenes); 
2 26 Feb 07 jari 1837 //       System.out.println(" "); 
2 26 Feb 07 jari 1838
2 26 Feb 07 jari 1839          }
2 26 Feb 07 jari 1840          else {
2 26 Feb 07 jari 1841
2 26 Feb 07 jari 1842          for (int i = 0; i < numberOfGenes; i++) {
2 26 Feb 07 jari 1843                  reducedGeneSet[used].add(new Integer(i));  
2 26 Feb 07 jari 1844        }
2 26 Feb 07 jari 1845
2 26 Feb 07 jari 1846              //  Construct selectedTrainingMatrix
2 26 Feb 07 jari 1847              selectedTrainingMatrix = new Matrix(numberOfGenes, numOfTrainingSamples);
2 26 Feb 07 jari 1848        selectedTrainingMatrix = trainingMatrix;
2 26 Feb 07 jari 1849
2 26 Feb 07 jari 1850              selectedExpMatrix = new Matrix(numberOfGenes, numberOfSamples);
2 26 Feb 07 jari 1851              selectedExpMatrix = expMatrix;
2 26 Feb 07 jari 1852          }
2 26 Feb 07 jari 1853
2 26 Feb 07 jari 1854
2 26 Feb 07 jari 1855          // 2. (InitialClassification) Dimension Reduction: Fit PLS (or PCA) to obtain PLS gene components matrix, 
2 26 Feb 07 jari 1856          //    compMatrix, of size NxK
2 26 Feb 07 jari 1857
2 26 Feb 07 jari 1858          // Run MPLS for dimension reduction to obtain components matrix: kValue x N
2 26 Feb 07 jari 1859          trainingRespMatrix = new Matrix(numOfTrainingSamples, numberOfClasses);
2 26 Feb 07 jari 1860    trainingRespMatrix = responseMatrix.getMatrix(trainingIndices, 0, numberOfClasses-1);
2 26 Feb 07 jari 1861
2 26 Feb 07 jari 1862          if (selectedExpMatrix.getRowDimension() == kValue) {
2 26 Feb 07 jari 1863              geneCompMatrix = selectedExpMatrix;
2 26 Feb 07 jari 1864          } else {
2 26 Feb 07 jari 1865              geneCompMatrix = mplsAlgorithm(selectedTrainingMatrix, trainingRespMatrix, selectedExpMatrix);
2 26 Feb 07 jari 1866          }
2 26 Feb 07 jari 1867
2 26 Feb 07 jari 1868          geneComponentMatrix = geneCompMatrix;
2 26 Feb 07 jari 1869
2 26 Feb 07 jari 1870 //         System.out.println(" ");
2 26 Feb 07 jari 1871 //         System.out.println("InitialClassification() - geneCompMatrix row = " + geneCompMatrix.getRowDimension());
2 26 Feb 07 jari 1872 //         System.out.println("InitialClassification() - geneCompMatrix column = " + geneCompMatrix.getColumnDimension());
2 26 Feb 07 jari 1873
2 26 Feb 07 jari 1874 //         System.out.println("InitialClassification() - geneCompMatrix: "); 
2 26 Feb 07 jari 1875  //        printMatrix(geneCompMatrix);
2 26 Feb 07 jari 1876
2 26 Feb 07 jari 1877
2 26 Feb 07 jari 1878          // obtain Gene Comp Matrix and Response Matrix only for the training data
2 26 Feb 07 jari 1879          trainingCompMatrix = geneCompMatrix.getMatrix(0, geneCompMatrix.getRowDimension()-1, trainingIndices);
2 26 Feb 07 jari 1880
2 26 Feb 07 jari 1881 //         System.out.println("InitialClassification() - trainingCompMatrix row = " + trainingCompMatrix.getRowDimension());
2 26 Feb 07 jari 1882 //         System.out.println("InitialClassification() - trainingCompMatrix column = " + trainingCompMatrix.getColumnDimension());
2 26 Feb 07 jari 1883
2 26 Feb 07 jari 1884  //        System.out.println("InitialClassification() - trainingCompMatrix :");
2 26 Feb 07 jari 1885 //         printMatrix(trainingCompMatrix);
2 26 Feb 07 jari 1886
2 26 Feb 07 jari 1887  //        System.out.println("InitialClassification() - trainingRespMatrix :");
2 26 Feb 07 jari 1888  //        printMatrix(trainingRespMatrix);
2 26 Feb 07 jari 1889
2 26 Feb 07 jari 1890          // Call MLE algorithm to calculate beta value from geneComponentMatrix
2 26 Feb 07 jari 1891    if (isPDA == true) {
2 26 Feb 07 jari 1892              try { 
2 26 Feb 07 jari 1893                 beta = mleAlgorithm(trainingCompMatrix, trainingRespMatrix);
2 26 Feb 07 jari 1894              } catch (AlgorithmException ex) {
2 26 Feb 07 jari 1895                 throw new AbortException();
2 26 Feb 07 jari 1896              }
2 26 Feb 07 jari 1897          } else {
2 26 Feb 07 jari 1898              try {
2 26 Feb 07 jari 1899                 calculateQDAParameters(trainingCompMatrix, trainingRespMatrix);
2 26 Feb 07 jari 1900              } catch (AlgorithmException ex) {
2 26 Feb 07 jari 1901                 throw new AbortException();
2 26 Feb 07 jari 1902              }
2 26 Feb 07 jari 1903          }
2 26 Feb 07 jari 1904
2 26 Feb 07 jari 1905
2 26 Feb 07 jari 1906    // 3. (InitialClassification) Classification:  classify all test data
2 26 Feb 07 jari 1907
2 26 Feb 07 jari 1908          // Obtain testMatrix to perform classification 
2 26 Feb 07 jari 1909        numOfGenes = geneCompMatrix.getRowDimension();
2 26 Feb 07 jari 1910        numOfTestSamples = geneCompMatrix.getColumnDimension();
2 26 Feb 07 jari 1911
2 26 Feb 07 jari 1912              int[] columnList = new int[numOfTestSamples-1];
2 26 Feb 07 jari 1913
2 26 Feb 07 jari 1914        testMatrix = new Matrix(numOfGenes, 1);
2 26 Feb 07 jari 1915
2 26 Feb 07 jari 1916 //       System.out.println("InitialClassification() - numOfGenes " + numOfGenes);
2 26 Feb 07 jari 1917 //       System.out.println("InitialClassification() - numOfTestSamples  " + numOfTestSamples);
2 26 Feb 07 jari 1918
2 26 Feb 07 jari 1919              int testSample;
2 26 Feb 07 jari 1920        for(int testSampleId=0; testSampleId<testIndices.length; testSampleId++) {
2 26 Feb 07 jari 1921
2 26 Feb 07 jari 1922                  testSample = testIndices[testSampleId]; 
2 26 Feb 07 jari 1923
2 26 Feb 07 jari 1924  //                System.out.println(" ");
2 26 Feb 07 jari 1925  //                System.out.println("InitialClassification() -- begin testSample = " + testSample);
2 26 Feb 07 jari 1926
2 26 Feb 07 jari 1927      testMatrix = geneCompMatrix.getMatrix(0, numOfGenes-1, testSample, testSample);
2 26 Feb 07 jari 1928
2 26 Feb 07 jari 1929 //     System.out.println(" "); 
2 26 Feb 07 jari 1930 //     System.out.println("InitialClassification() - testMatrix: "); 
2 26 Feb 07 jari 1931 //     printMatrix(testMatrix);
2 26 Feb 07 jari 1932
2 26 Feb 07 jari 1933      if (isPDA == true) {
2 26 Feb 07 jari 1934          probFunction[testSample] = pdaAlgorithm(testMatrix);
2 26 Feb 07 jari 1935
2 26 Feb 07 jari 1936      } else {
2 26 Feb 07 jari 1937          probFunction[testSample] = qdaAlgorithm(testMatrix);
2 26 Feb 07 jari 1938      }
2 26 Feb 07 jari 1939
2 26 Feb 07 jari 1940       } 
2 26 Feb 07 jari 1941
2 26 Feb 07 jari 1942 //      System.out.println(" ");
2 26 Feb 07 jari 1943 //      System.out.println("InitialClassification() - Probability Function: " );
2 26 Feb 07 jari 1944
2 26 Feb 07 jari 1945       for(int testSampleId=0; testSampleId<testIndices.length; testSampleId++) {
2 26 Feb 07 jari 1946
2 26 Feb 07 jari 1947                 testSample = testIndices[testSampleId]; 
2 26 Feb 07 jari 1948
2 26 Feb 07 jari 1949 //          for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1950 //        System.out.print(probFunction[testSample][classId] + ",  ");
2 26 Feb 07 jari 1951  //               }
2 26 Feb 07 jari 1952 //          System.out.println(" ");
2 26 Feb 07 jari 1953       }
2 26 Feb 07 jari 1954 //      System.out.println(" ");
2 26 Feb 07 jari 1955
2 26 Feb 07 jari 1956
2 26 Feb 07 jari 1957             double max=0.0;
2 26 Feb 07 jari 1958             int maxClassId = 0;
2 26 Feb 07 jari 1959
2 26 Feb 07 jari 1960       for(int testSampleId=0; testSampleId<testIndices.length; testSampleId++) {
2 26 Feb 07 jari 1961
2 26 Feb 07 jari 1962                 testSample = testIndices[testSampleId]; 
2 26 Feb 07 jari 1963     maxClassId = 0;
2 26 Feb 07 jari 1964
2 26 Feb 07 jari 1965     if (!Double.isNaN(probFunction[testSample][0])) {
2 26 Feb 07 jari 1966
2 26 Feb 07 jari 1967                     max = probFunction[testSample][0]; 
2 26 Feb 07 jari 1968         maxClassId = 1;
2 26 Feb 07 jari 1969
2 26 Feb 07 jari 1970         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 1971       if (!Double.isNaN(probFunction[testSample][classId])) {
2 26 Feb 07 jari 1972           if (probFunction[testSample][classId] > max) {
2 26 Feb 07 jari 1973              max = probFunction[testSample][classId];
2 26 Feb 07 jari 1974              maxClassId = classId+1;
2 26 Feb 07 jari 1975           }
2 26 Feb 07 jari 1976                         }
2 26 Feb 07 jari 1977                         else {
2 26 Feb 07 jari 1978                             maxClassId = 0;
2 26 Feb 07 jari 1979                             classId = numberOfClasses;
2 26 Feb 07 jari 1980                         }
2 26 Feb 07 jari 1981         }
2 26 Feb 07 jari 1982                 } 
2 26 Feb 07 jari 1983                 else {
2 26 Feb 07 jari 1984                    maxClassId = 0;
2 26 Feb 07 jari 1985                 }
2 26 Feb 07 jari 1986
2 26 Feb 07 jari 1987            classified[maxClassId].add(new Integer(testSample));
2 26 Feb 07 jari 1988
2 26 Feb 07 jari 1989  //     System.out.println("InitialClassification() - Sample # " + testSample + " is in class " + maxClassId);
2 26 Feb 07 jari 1990
2 26 Feb 07 jari 1991       }
2 26 Feb 07 jari 1992 //      System.out.println(" ");
2 26 Feb 07 jari 1993
2 26 Feb 07 jari 1994 //         System.out.println(" "); 
2 26 Feb 07 jari 1995 //         System.out.println("******************** End InitialClassification() ********************"); 
2 26 Feb 07 jari 1996 //         System.out.println(" ");
2 26 Feb 07 jari 1997
2 26 Feb 07 jari 1998          return new Matrix(probFunction);
2 26 Feb 07 jari 1999
2 26 Feb 07 jari 2000     }  // end of initialClassification
2 26 Feb 07 jari 2001
2 26 Feb 07 jari 2002
2 26 Feb 07 jari 2003
2 26 Feb 07 jari 2004     /** 
2 26 Feb 07 jari 2005      * A0Algorithm(): 
2 26 Feb 07 jari 2006      * 1. Select Genes - Select a set, geneSet, of p* genes giving an expression matrix expMatrix of size Nxp*
2 26 Feb 07 jari 2007      * 2. Dimension Reduction: Fit PLS (or PCA) to obtain PLS gene components matrix, compMatrix, of size NxK
2 26 Feb 07 jari 2008      * 3. Classification/Prediction: Classification is based on LOOCV:
2 26 Feb 07 jari 2009      *     For i=1 to N DO
2 26 Feb 07 jari 2010      *       Leave out sample (row) i of compMatrix. Fit classifier to the remaining N-1 
2 26 Feb 07 jari 2011      *       samples and use the fitted classifier to predict left out sample i
2 26 Feb 07 jari 2012      *     End         
2 26 Feb 07 jari 2013      *
2 26 Feb 07 jari 2014      */
2 26 Feb 07 jari 2015     public Matrix A0Algorithm(Matrix expMatrix) throws AlgorithmException{
2 26 Feb 07 jari 2016  
2 26 Feb 07 jari 2017 //   System.out.println(" ");
2 26 Feb 07 jari 2018 //     System.out.println("******************** Begin A0Algorithm() ********************");
2 26 Feb 07 jari 2019 //   System.out.println(" ");
2 26 Feb 07 jari 2020
2 26 Feb 07 jari 2021          Matrix selectedExpMatrix, geneCompMatrix;
2 26 Feb 07 jari 2022
2 26 Feb 07 jari 2023          Matrix subGeneCompMatrix, subResponseMatrix, testMatrix;
2 26 Feb 07 jari 2024  
2 26 Feb 07 jari 2025          double[][] probFunction = new double[numberOfSamples][numberOfClasses];
2 26 Feb 07 jari 2026
2 26 Feb 07 jari 2027          int numOfSelectedGenes=1;
2 26 Feb 07 jari 2028
2 26 Feb 07 jari 2029          double [] selectedGenesArray = new double[numberOfGenes];
2 26 Feb 07 jari 2030
2 26 Feb 07 jari 2031          int numOfGenes;
2 26 Feb 07 jari 2032          int numOfTrainingSamples = trainingMatrix.getColumnDimension();
2 26 Feb 07 jari 2033          int numOfTestSamples;
2 26 Feb 07 jari 2034
2 26 Feb 07 jari 2035          AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, numberOfGenes, "Classification Algorithm : A0\n");
2 26 Feb 07 jari 2036          fireValueChanged(event);
2 26 Feb 07 jari 2037          event.setId(AlgorithmEvent.PROGRESS_VALUE);            
2 26 Feb 07 jari 2038
2 26 Feb 07 jari 2039          
2 26 Feb 07 jari 2040          //  1. (A0Algorithm) Select Genes - Select a set, geneSet, of p* genes giving an expression 
2 26 Feb 07 jari 2041          //     matrix expMatrix of size p* x N
2 26 Feb 07 jari 2042
2 26 Feb 07 jari 2043          // Perform gene selection from the training matrix if required
2 26 Feb 07 jari 2044          if (preSelectGenes && (expMatrix.getRowDimension() > kValue)) {
2 26 Feb 07 jari 2045
2 26 Feb 07 jari 2046              selectedGeneIndices = geneSelection(expMatrix);
2 26 Feb 07 jari 2047
2 26 Feb 07 jari 2048              if (selectedGeneIndices == null) return null;
2 26 Feb 07 jari 2049              
2 26 Feb 07 jari 2050              if(highestGeneRank > 1) {
2 26 Feb 07 jari 2051                  numOfSelectedGenes = selectedGeneIndices[highestGeneRank-1].size() + 
2 26 Feb 07 jari 2052                                       selectedGeneIndices[highestGeneRank].size();
2 26 Feb 07 jari 2053              } else if (highestGeneRank == 1){
2 26 Feb 07 jari 2054                  numOfSelectedGenes = selectedGeneIndices[highestGeneRank].size();
2 26 Feb 07 jari 2055              } else if (highestGeneRank == 0) {
2 26 Feb 07 jari 2056                  numOfSelectedGenes = selectedGeneIndices[0].size();
2 26 Feb 07 jari 2057              }
2 26 Feb 07 jari 2058
2 26 Feb 07 jari 2059          for (int i = 0; i < selectedGeneIndices[highestGeneRank].size(); i++) {
2 26 Feb 07 jari 2060                  reducedGeneSet[used].add(selectedGeneIndices[highestGeneRank].get(i));  
2 26 Feb 07 jari 2061        }
2 26 Feb 07 jari 2062
2 26 Feb 07 jari 2063              if (highestGeneRank > 1) {
2 26 Feb 07 jari 2064      for (int i = 0; i < selectedGeneIndices[highestGeneRank-1].size(); i++) {
2 26 Feb 07 jari 2065          reducedGeneSet[used].add(selectedGeneIndices[highestGeneRank-1].get(i));  
2 26 Feb 07 jari 2066      }
2 26 Feb 07 jari 2067              }
2 26 Feb 07 jari 2068             
2 26 Feb 07 jari 2069              int j=0;
2 26 Feb 07 jari 2070          for (int i = 0; i < numberOfGenes; i++) {
2 26 Feb 07 jari 2071
2 26 Feb 07 jari 2072      if (!isFoundInVector(i, reducedGeneSet[used])) {
2 26 Feb 07 jari 2073                        reducedGeneSet[unused].add(new Integer(i));
2 26 Feb 07 jari 2074                  }
2 26 Feb 07 jari 2075        }
2 26 Feb 07 jari 2076
2 26 Feb 07 jari 2077 /*
2 26 Feb 07 jari 2078
2 26 Feb 07 jari 2079              System.out.println("DAM.java: usedGeneIndices size: " + usedGeneIndices.length); 
2 26 Feb 07 jari 2080              for(int i=0; i< usedGeneIndices.length; i++) {
2 26 Feb 07 jari 2081                  System.out.print(usedGeneIndices[i] + ", "); 
2 26 Feb 07 jari 2082              }
2 26 Feb 07 jari 2083              System.out.println(" ");
2 26 Feb 07 jari 2084              System.out.println(" ");
2 26 Feb 07 jari 2085              System.out.println("DAM.java: unusedGeneIndices size: " + unusedGeneIndices.length); 
2 26 Feb 07 jari 2086              for(int i=0; i< unusedGeneIndices.length; i++) {
2 26 Feb 07 jari 2087                  System.out.print(unusedGeneIndices[i] + ", "); 
2 26 Feb 07 jari 2088              }
2 26 Feb 07 jari 2089              System.out.println(" ");
2 26 Feb 07 jari 2090              System.out.println(" ");
2 26 Feb 07 jari 2091 */
2 26 Feb 07 jari 2092
2 26 Feb 07 jari 2093              // Obtain selectedGenesArray
2 26 Feb 07 jari 2094        for (int gene=0; gene<numberOfGenes; gene++) {
2 26 Feb 07 jari 2095      if (isFoundInVector(gene, selectedGeneIndices[highestGeneRank]) || 
2 26 Feb 07 jari 2096                      isFoundInVector(gene, selectedGeneIndices[highestGeneRank-1])  ) {
2 26 Feb 07 jari 2097          selectedGenesArray[gene] = 1;
2 26 Feb 07 jari 2098      } else {
2 26 Feb 07 jari 2099          selectedGenesArray[gene] = 0;
2 26 Feb 07 jari 2100      }
2 26 Feb 07 jari 2101        }
2 26 Feb 07 jari 2102
2 26 Feb 07 jari 2103              //  Construct selectedExpMatrix for MPLS algorithm
2 26 Feb 07 jari 2104              selectedExpMatrix = new Matrix(numOfSelectedGenes, numberOfSamples);
2 26 Feb 07 jari 2105
2 26 Feb 07 jari 2106        int gene=0;
2 26 Feb 07 jari 2107        for (int selected=0; selected<numOfSelectedGenes; selected++) {
2 26 Feb 07 jari 2108
2 26 Feb 07 jari 2109      if (selectedGenesArray[gene] == 1) {
2 26 Feb 07 jari 2110                      for (int sample=0; sample<numberOfSamples; sample++) {
2 26 Feb 07 jari 2111        selectedExpMatrix.set(selected, sample, expMatrix.get(gene, sample));
2 26 Feb 07 jari 2112          } 
2 26 Feb 07 jari 2113      }
2 26 Feb 07 jari 2114      gene++;
2 26 Feb 07 jari 2115        }
2 26 Feb 07 jari 2116
2 26 Feb 07 jari 2117 //       System.out.println(" "); 
2 26 Feb 07 jari 2118 //       System.out.println("A0Algorithm() - numOfSelectedGenes = " + numOfSelectedGenes); 
2 26 Feb 07 jari 2119 //       System.out.println(" "); 
2 26 Feb 07 jari 2120
2 26 Feb 07 jari 2121          }
2 26 Feb 07 jari 2122          else {
2 26 Feb 07 jari 2123
2 26 Feb 07 jari 2124          for (int i = 0; i < numberOfGenes; i++) {
2 26 Feb 07 jari 2125                  reducedGeneSet[used].add(new Integer(i));  
2 26 Feb 07 jari 2126        }
2 26 Feb 07 jari 2127
2 26 Feb 07 jari 2128              //  Construct selectedExpMatrix
2 26 Feb 07 jari 2129              selectedExpMatrix = new Matrix(numberOfGenes, numberOfSamples);
2 26 Feb 07 jari 2130              selectedExpMatrix = expMatrix;
2 26 Feb 07 jari 2131          }
2 26 Feb 07 jari 2132
2 26 Feb 07 jari 2133
2 26 Feb 07 jari 2134          // 2. (A0Algorithm) Dimension Reduction: Fit PLS (or PCA) to obtain PLS gene components matrix, 
2 26 Feb 07 jari 2135          //    geneCompMatrix, of size NxK
2 26 Feb 07 jari 2136
2 26 Feb 07 jari 2137          // Run MPLS for dimension reduction to obtain components matrix: kValue x N
2 26 Feb 07 jari 2138
2 26 Feb 07 jari 2139          if (selectedExpMatrix.getRowDimension() == kValue) {
2 26 Feb 07 jari 2140              geneCompMatrix = selectedExpMatrix;
2 26 Feb 07 jari 2141          } else {
2 26 Feb 07 jari 2142              geneCompMatrix = mplsAlgorithm(selectedExpMatrix, responseMatrix, selectedExpMatrix);
2 26 Feb 07 jari 2143          }
2 26 Feb 07 jari 2144
2 26 Feb 07 jari 2145          geneComponentMatrix = geneCompMatrix;
2 26 Feb 07 jari 2146
2 26 Feb 07 jari 2147 //         System.out.println(" ");
2 26 Feb 07 jari 2148  //        System.out.println("A0Algorithm() - geneCompMatrix row = " + geneCompMatrix.getRowDimension());
2 26 Feb 07 jari 2149 //         System.out.println("A0Algorithm() - geneCompMatrix column = " + geneCompMatrix.getColumnDimension());
2 26 Feb 07 jari 2150
2 26 Feb 07 jari 2151  //        System.out.println(" "); 
2 26 Feb 07 jari 2152  //        System.out.println("A0Algorithm() - geneCompMatrix: "); 
2 26 Feb 07 jari 2153    //      for (int row = 0; row < geneCompMatrix.getRowDimension(); row++) {
2 26 Feb 07 jari 2154    //          for (int column = 0; column < geneCompMatrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 2155    //              System.out.print(geneCompMatrix.get(row, column) + ", ");
2 26 Feb 07 jari 2156     //         }
2 26 Feb 07 jari 2157   //           System.out.println(" "); 
2 26 Feb 07 jari 2158   //           System.out.println(" "); 
2 26 Feb 07 jari 2159    //      }
2 26 Feb 07 jari 2160
2 26 Feb 07 jari 2161
2 26 Feb 07 jari 2162    // 3. (A0Algorithm) Classification/Prediction: Classification is based on LOOCV:
2 26 Feb 07 jari 2163    //     For i=1 to N DO
2 26 Feb 07 jari 2164    //       Leave out sample (column) i of geneCompMatrix. Fit classifier to the remaining N-1 
2 26 Feb 07 jari 2165    //       samples and use the fitted classifier to predict left out sample i
2 26 Feb 07 jari 2166
2 26 Feb 07 jari 2167          // Obtain testMatrix to perform LOOCV classification on componmentsMatrix
2 26 Feb 07 jari 2168          if (performLOOCV == true) {
2 26 Feb 07 jari 2169        numOfGenes = geneCompMatrix.getRowDimension();
2 26 Feb 07 jari 2170        numOfTestSamples = geneCompMatrix.getColumnDimension();
2 26 Feb 07 jari 2171
2 26 Feb 07 jari 2172              int[] columnList = new int[numOfTestSamples-1];
2 26 Feb 07 jari 2173
2 26 Feb 07 jari 2174        testMatrix = new Matrix(numOfGenes, 1);
2 26 Feb 07 jari 2175        subGeneCompMatrix = new Matrix(numOfGenes, numOfTestSamples-1);
2 26 Feb 07 jari 2176
2 26 Feb 07 jari 2177        subResponseMatrix = new Matrix(numOfTestSamples-1, numberOfClasses);
2 26 Feb 07 jari 2178
2 26 Feb 07 jari 2179 //       System.out.println("A0Algorithm() - numOfGenes " + numOfGenes);
2 26 Feb 07 jari 2180 //       System.out.println("A0Algorithm() - numOfTestSamples  " + numOfTestSamples);
2 26 Feb 07 jari 2181
2 26 Feb 07 jari 2182              //3. Leave out one sample from geneCompMatrix
2 26 Feb 07 jari 2183        for (int leaveOutSample=0; leaveOutSample<numOfTestSamples; leaveOutSample++) {
2 26 Feb 07 jari 2184
2 26 Feb 07 jari 2185 //                 System.out.println("A0Algorithm() -- begin leaveOutSample = " + leaveOutSample);
2 26 Feb 07 jari 2186
2 26 Feb 07 jari 2187                  singularMatrix[leaveOutSample] = false;
2 26 Feb 07 jari 2188
2 26 Feb 07 jari 2189 /*
2 26 Feb 07 jari 2190      System.out.println(" "); 
2 26 Feb 07 jari 2191      System.out.println("A0Algorithm() - geneCompMatrix: "); 
2 26 Feb 07 jari 2192      for (int row = 0; row < geneCompMatrix.getRowDimension(); row++) {
2 26 Feb 07 jari 2193          for (int column = 0; column < geneCompMatrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 2194        System.out.print(geneCompMatrix.get(row, column) + ", ");
2 26 Feb 07 jari 2195          }
2 26 Feb 07 jari 2196          System.out.println(" "); 
2 26 Feb 07 jari 2197          System.out.println(" "); 
2 26 Feb 07 jari 2198      }
2 26 Feb 07 jari 2199 */
2 26 Feb 07 jari 2200
2 26 Feb 07 jari 2201      for (int sample=0; sample<numOfTestSamples; sample++) {
2 26 Feb 07 jari 2202          if (sample < leaveOutSample) {
2 26 Feb 07 jari 2203       columnList[sample] = sample; 
2 26 Feb 07 jari 2204                 } else if ( sample > leaveOutSample ) {
2 26 Feb 07 jari 2205       columnList[sample-1] = sample; 
2 26 Feb 07 jari 2206          }
2 26 Feb 07 jari 2207      }
2 26 Feb 07 jari 2208
2 26 Feb 07 jari 2209 /*     System.out.println(" "); 
2 26 Feb 07 jari 2210      System.out.println("A0Algorithm() - columnList: "); 
2 26 Feb 07 jari 2211      for (int column = 0; column < columnList.length ; column++) {
2 26 Feb 07 jari 2212        System.out.print(columnList[column] + ", ");
2 26 Feb 07 jari 2213      }
2 26 Feb 07 jari 2214      System.out.println(" "); 
2 26 Feb 07 jari 2215      System.out.println(" "); 
2 26 Feb 07 jari 2216 */
2 26 Feb 07 jari 2217      subGeneCompMatrix = geneCompMatrix.getMatrix(0, numOfGenes-1, columnList); 
2 26 Feb 07 jari 2218      subResponseMatrix = responseMatrix.getMatrix(columnList, 0, numberOfClasses-1); 
2 26 Feb 07 jari 2219
2 26 Feb 07 jari 2220 //     System.out.println(" "); 
2 26 Feb 07 jari 2221 //     System.out.println("A0Algorithm() - subGeneCompMatrix: " + subGeneCompMatrix.getRowDimension() + " X " + 
2 26 Feb 07 jari 2222  //                      subGeneCompMatrix.getColumnDimension() ); 
2 26 Feb 07 jari 2223  //                printMatrix(subGeneCompMatrix);
2 26 Feb 07 jari 2224
2 26 Feb 07 jari 2225              testMatrix = geneCompMatrix.getMatrix(0, numOfGenes-1, leaveOutSample, leaveOutSample);
2 26 Feb 07 jari 2226
2 26 Feb 07 jari 2227 //     System.out.println(" "); 
2 26 Feb 07 jari 2228 //     System.out.println("A0Algorithm() - testMatrix: " + testMatrix.getRowDimension() + " X " + testMatrix.getColumnDimension() ); 
2 26 Feb 07 jari 2229 //     System.out.println("A0Algorithm() - testMatrix: "); 
2 26 Feb 07 jari 2230 //                 printMatrix(testMatrix);
2 26 Feb 07 jari 2231
2 26 Feb 07 jari 2232                  // (A0Algorithm) Call PDA or QDA algorithm
2 26 Feb 07 jari 2233      if (isPDA == true) {
2 26 Feb 07 jari 2234                      try { 
2 26 Feb 07 jari 2235                          beta = mleAlgorithm(subGeneCompMatrix, subResponseMatrix);
2 26 Feb 07 jari 2236                      } catch (AlgorithmException ex) {
2 26 Feb 07 jari 2237                          singularMatrix[leaveOutSample] = true;
2 26 Feb 07 jari 2238                          continue;
2 26 Feb 07 jari 2239                      }
2 26 Feb 07 jari 2240          probFunction[leaveOutSample] = pdaAlgorithm(testMatrix);
2 26 Feb 07 jari 2241
2 26 Feb 07 jari 2242      } else {
2 26 Feb 07 jari 2243                      try {
2 26 Feb 07 jari 2244                          calculateQDAParameters(subGeneCompMatrix, subResponseMatrix);
2 26 Feb 07 jari 2245                      } catch (AlgorithmException ex) {
2 26 Feb 07 jari 2246                          singularMatrix[leaveOutSample] = true;
2 26 Feb 07 jari 2247                          continue;
2 26 Feb 07 jari 2248                      }
2 26 Feb 07 jari 2249                      
2 26 Feb 07 jari 2250          probFunction[leaveOutSample] = qdaAlgorithm(testMatrix);
2 26 Feb 07 jari 2251      }
2 26 Feb 07 jari 2252
2 26 Feb 07 jari 2253       } 
2 26 Feb 07 jari 2254
2 26 Feb 07 jari 2255 //      System.out.println(" ");
2 26 Feb 07 jari 2256 //      System.out.println("A0Algorithm() - Probability Function: " );
2 26 Feb 07 jari 2257 //      for (int leaveOutSample=0; leaveOutSample<numOfTestSamples; leaveOutSample++) {
2 26 Feb 07 jari 2258 //          for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 2259 //        System.out.print(probFunction[leaveOutSample][classId] + ",  ");
2 26 Feb 07 jari 2260 //                }
2 26 Feb 07 jari 2261 //          System.out.println(" ");
2 26 Feb 07 jari 2262 //      }
2 26 Feb 07 jari 2263 //      System.out.println(" ");
2 26 Feb 07 jari 2264
2 26 Feb 07 jari 2265
2 26 Feb 07 jari 2266             double max=0.0;
2 26 Feb 07 jari 2267             int maxClassId = 0;
2 26 Feb 07 jari 2268
2 26 Feb 07 jari 2269       for (int leaveOutSample=0; leaveOutSample<numOfTestSamples; leaveOutSample++) {
2 26 Feb 07 jari 2270
2 26 Feb 07 jari 2271                 if (singularMatrix[leaveOutSample] == true) continue;
2 26 Feb 07 jari 2272
2 26 Feb 07 jari 2273                 max=probFunction[leaveOutSample][0];
2 26 Feb 07 jari 2274                 maxClassId = 0;
2 26 Feb 07 jari 2275
2 26 Feb 07 jari 2276     if (!Double.isNaN(probFunction[leaveOutSample][0])) {
2 26 Feb 07 jari 2277
2 26 Feb 07 jari 2278                     max = probFunction[leaveOutSample][0]; 
2 26 Feb 07 jari 2279         maxClassId = 1;
2 26 Feb 07 jari 2280
2 26 Feb 07 jari 2281         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 2282       if (!Double.isNaN(probFunction[leaveOutSample][classId])) {
2 26 Feb 07 jari 2283           if (probFunction[leaveOutSample][classId] > max) {
2 26 Feb 07 jari 2284              max = probFunction[leaveOutSample][classId];
2 26 Feb 07 jari 2285              maxClassId = classId+1;
2 26 Feb 07 jari 2286           }
2 26 Feb 07 jari 2287                         }
2 26 Feb 07 jari 2288                         else {
2 26 Feb 07 jari 2289                             maxClassId = 0;
2 26 Feb 07 jari 2290                             classId = numberOfClasses;
2 26 Feb 07 jari 2291                         }
2 26 Feb 07 jari 2292         }
2 26 Feb 07 jari 2293                 } 
2 26 Feb 07 jari 2294                 else {
2 26 Feb 07 jari 2295                    maxClassId = 0;
2 26 Feb 07 jari 2296                 }
2 26 Feb 07 jari 2297
2 26 Feb 07 jari 2298            classified[maxClassId].add(new Integer(leaveOutSample));
2 26 Feb 07 jari 2299
2 26 Feb 07 jari 2300   //              System.out.println("A0Algorithm() - Sample # " + leaveOutSample + " is in class # " + maxClassId);
2 26 Feb 07 jari 2301
2 26 Feb 07 jari 2302       }
2 26 Feb 07 jari 2303 //      System.out.println(" ");
2 26 Feb 07 jari 2304
2 26 Feb 07 jari 2305          }
2 26 Feb 07 jari 2306
2 26 Feb 07 jari 2307    //      System.out.println("DAM.java - A0Algorithm() -  classified.length = " + classified.length);
2 26 Feb 07 jari 2308
2 26 Feb 07 jari 2309   // for (int i = 0; i < classified.length; i++) {
2 26 Feb 07 jari 2310    //          System.out.println("DAM.java - A0Algorithm() -  classified[" + i + "].size() = " + classified[i].size());
2 26 Feb 07 jari 2311   // }
2 26 Feb 07 jari 2312
2 26 Feb 07 jari 2313   // System.out.println(" ");
2 26 Feb 07 jari 2314
2 26 Feb 07 jari 2315   // System.out.println(" ");
2 26 Feb 07 jari 2316      //System.out.println("******************** End A0Algorithm() ********************");
2 26 Feb 07 jari 2317    //System.out.println(" ");
2 26 Feb 07 jari 2318          
2 26 Feb 07 jari 2319          return new Matrix(probFunction);
2 26 Feb 07 jari 2320     } // end of A0Algorithm
2 26 Feb 07 jari 2321
2 26 Feb 07 jari 2322
2 26 Feb 07 jari 2323
2 26 Feb 07 jari 2324     /**
2 26 Feb 07 jari 2325      * A1Algorithm(): 
2 26 Feb 07 jari 2326      * 1. Select Genes: Select a set, geneSet, of p* genes giving an expression matrix, expMatrix of size Nxp*
2 26 Feb 07 jari 2327      *    For i=1 to N Do
2 26 Feb 07 jari 2328      *        Leave out sample (row) i of expression matrix expMatrix, say subTrainingMatrix
2 26 Feb 07 jari 2329      * 2. Dimension Reduction: Fit PLS using subTrainingMatrix to obtain PLS gene component matrix, geneCompMatrix
2 26 Feb 07 jari 2330      * 3. Classification/Prediction: Fit classifier to the remaining N-1 samples. 
2 26 Feb 07 jari 2331      *    i.e. using geneCompMatrix. Use the fitted classifier to predict left out sample i.
2 26 Feb 07 jari 2332      */
2 26 Feb 07 jari 2333     public Matrix A1Algorithm(Matrix expMatrix) throws AlgorithmException{
2 26 Feb 07 jari 2334
2 26 Feb 07 jari 2335   // System.out.println(" ");
2 26 Feb 07 jari 2336   //   System.out.println("******************** Begin A1Algorithm() ********************");
2 26 Feb 07 jari 2337   // System.out.println(" ");
2 26 Feb 07 jari 2338
2 26 Feb 07 jari 2339          Matrix selectedExpMatrix;
2 26 Feb 07 jari 2340
2 26 Feb 07 jari 2341          Matrix subSelectedExpMatrix, subResponseMatrix;
2 26 Feb 07 jari 2342
2 26 Feb 07 jari 2343          Matrix geneCompMatrix, subGeneCompMatrix, testMatrix;
2 26 Feb 07 jari 2344
2 26 Feb 07 jari 2345          double[][] probFunction = new double[numberOfSamples][numberOfClasses];
2 26 Feb 07 jari 2346
2 26 Feb 07 jari 2347          double [] selectedGenesArray = new double[numberOfGenes];
2 26 Feb 07 jari 2348
2 26 Feb 07 jari 2349          int numOfSelectedGenes=1;
2 26 Feb 07 jari 2350
2 26 Feb 07 jari 2351          AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, numberOfGenes, "Classification Algorithm : A1\n");
2 26 Feb 07 jari 2352          fireValueChanged(event);
2 26 Feb 07 jari 2353          event.setId(AlgorithmEvent.PROGRESS_VALUE);            
2 26 Feb 07 jari 2354
2 26 Feb 07 jari 2355          //1. (A1Algorithm) Select Genes: Select a set, geneSet, of p* genes giving a 
2 26 Feb 07 jari 2356          // expression matrix of size p* X N
2 26 Feb 07 jari 2357
2 26 Feb 07 jari 2358          // Perform gene selection from the training matrix if required
2 26 Feb 07 jari 2359          if (preSelectGenes && (expMatrix.getRowDimension() > kValue)) {
2 26 Feb 07 jari 2360
2 26 Feb 07 jari 2361              selectedGeneIndices = geneSelection(expMatrix);
2 26 Feb 07 jari 2362
2 26 Feb 07 jari 2363              if (selectedGeneIndices == null) return null;
2 26 Feb 07 jari 2364              
2 26 Feb 07 jari 2365              if(highestGeneRank > 1) {
2 26 Feb 07 jari 2366                  numOfSelectedGenes = selectedGeneIndices[highestGeneRank-1].size() + 
2 26 Feb 07 jari 2367                                       selectedGeneIndices[highestGeneRank].size();
2 26 Feb 07 jari 2368              } else if (highestGeneRank == 1){
2 26 Feb 07 jari 2369                  numOfSelectedGenes = selectedGeneIndices[highestGeneRank].size();
2 26 Feb 07 jari 2370              } else if (highestGeneRank == 0) {
2 26 Feb 07 jari 2371                  numOfSelectedGenes = selectedGeneIndices[0].size();
2 26 Feb 07 jari 2372              }
2 26 Feb 07 jari 2373
2 26 Feb 07 jari 2374          for (int i = 0; i < selectedGeneIndices[highestGeneRank].size(); i++) {
2 26 Feb 07 jari 2375                  reducedGeneSet[used].add(selectedGeneIndices[highestGeneRank].get(i));  
2 26 Feb 07 jari 2376        }
2 26 Feb 07 jari 2377
2 26 Feb 07 jari 2378              if (highestGeneRank > 1) {
2 26 Feb 07 jari 2379      for (int i = 0; i < selectedGeneIndices[highestGeneRank-1].size(); i++) {
2 26 Feb 07 jari 2380          reducedGeneSet[used].add(selectedGeneIndices[highestGeneRank-1].get(i));  
2 26 Feb 07 jari 2381      }
2 26 Feb 07 jari 2382              }
2 26 Feb 07 jari 2383             
2 26 Feb 07 jari 2384  
2 26 Feb 07 jari 2385              int j=0;
2 26 Feb 07 jari 2386          for (int i = 0; i < numberOfGenes; i++) {
2 26 Feb 07 jari 2387
2 26 Feb 07 jari 2388      if (!isFoundInVector(i, reducedGeneSet[used])) {
2 26 Feb 07 jari 2389                        reducedGeneSet[unused].add(new Integer(i));
2 26 Feb 07 jari 2390                  }
2 26 Feb 07 jari 2391        }
2 26 Feb 07 jari 2392
2 26 Feb 07 jari 2393
2 26 Feb 07 jari 2394        for (int gene=0; gene<numberOfGenes; gene++) {
2 26 Feb 07 jari 2395      if (isFoundInVector(gene, selectedGeneIndices[highestGeneRank]) || 
2 26 Feb 07 jari 2396                      isFoundInVector(gene, selectedGeneIndices[highestGeneRank-1])  ) {
2 26 Feb 07 jari 2397          selectedGenesArray[gene] = 1;
2 26 Feb 07 jari 2398      } else {
2 26 Feb 07 jari 2399          selectedGenesArray[gene] = 0;
2 26 Feb 07 jari 2400      }
2 26 Feb 07 jari 2401        }
2 26 Feb 07 jari 2402
2 26 Feb 07 jari 2403              //  Construct selectedExpMatrix for MPLS algorithm
2 26 Feb 07 jari 2404              selectedExpMatrix = new Matrix(numOfSelectedGenes, numberOfSamples);
2 26 Feb 07 jari 2405
2 26 Feb 07 jari 2406        int gene=0;
2 26 Feb 07 jari 2407        for (int selected=0; selected<numOfSelectedGenes; selected++) {
2 26 Feb 07 jari 2408
2 26 Feb 07 jari 2409      if (selectedGenesArray[gene] == 1) {
2 26 Feb 07 jari 2410
2 26 Feb 07 jari 2411                      for (int sample=0; sample<numberOfSamples; sample++) {
2 26 Feb 07 jari 2412        selectedExpMatrix.set(selected, sample, expMatrix.get(gene, sample));
2 26 Feb 07 jari 2413          } 
2 26 Feb 07 jari 2414      }
2 26 Feb 07 jari 2415      gene++;
2 26 Feb 07 jari 2416        }
2 26 Feb 07 jari 2417
2 26 Feb 07 jari 2418    //    System.out.println(" "); 
2 26 Feb 07 jari 2419    //    System.out.println("A1Algorithm() - numOfSelectedGenes = " + numOfSelectedGenes); 
2 26 Feb 07 jari 2420    //    System.out.println(" "); 
2 26 Feb 07 jari 2421
2 26 Feb 07 jari 2422          } else {
2 26 Feb 07 jari 2423
2 26 Feb 07 jari 2424          for (int i = 0; i < numberOfGenes; i++) {
2 26 Feb 07 jari 2425                  reducedGeneSet[used].add(new Integer(i));  
2 26 Feb 07 jari 2426        }
2 26 Feb 07 jari 2427
2 26 Feb 07 jari 2428              //  Construct selectedExpMatrix
2 26 Feb 07 jari 2429              selectedExpMatrix = new Matrix(numberOfGenes, numberOfSamples);
2 26 Feb 07 jari 2430              selectedExpMatrix = expMatrix;
2 26 Feb 07 jari 2431          }
2 26 Feb 07 jari 2432
2 26 Feb 07 jari 2433
2 26 Feb 07 jari 2434          // 2&3. (A1Algorithm) For i=1 to N Do
2 26 Feb 07 jari 2435          //    Leave out sample (row) i of selectedExpMatrix
2 26 Feb 07 jari 2436
2 26 Feb 07 jari 2437          for (int leaveOutSample = 0; leaveOutSample < numberOfSamples; leaveOutSample++) {
2 26 Feb 07 jari 2438
2 26 Feb 07 jari 2439              singularMatrix[leaveOutSample] = false;
2 26 Feb 07 jari 2440
2 26 Feb 07 jari 2441        int gene=0;
2 26 Feb 07 jari 2442
2 26 Feb 07 jari 2443         //     System.out.println("A1Algorithm() -- begin leaveOutSample = " + leaveOutSample);
2 26 Feb 07 jari 2444
2 26 Feb 07 jari 2445        int[] columnList = new int[numberOfSamples-1];
2 26 Feb 07 jari 2446
2 26 Feb 07 jari 2447        for (int sample=0; sample<numberOfSamples; sample++) {
2 26 Feb 07 jari 2448      if (sample < leaveOutSample) {
2 26 Feb 07 jari 2449         columnList[sample] = sample; 
2 26 Feb 07 jari 2450      } else if (sample > leaveOutSample) {
2 26 Feb 07 jari 2451         columnList[sample-1] = sample; 
2 26 Feb 07 jari 2452      }
2 26 Feb 07 jari 2453        }
2 26 Feb 07 jari 2454
2 26 Feb 07 jari 2455   //     System.out.println(" "); 
2 26 Feb 07 jari 2456   //     System.out.println("A1Algorithm() - columnList: "); 
2 26 Feb 07 jari 2457   //     for (int column = 0; column < columnList.length ; column++) {
2 26 Feb 07 jari 2458   //       System.out.print(columnList[column] + ", ");
2 26 Feb 07 jari 2459     //         }
2 26 Feb 07 jari 2460   //     System.out.println(" "); 
2 26 Feb 07 jari 2461    //    System.out.println(" "); 
2 26 Feb 07 jari 2462
2 26 Feb 07 jari 2463
2 26 Feb 07 jari 2464              //  2. (A1Algorithm) Dimension Reduction: Fit PLS using selectedExpMatrix to obtain PLS gene component matrix, geneCompMatrix
2 26 Feb 07 jari 2465
2 26 Feb 07 jari 2466              // obtain subSelectedExpMatrix and subResponseMatrix for MPLS algorithm
2 26 Feb 07 jari 2467           subSelectedExpMatrix = selectedExpMatrix.getMatrix(0, numOfSelectedGenes-1, columnList); 
2 26 Feb 07 jari 2468          subResponseMatrix = responseMatrix.getMatrix(columnList, 0, numberOfClasses-1); 
2 26 Feb 07 jari 2469
2 26 Feb 07 jari 2470        // Run MPLS for dimension reduction to obtain components matrix: kValue x N
2 26 Feb 07 jari 2471        if (selectedExpMatrix.getRowDimension() == kValue) {
2 26 Feb 07 jari 2472      geneCompMatrix = selectedExpMatrix;
2 26 Feb 07 jari 2473        } else {
2 26 Feb 07 jari 2474      geneCompMatrix = mplsAlgorithm(subSelectedExpMatrix, subResponseMatrix, selectedExpMatrix);
2 26 Feb 07 jari 2475        }
2 26 Feb 07 jari 2476
2 26 Feb 07 jari 2477 //      geneCompMatrix = mplsAlgorithm(subSelectedExpMatrix, subResponseMatrix, selectedExpMatrix);
2 26 Feb 07 jari 2478
2 26 Feb 07 jari 2479              geneComponentMatrix = geneCompMatrix;
2 26 Feb 07 jari 2480
2 26 Feb 07 jari 2481 /*
2 26 Feb 07 jari 2482        System.out.println(" "); 
2 26 Feb 07 jari 2483        System.out.println("A1Algorithm() - geneCompMatrix obtained from mplsAlgorithm(): "); 
2 26 Feb 07 jari 2484        for (int row = 0; row < geneCompMatrix.getRowDimension(); row++) {
2 26 Feb 07 jari 2485      for (int column = 0; column < geneCompMatrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 2486          System.out.print(geneCompMatrix.get(row, column) + ", ");
2 26 Feb 07 jari 2487      }
2 26 Feb 07 jari 2488      System.out.println(" "); 
2 26 Feb 07 jari 2489      System.out.println(" "); 
2 26 Feb 07 jari 2490        }
2 26 Feb 07 jari 2491 */
2 26 Feb 07 jari 2492
2 26 Feb 07 jari 2493
2 26 Feb 07 jari 2494        //  3. Classification/Prediction: Fit classifier to the remaining N-1 samples. 
2 26 Feb 07 jari 2495              //     i.e. using geneCompMatrix. Use the fitted classifier to predict left out sample i.
2 26 Feb 07 jari 2496
2 26 Feb 07 jari 2497        // Obtain testMatrix 
2 26 Feb 07 jari 2498  
2 26 Feb 07 jari 2499              int numOfGenes = geneCompMatrix.getRowDimension();
2 26 Feb 07 jari 2500
2 26 Feb 07 jari 2501        testMatrix = new Matrix(numOfGenes, 1);
2 26 Feb 07 jari 2502        testMatrix = geneCompMatrix.getMatrix(0, numOfGenes-1, leaveOutSample, leaveOutSample);
2 26 Feb 07 jari 2503
2 26 Feb 07 jari 2504        subGeneCompMatrix = new Matrix(numOfGenes, numberOfSamples-1);
2 26 Feb 07 jari 2505        subGeneCompMatrix = geneCompMatrix.getMatrix(0, numOfGenes-1, columnList); 
2 26 Feb 07 jari 2506
2 26 Feb 07 jari 2507 /*
2 26 Feb 07 jari 2508        System.out.println(" "); 
2 26 Feb 07 jari 2509        System.out.println("A1Algorithm() - testMatrix: "); 
2 26 Feb 07 jari 2510        for (int row = 0; row < testMatrix.getRowDimension(); row++) {
2 26 Feb 07 jari 2511      for (int column = 0; column < testMatrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 2512          System.out.print(testMatrix.get(row, column) + ", ");
2 26 Feb 07 jari 2513      }
2 26 Feb 07 jari 2514      System.out.println(" "); 
2 26 Feb 07 jari 2515      System.out.println(" "); 
2 26 Feb 07 jari 2516        }
2 26 Feb 07 jari 2517
2 26 Feb 07 jari 2518        System.out.println(" "); 
2 26 Feb 07 jari 2519        System.out.println("A1Algorithm() - subGeneCompMatrix: "); 
2 26 Feb 07 jari 2520        for (int row = 0; row < subGeneCompMatrix.getRowDimension(); row++) {
2 26 Feb 07 jari 2521      for (int column = 0; column < subGeneCompMatrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 2522          System.out.print(subGeneCompMatrix.get(row, column) + ", ");
2 26 Feb 07 jari 2523      }
2 26 Feb 07 jari 2524      System.out.println(" "); 
2 26 Feb 07 jari 2525      System.out.println(" "); 
2 26 Feb 07 jari 2526        }
2 26 Feb 07 jari 2527 */
2 26 Feb 07 jari 2528
2 26 Feb 07 jari 2529              // (A1Algorithm) Call PDA or QDA algorithm
2 26 Feb 07 jari 2530        if (isPDA == true) {
2 26 Feb 07 jari 2531                  try { 
2 26 Feb 07 jari 2532                      beta = mleAlgorithm(subGeneCompMatrix, subResponseMatrix);
2 26 Feb 07 jari 2533                  } catch (AlgorithmException ex) {
2 26 Feb 07 jari 2534                      singularMatrix[leaveOutSample] = true;
2 26 Feb 07 jari 2535                      continue;
2 26 Feb 07 jari 2536                  }
2 26 Feb 07 jari 2537              probFunction[leaveOutSample] = pdaAlgorithm(testMatrix);
2 26 Feb 07 jari 2538
2 26 Feb 07 jari 2539        } else {
2 26 Feb 07 jari 2540                  try {
2 26 Feb 07 jari 2541                      calculateQDAParameters(subGeneCompMatrix, subResponseMatrix);
2 26 Feb 07 jari 2542                  } catch (AlgorithmException ex) {
2 26 Feb 07 jari 2543                      singularMatrix[leaveOutSample] = true;
2 26 Feb 07 jari 2544                      continue;
2 26 Feb 07 jari 2545                  }
2 26 Feb 07 jari 2546            probFunction[leaveOutSample] = qdaAlgorithm(testMatrix);
2 26 Feb 07 jari 2547        }
2 26 Feb 07 jari 2548          } 
2 26 Feb 07 jari 2549
2 26 Feb 07 jari 2550    /* System.out.println(" ");
2 26 Feb 07 jari 2551   // System.out.println("A1Algorithm() - Probablity Function: " );
2 26 Feb 07 jari 2552    for (int leaveOutSample=0; leaveOutSample<numberOfSamples; leaveOutSample++) {
2 26 Feb 07 jari 2553       for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 2554     System.out.print(probFunction[leaveOutSample][classId] + ",  ");
2 26 Feb 07 jari 2555       }
2 26 Feb 07 jari 2556       System.out.println(" ");
2 26 Feb 07 jari 2557    }
2 26 Feb 07 jari 2558    System.out.println(" ");
2 26 Feb 07 jari 2559 */
2 26 Feb 07 jari 2560
2 26 Feb 07 jari 2561             double max=0.0;
2 26 Feb 07 jari 2562             int maxClassId = 0;
2 26 Feb 07 jari 2563
2 26 Feb 07 jari 2564       for (int leaveOutSample=0; leaveOutSample<numberOfSamples; leaveOutSample++) {
2 26 Feb 07 jari 2565
2 26 Feb 07 jari 2566                 if (singularMatrix[leaveOutSample] == true) continue;
2 26 Feb 07 jari 2567
2 26 Feb 07 jari 2568                 max=probFunction[leaveOutSample][0];
2 26 Feb 07 jari 2569                 maxClassId = 0;
2 26 Feb 07 jari 2570
2 26 Feb 07 jari 2571     if (!Double.isNaN(probFunction[leaveOutSample][0])) {
2 26 Feb 07 jari 2572
2 26 Feb 07 jari 2573                     max = probFunction[leaveOutSample][0]; 
2 26 Feb 07 jari 2574         maxClassId = 1;
2 26 Feb 07 jari 2575
2 26 Feb 07 jari 2576         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 2577       if (!Double.isNaN(probFunction[leaveOutSample][classId])) {
2 26 Feb 07 jari 2578           if (probFunction[leaveOutSample][classId] > max) {
2 26 Feb 07 jari 2579              max = probFunction[leaveOutSample][classId];
2 26 Feb 07 jari 2580              maxClassId = classId+1;
2 26 Feb 07 jari 2581           }
2 26 Feb 07 jari 2582                         }
2 26 Feb 07 jari 2583                         else {
2 26 Feb 07 jari 2584                             maxClassId = 0;
2 26 Feb 07 jari 2585                             classId = numberOfClasses;
2 26 Feb 07 jari 2586                         }
2 26 Feb 07 jari 2587         }
2 26 Feb 07 jari 2588                 } 
2 26 Feb 07 jari 2589                 else {
2 26 Feb 07 jari 2590                    maxClassId = 0;
2 26 Feb 07 jari 2591                 }
2 26 Feb 07 jari 2592
2 26 Feb 07 jari 2593            classified[maxClassId].add(new Integer(leaveOutSample));
2 26 Feb 07 jari 2594
2 26 Feb 07 jari 2595            //     System.out.println("A1Algorithm() - Sample # " + leaveOutSample + " is in class # " + maxClassId);
2 26 Feb 07 jari 2596
2 26 Feb 07 jari 2597       }
2 26 Feb 07 jari 2598     //  System.out.println(" ");
2 26 Feb 07 jari 2599
2 26 Feb 07 jari 2600       //   System.out.println("DAM.java - A1Algorithm() -  classified.length = " + classified.length);
2 26 Feb 07 jari 2601
2 26 Feb 07 jari 2602   // for (int i = 0; i < classified.length; i++) {
2 26 Feb 07 jari 2603     //         System.out.println("DAM.java - A1Algorithm() -  classified[" + i + "].size() = " + classified[i].size());
2 26 Feb 07 jari 2604   // }
2 26 Feb 07 jari 2605
2 26 Feb 07 jari 2606   // System.out.println(" ");
2 26 Feb 07 jari 2607
2 26 Feb 07 jari 2608   // System.out.println(" ");
2 26 Feb 07 jari 2609     // System.out.println("******************** End A1Algorithm() ********************");
2 26 Feb 07 jari 2610   // System.out.println(" ");
2 26 Feb 07 jari 2611
2 26 Feb 07 jari 2612          return new Matrix(probFunction);
2 26 Feb 07 jari 2613     } // end of A1Algorithm
2 26 Feb 07 jari 2614
2 26 Feb 07 jari 2615
2 26 Feb 07 jari 2616
2 26 Feb 07 jari 2617     /*
2 26 Feb 07 jari 2618      *  A2Algorithm(): 
2 26 Feb 07 jari 2619      *  For i=1 to N Do
2 26 Feb 07 jari 2620      *     Leave out sample (row) i of expression matrix expMatrix to obtain subExpMatrix
2 26 Feb 07 jari 2621      *  1. Select Genes: Select a set, geneSet, of p* genes giving an expression matrix, subTrainingMatrix of size N-1xp*
2 26 Feb 07 jari 2622      *  2. Dimension Reduction: Fit PLS using subTrainingMatrix to obtain PLS gene component matrix, geneCompMatrix
2 26 Feb 07 jari 2623      *  3. Classification/Prediction: Fit classifier to the remaining N-1 samples. 
2 26 Feb 07 jari 2624      *     i.e. using geneCompMatrix. Use the fitted classifier to predict left out sample i.
2 26 Feb 07 jari 2625      */
2 26 Feb 07 jari 2626     public Matrix A2Algorithm(Matrix expMatrix) throws AlgorithmException{
2 26 Feb 07 jari 2627
2 26 Feb 07 jari 2628   // System.out.println(" ");
2 26 Feb 07 jari 2629     // System.out.println("******************** Begin A2Algorithm() ********************");
2 26 Feb 07 jari 2630   // System.out.println(" ");
2 26 Feb 07 jari 2631
2 26 Feb 07 jari 2632          Matrix subExpMatrix, subResponseMatrix;
2 26 Feb 07 jari 2633
2 26 Feb 07 jari 2634          Matrix selectedExpMatrix, subSelectedExpMatrix;
2 26 Feb 07 jari 2635
2 26 Feb 07 jari 2636          Matrix geneCompMatrix, subGeneCompMatrix, testMatrix;
2 26 Feb 07 jari 2637
2 26 Feb 07 jari 2638          double[][] probFunction = new double[numberOfSamples][numberOfClasses];
2 26 Feb 07 jari 2639
2 26 Feb 07 jari 2640          double [] selectedGenesArray = new double[numberOfGenes];
2 26 Feb 07 jari 2641
2 26 Feb 07 jari 2642          int numOfSelectedGenes=1;
2 26 Feb 07 jari 2643
2 26 Feb 07 jari 2644          AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, numberOfGenes, "Classification Algorithm : A2\n");
2 26 Feb 07 jari 2645          fireValueChanged(event);
2 26 Feb 07 jari 2646          event.setId(AlgorithmEvent.PROGRESS_VALUE);            
2 26 Feb 07 jari 2647
2 26 Feb 07 jari 2648    // For i=1 to N Do
2 26 Feb 07 jari 2649    //    Leave out sample (row) i of expMatrix to obtain a subExpMatrix
2 26 Feb 07 jari 2650
2 26 Feb 07 jari 2651         // System.out.println("A2Algorithm() Start ---------");
2 26 Feb 07 jari 2652
2 26 Feb 07 jari 2653
2 26 Feb 07 jari 2654          for (int leaveOutSample=0; leaveOutSample<numberOfSamples; leaveOutSample++) {
2 26 Feb 07 jari 2655
2 26 Feb 07 jari 2656              singularMatrix[leaveOutSample] = false;
2 26 Feb 07 jari 2657
2 26 Feb 07 jari 2658         //     System.out.println("A2Algorithm() -- For leaveOutSamp = " + leaveOutSample + " Start: ");
2 26 Feb 07 jari 2659
2 26 Feb 07 jari 2660        // 1. Select Genes: Select a set, geneSet, of p* genes giving  expression matrix, 
2 26 Feb 07 jari 2661              //    TrainingsubMatrix of size (N-1) x p*
2 26 Feb 07 jari 2662
2 26 Feb 07 jari 2663              int[] columnList = new int[numberOfSamples-1];
2 26 Feb 07 jari 2664
2 26 Feb 07 jari 2665              for (int sample=0; sample<numberOfSamples; sample++) { 
2 26 Feb 07 jari 2666      if (sample < leaveOutSample) {
2 26 Feb 07 jari 2667         columnList[sample] = sample; 
2 26 Feb 07 jari 2668             } else if (sample > leaveOutSample) {
2 26 Feb 07 jari 2669         columnList[sample-1] = sample; 
2 26 Feb 07 jari 2670      }
2 26 Feb 07 jari 2671        }
2 26 Feb 07 jari 2672
2 26 Feb 07 jari 2673   /*     System.out.println(" "); 
2 26 Feb 07 jari 2674        System.out.println("A2Algorithm() - columnList: "); 
2 26 Feb 07 jari 2675        for (int column = 0; column < columnList.length ; column++) {
2 26 Feb 07 jari 2676          System.out.print(columnList[column] + ", ");
2 26 Feb 07 jari 2677              }
2 26 Feb 07 jari 2678        System.out.println(" "); 
2 26 Feb 07 jari 2679        System.out.println(" "); 
2 26 Feb 07 jari 2680 */
2 26 Feb 07 jari 2681        subExpMatrix = expMatrix.getMatrix(0, numberOfGenes-1, columnList); 
2 26 Feb 07 jari 2682          subResponseMatrix = responseMatrix.getMatrix(columnList, 0, numberOfClasses-1); 
2 26 Feb 07 jari 2683          
2 26 Feb 07 jari 2684        // Perform gene selection if required
2 26 Feb 07 jari 2685              if (preSelectGenes && (expMatrix.getRowDimension() > kValue)) {
2 26 Feb 07 jari 2686
2 26 Feb 07 jari 2687      selectedGeneIndices = geneSelection(subExpMatrix);
2 26 Feb 07 jari 2688
2 26 Feb 07 jari 2689      if (selectedGeneIndices == null) return null;
2 26 Feb 07 jari 2690
2 26 Feb 07 jari 2691      if(highestGeneRank > 1) {
2 26 Feb 07 jari 2692          numOfSelectedGenes = selectedGeneIndices[highestGeneRank-1].size() + 
2 26 Feb 07 jari 2693             selectedGeneIndices[highestGeneRank].size();
2 26 Feb 07 jari 2694      } else if (highestGeneRank == 1){
2 26 Feb 07 jari 2695          numOfSelectedGenes = selectedGeneIndices[highestGeneRank].size();
2 26 Feb 07 jari 2696      } else if (highestGeneRank == 0) {
2 26 Feb 07 jari 2697          numOfSelectedGenes = selectedGeneIndices[0].size();
2 26 Feb 07 jari 2698      }
2 26 Feb 07 jari 2699
2 26 Feb 07 jari 2700
2 26 Feb 07 jari 2701      for (int i = 0; i < selectedGeneIndices[highestGeneRank].size(); i++) {
2 26 Feb 07 jari 2702          reducedGeneSetForA2[leaveOutSample][used].add(selectedGeneIndices[highestGeneRank].get(i));  
2 26 Feb 07 jari 2703      }
2 26 Feb 07 jari 2704
2 26 Feb 07 jari 2705      if (highestGeneRank > 1) {
2 26 Feb 07 jari 2706          for (int i = 0; i < selectedGeneIndices[highestGeneRank-1].size(); i++) {
2 26 Feb 07 jari 2707        reducedGeneSetForA2[leaveOutSample][used].add(selectedGeneIndices[highestGeneRank-1].get(i));  
2 26 Feb 07 jari 2708          }
2 26 Feb 07 jari 2709      }
2 26 Feb 07 jari 2710
2 26 Feb 07 jari 2711
2 26 Feb 07 jari 2712      int j=0;
2 26 Feb 07 jari 2713      for (int i = 0; i < numberOfGenes; i++) {
2 26 Feb 07 jari 2714
2 26 Feb 07 jari 2715          if (!isFoundInVector(i, reducedGeneSetForA2[leaveOutSample][used])) {
2 26 Feb 07 jari 2716          reducedGeneSetForA2[leaveOutSample][unused].add(new Integer(i));
2 26 Feb 07 jari 2717          }
2 26 Feb 07 jari 2718      }
2 26 Feb 07 jari 2719
2 26 Feb 07 jari 2720
2 26 Feb 07 jari 2721      for (int gene=0; gene<numberOfGenes; gene++) {
2 26 Feb 07 jari 2722          if ( (isFoundInVector(gene, selectedGeneIndices[highestGeneRank])) || 
2 26 Feb 07 jari 2723                           ((highestGeneRank > 0) && isFoundInVector(gene, selectedGeneIndices[highestGeneRank-1])) )  
2 26 Feb 07 jari 2724                      {
2 26 Feb 07 jari 2725        selectedGenesArray[gene] = 1;
2 26 Feb 07 jari 2726          } else {
2 26 Feb 07 jari 2727        selectedGenesArray[gene] = 0;
2 26 Feb 07 jari 2728          }
2 26 Feb 07 jari 2729      }
2 26 Feb 07 jari 2730
2 26 Feb 07 jari 2731   //   System.out.println(" "); 
2 26 Feb 07 jari 2732   //   System.out.println("A2Algorithm() - numOfSelectedGenes = " + numOfSelectedGenes); 
2 26 Feb 07 jari 2733   //   System.out.println(" "); 
2 26 Feb 07 jari 2734
2 26 Feb 07 jari 2735      //  Construct selectedExpMatrix for MPLS algorithm
2 26 Feb 07 jari 2736      selectedExpMatrix = new Matrix(numOfSelectedGenes, numberOfSamples);
2 26 Feb 07 jari 2737
2 26 Feb 07 jari 2738      int gene=0;
2 26 Feb 07 jari 2739      for (int selected=0; selected<numOfSelectedGenes; selected++) {
2 26 Feb 07 jari 2740
2 26 Feb 07 jari 2741          if (selectedGenesArray[gene] == 1) {
2 26 Feb 07 jari 2742        for (int sample=0; sample<numberOfSamples; sample++) {
2 26 Feb 07 jari 2743            selectedExpMatrix.set(selected, sample, expMatrix.get(gene, sample));
2 26 Feb 07 jari 2744        }
2 26 Feb 07 jari 2745          }
2 26 Feb 07 jari 2746          gene++;
2 26 Feb 07 jari 2747      }
2 26 Feb 07 jari 2748
2 26 Feb 07 jari 2749   //   System.out.println(" "); 
2 26 Feb 07 jari 2750   //   System.out.println("A2Algorithm() - numOfSelectedGenes = " + numOfSelectedGenes); 
2 26 Feb 07 jari 2751   //   System.out.println(" "); 
2 26 Feb 07 jari 2752                  
2 26 Feb 07 jari 2753          } else {
2 26 Feb 07 jari 2754
2 26 Feb 07 jari 2755                for (int i = 0; i < numberOfGenes; i++) {
2 26 Feb 07 jari 2756                      reducedGeneSetForA2[leaveOutSample][used].add(new Integer(i));  
2 26 Feb 07 jari 2757            }
2 26 Feb 07 jari 2758
2 26 Feb 07 jari 2759      //  Construct selectedExpMatrix
2 26 Feb 07 jari 2760      selectedExpMatrix = new Matrix(numberOfGenes, numberOfSamples);
2 26 Feb 07 jari 2761      selectedExpMatrix = expMatrix;
2 26 Feb 07 jari 2762        }
2 26 Feb 07 jari 2763
2 26 Feb 07 jari 2764        //  2. (A2Algorithm) Dimension Reduction: Fit PLS using subExpMatrix to obtain PLS gene component matrix, geneCompMatrix
2 26 Feb 07 jari 2765
2 26 Feb 07 jari 2766        //  Construct subSelectedExpMatrix for MPLS algorithm
2 26 Feb 07 jari 2767        subSelectedExpMatrix = new Matrix(numOfSelectedGenes, numberOfSamples-1);
2 26 Feb 07 jari 2768        subSelectedExpMatrix = selectedExpMatrix.getMatrix(0, numOfSelectedGenes-1, columnList);
2 26 Feb 07 jari 2769
2 26 Feb 07 jari 2770        // Run MPLS for dimension reduction to obtain components matrix: N x kValue
2 26 Feb 07 jari 2771        if (selectedExpMatrix.getRowDimension() == kValue) {
2 26 Feb 07 jari 2772      geneCompMatrix = selectedExpMatrix;
2 26 Feb 07 jari 2773        } else {
2 26 Feb 07 jari 2774                  geneCompMatrix = mplsAlgorithm(subSelectedExpMatrix, subResponseMatrix, selectedExpMatrix);
2 26 Feb 07 jari 2775        }
2 26 Feb 07 jari 2776 //             geneCompMatrix = mplsAlgorithm(subSelectedExpMatrix, subResponseMatrix, selectedExpMatrix);
2 26 Feb 07 jari 2777
2 26 Feb 07 jari 2778              geneComponentMatrix = geneCompMatrix;
2 26 Feb 07 jari 2779
2 26 Feb 07 jari 2780 /*
2 26 Feb 07 jari 2781        System.out.println(" "); 
2 26 Feb 07 jari 2782        System.out.println("A2Algorithm() - geneCompMatrix: "); 
2 26 Feb 07 jari 2783        for (int row = 0; row < geneCompMatrix.getRowDimension(); row++) {
2 26 Feb 07 jari 2784      for (int column = 0; column < geneCompMatrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 2785          System.out.print(geneCompMatrix.get(row, column) + ", ");
2 26 Feb 07 jari 2786      }
2 26 Feb 07 jari 2787      System.out.println(" "); 
2 26 Feb 07 jari 2788      System.out.println(" "); 
2 26 Feb 07 jari 2789        }
2 26 Feb 07 jari 2790 */
2 26 Feb 07 jari 2791
2 26 Feb 07 jari 2792
2 26 Feb 07 jari 2793        //  3. Classification/Prediction: Fit classifier to the remaining N-1 samples. 
2 26 Feb 07 jari 2794          //     i.e. using geneCompMatrix. Use the fitted classifier to predict left out sample i.
2 26 Feb 07 jari 2795
2 26 Feb 07 jari 2796        // Obtain testMatrix
2 26 Feb 07 jari 2797              int numOfGenes = geneCompMatrix.getRowDimension();
2 26 Feb 07 jari 2798        testMatrix = new Matrix(numOfGenes, 1);
2 26 Feb 07 jari 2799        testMatrix = geneCompMatrix.getMatrix(0, numOfGenes-1, leaveOutSample, leaveOutSample);
2 26 Feb 07 jari 2800
2 26 Feb 07 jari 2801           subGeneCompMatrix = new Matrix(numOfGenes, numberOfSamples-1);
2 26 Feb 07 jari 2802          subGeneCompMatrix = geneCompMatrix.getMatrix(0, numOfGenes-1, columnList); 
2 26 Feb 07 jari 2803
2 26 Feb 07 jari 2804
2 26 Feb 07 jari 2805 /*
2 26 Feb 07 jari 2806        System.out.println(" "); 
2 26 Feb 07 jari 2807        System.out.println("A2Algorithm() - testMatrix: "); 
2 26 Feb 07 jari 2808        for (int row = 0; row < testMatrix.getRowDimension(); row++) {
2 26 Feb 07 jari 2809      for (int column = 0; column < testMatrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 2810          System.out.print(testMatrix.get(row, column) + ", ");
2 26 Feb 07 jari 2811      }
2 26 Feb 07 jari 2812      System.out.println(" "); 
2 26 Feb 07 jari 2813      System.out.println(" "); 
2 26 Feb 07 jari 2814        }
2 26 Feb 07 jari 2815
2 26 Feb 07 jari 2816        System.out.println(" "); 
2 26 Feb 07 jari 2817        System.out.println("A2Algorithm() - subGeneCompMatrix: "); 
2 26 Feb 07 jari 2818        for (int row = 0; row < subGeneCompMatrix.getRowDimension(); row++) {
2 26 Feb 07 jari 2819      for (int column = 0; column < subGeneCompMatrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 2820          System.out.print(subGeneCompMatrix.get(row, column) + ", ");
2 26 Feb 07 jari 2821      }
2 26 Feb 07 jari 2822      System.out.println(" "); 
2 26 Feb 07 jari 2823      System.out.println(" "); 
2 26 Feb 07 jari 2824        }
2 26 Feb 07 jari 2825 */
2 26 Feb 07 jari 2826
2 26 Feb 07 jari 2827
2 26 Feb 07 jari 2828              // (A2Algorithm) Call PDA or QDA algorithm
2 26 Feb 07 jari 2829        if (isPDA == true) {
2 26 Feb 07 jari 2830                  try { 
2 26 Feb 07 jari 2831                      beta = mleAlgorithm(subGeneCompMatrix, subResponseMatrix);
2 26 Feb 07 jari 2832                  } catch (AlgorithmException ex) {
2 26 Feb 07 jari 2833                      singularMatrix[leaveOutSample] = true;
2 26 Feb 07 jari 2834                      continue;
2 26 Feb 07 jari 2835                  }
2 26 Feb 07 jari 2836              probFunction[leaveOutSample] = pdaAlgorithm(testMatrix);
2 26 Feb 07 jari 2837
2 26 Feb 07 jari 2838        } else {
2 26 Feb 07 jari 2839                  try {
2 26 Feb 07 jari 2840                      calculateQDAParameters(subGeneCompMatrix, subResponseMatrix);
2 26 Feb 07 jari 2841                  } catch (AlgorithmException ex) {
2 26 Feb 07 jari 2842                      singularMatrix[leaveOutSample] = true;
2 26 Feb 07 jari 2843                      continue;
2 26 Feb 07 jari 2844                  }
2 26 Feb 07 jari 2845            probFunction[leaveOutSample] = qdaAlgorithm(testMatrix);
2 26 Feb 07 jari 2846        }
2 26 Feb 07 jari 2847
2 26 Feb 07 jari 2848      //        System.out.println("A2Algorithm() -- For leaveOutSamp = " + leaveOutSample + " End: ");
2 26 Feb 07 jari 2849       //       System.out.println(" ");
2 26 Feb 07 jari 2850       //       System.out.println(" ");
2 26 Feb 07 jari 2851
2 26 Feb 07 jari 2852          } 
2 26 Feb 07 jari 2853
2 26 Feb 07 jari 2854   /* System.out.println(" ");
2 26 Feb 07 jari 2855   // System.out.println("A2Algorithm() - Probablity Function: " );
2 26 Feb 07 jari 2856    for (int leaveOutSample=0; leaveOutSample<numberOfSamples; leaveOutSample++) {
2 26 Feb 07 jari 2857       for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 2858     System.out.print(probFunction[leaveOutSample][classId] + ",  ");
2 26 Feb 07 jari 2859       }
2 26 Feb 07 jari 2860       System.out.println(" ");
2 26 Feb 07 jari 2861    }
2 26 Feb 07 jari 2862    System.out.println(" ");
2 26 Feb 07 jari 2863 */
2 26 Feb 07 jari 2864
2 26 Feb 07 jari 2865             double max=0.0;
2 26 Feb 07 jari 2866             int maxClassId = 0;
2 26 Feb 07 jari 2867
2 26 Feb 07 jari 2868       for (int leaveOutSample=0; leaveOutSample<numberOfSamples; leaveOutSample++) {
2 26 Feb 07 jari 2869
2 26 Feb 07 jari 2870                 if (singularMatrix[leaveOutSample] == true) continue;
2 26 Feb 07 jari 2871
2 26 Feb 07 jari 2872                 max=probFunction[leaveOutSample][0];
2 26 Feb 07 jari 2873                 maxClassId = 0;
2 26 Feb 07 jari 2874
2 26 Feb 07 jari 2875     if (!Double.isNaN(probFunction[leaveOutSample][0])) {
2 26 Feb 07 jari 2876
2 26 Feb 07 jari 2877                     max = probFunction[leaveOutSample][0]; 
2 26 Feb 07 jari 2878         maxClassId = 1;
2 26 Feb 07 jari 2879
2 26 Feb 07 jari 2880         for (int classId=0; classId<numberOfClasses; classId++) {
2 26 Feb 07 jari 2881       if (!Double.isNaN(probFunction[leaveOutSample][classId])) {
2 26 Feb 07 jari 2882           if (probFunction[leaveOutSample][classId] > max) {
2 26 Feb 07 jari 2883              max = probFunction[leaveOutSample][classId];
2 26 Feb 07 jari 2884              maxClassId = classId+1;
2 26 Feb 07 jari 2885           }
2 26 Feb 07 jari 2886                         }
2 26 Feb 07 jari 2887                         else {
2 26 Feb 07 jari 2888                             maxClassId = 0;
2 26 Feb 07 jari 2889                             classId = numberOfClasses;
2 26 Feb 07 jari 2890                         }
2 26 Feb 07 jari 2891         }
2 26 Feb 07 jari 2892                 } 
2 26 Feb 07 jari 2893                 else {
2 26 Feb 07 jari 2894                    maxClassId = 0;
2 26 Feb 07 jari 2895                 }
2 26 Feb 07 jari 2896
2 26 Feb 07 jari 2897            classified[maxClassId].add(new Integer(leaveOutSample));
2 26 Feb 07 jari 2898
2 26 Feb 07 jari 2899           //      System.out.println("A2Algorithm() - Sample # " + leaveOutSample + " is in class # " + maxClassId);
2 26 Feb 07 jari 2900
2 26 Feb 07 jari 2901       }
2 26 Feb 07 jari 2902   //    System.out.println(" ");
2 26 Feb 07 jari 2903
2 26 Feb 07 jari 2904
2 26 Feb 07 jari 2905    //      System.out.println("DAM.java - A2Algorithm() -  classified.length = " + classified.length);
2 26 Feb 07 jari 2906
2 26 Feb 07 jari 2907   // for (int i = 0; i < classified.length; i++) {
2 26 Feb 07 jari 2908   //           System.out.println("DAM.java - A2Algorithm() -  classified[" + i + "].size() = " + classified[i].size());
2 26 Feb 07 jari 2909   // }
2 26 Feb 07 jari 2910
2 26 Feb 07 jari 2911   // System.out.println(" ");
2 26 Feb 07 jari 2912
2 26 Feb 07 jari 2913   // System.out.println(" ");
2 26 Feb 07 jari 2914     // System.out.println("******************** End A2Algorithm() ********************");
2 26 Feb 07 jari 2915   // System.out.println(" ");
2 26 Feb 07 jari 2916
2 26 Feb 07 jari 2917          return new Matrix(probFunction);
2 26 Feb 07 jari 2918     } // end of A2algorithm
2 26 Feb 07 jari 2919
2 26 Feb 07 jari 2920
2 26 Feb 07 jari 2921     /**
2 26 Feb 07 jari 2922      * getTValue(): obtain the t value from T table
2 26 Feb 07 jari 2923      *
2 26 Feb 07 jari 2924      */
2 26 Feb 07 jari 2925     public double getTValue(int degreeOfFreedom, double alpha) {
2 26 Feb 07 jari 2926
2 26 Feb 07 jari 2927       double ttbl[][] = {
2 26 Feb 07 jari 2928
2 26 Feb 07 jari 2929     { 1 , 6.314,   12.706,  31.821,  63.657, 636.619},
2 26 Feb 07 jari 2930     { 2 , 2.920 ,  4.303 ,  6.965 ,  9.925 , 31.598 },
2 26 Feb 07 jari 2931     { 3 , 2.353 ,  3.182 ,  4.541 ,  5.841 , 12.941 },
2 26 Feb 07 jari 2932     { 4 , 2.132 ,  2.776 ,  3.747 ,  4.604 ,  8.610 },
2 26 Feb 07 jari 2933     { 5 , 2.015 ,  2.571 ,  3.365 ,  4.032 ,  6.859 },
2 26 Feb 07 jari 2934     { 6 , 1.943 ,  2.447 ,  3.143 ,  3.707 ,  5.959 },
2 26 Feb 07 jari 2935     { 7 , 1.895 ,  2.365 ,  2.998 ,  3.499 ,  5.405 },
2 26 Feb 07 jari 2936     { 8 , 1.860 ,  2.306 ,  2.896 ,  3.355 ,  5.041 },
2 26 Feb 07 jari 2937     { 9 , 1.833 ,  2.262 ,  2.821 ,  3.250 ,  4.781 },
2 26 Feb 07 jari 2938     {10 , 1.812 ,  2.228 ,  2.764 ,  3.169 ,  4.587 },
2 26 Feb 07 jari 2939     {11 , 1.796 ,  2.201 ,  2.718 ,  3.106 ,  4.437 },
2 26 Feb 07 jari 2940     {12 , 1.782 ,  2.179 ,  2.681 ,  3.055 ,  4.318 },
2 26 Feb 07 jari 2941     {13 , 1.771 ,  2.160 ,  2.650 ,  3.012 ,  4.221 },
2 26 Feb 07 jari 2942     {14 , 1.761 ,  2.145 ,  2.624 ,  2.997 ,  4.140 },
2 26 Feb 07 jari 2943     {15 , 1.753 ,  2.131 ,  2.602 ,  2.947 ,  4.073 },
2 26 Feb 07 jari 2944     {16 , 1.746 ,  2.120 ,  2.583 ,  2.921 ,  4.015 },
2 26 Feb 07 jari 2945     {17 , 1.740 ,  2.110 ,  2.567 ,  2.898 ,  3.965 },
2 26 Feb 07 jari 2946     {18 , 1.734 ,  2.101 ,  2.552 ,  2.878 ,  3.922 },
2 26 Feb 07 jari 2947     {19 , 1.729 ,  2.093 ,  2.539 ,  2.861 ,  3.883 },
2 26 Feb 07 jari 2948     {20 , 1.725 ,  2.086 ,  2.528 ,  2.845 ,  3.850 },
2 26 Feb 07 jari 2949     {21 , 1.721 ,  2.080 ,  2.518 ,  2.831 ,  3.819 },
2 26 Feb 07 jari 2950     {22 , 1.717 ,  2.074 ,  2.508 ,  2.819 ,  3.792 },
2 26 Feb 07 jari 2951     {23 , 1.714 ,  2.069 ,  2.500 ,  2.807 ,  3.767 },
2 26 Feb 07 jari 2952     {24 , 1.711 ,  2.064 ,  2.492 ,  2.797 ,  3.745 },
2 26 Feb 07 jari 2953     {25 , 1.708 ,  2.060 ,  2.485 ,  2.787 ,  3.725 },
2 26 Feb 07 jari 2954     {26 , 1.706 ,  2.056 ,  2.479 ,  2.779 ,  3.707 },
2 26 Feb 07 jari 2955     {27 , 1.703 ,  2.052 ,  2.473 ,  2.771 ,  3.690 },
2 26 Feb 07 jari 2956     {28 , 1.701 ,  2.048 ,  2.467 ,  2.763 ,  3.674 },
2 26 Feb 07 jari 2957     {29 , 1.699 ,  2.045 ,  2.462 ,  2.756 ,  3.659 },
2 26 Feb 07 jari 2958     {30 , 1.697 ,  2.042 ,  2.457 ,  2.750 ,  3.646 },
2 26 Feb 07 jari 2959     {40 , 1.684 ,  2.021 ,  2.423 ,  2.704 ,  3.551 },
2 26 Feb 07 jari 2960     {60 , 1.671 ,  2.000 ,  2.390 ,  2.660 ,  3.460 },
2 26 Feb 07 jari 2961     {120 , 1.658 , 1.980 ,  2.358 ,  2.617 ,  3.373 },
2 26 Feb 07 jari 2962     {1000 , 1.645 , 1.960 , 2.326 ,  2.576 ,  3.291 } 
2 26 Feb 07 jari 2963       };
2 26 Feb 07 jari 2964  
2 26 Feb 07 jari 2965
2 26 Feb 07 jari 2966       int row=0; 
2 26 Feb 07 jari 2967       int column =0;
2 26 Feb 07 jari 2968       double error = 1E-5;
2 26 Feb 07 jari 2969
2 26 Feb 07 jari 2970       if (Math.abs(alpha-0.10) < error) {
2 26 Feb 07 jari 2971          column =1;
2 26 Feb 07 jari 2972       } else if (Math.abs(alpha-0.05) < error)  {
2 26 Feb 07 jari 2973          column =2;
2 26 Feb 07 jari 2974       } else if (Math.abs(alpha-0.02) < error)  {
2 26 Feb 07 jari 2975          column =3;
2 26 Feb 07 jari 2976       } else if (Math.abs(alpha-0.01) < error)  {
2 26 Feb 07 jari 2977          column =4;
2 26 Feb 07 jari 2978       } else if (Math.abs(alpha-0.001) < error)  {
2 26 Feb 07 jari 2979          column =5;
2 26 Feb 07 jari 2980       } else {
2 26 Feb 07 jari 2981          return -2.0;
2 26 Feb 07 jari 2982       }
2 26 Feb 07 jari 2983
2 26 Feb 07 jari 2984      
2 26 Feb 07 jari 2985       if (degreeOfFreedom <= 0) {
2 26 Feb 07 jari 2986     //     System.out.println("DAM.java: getTValue() - degreeOfFreedom <= 0, return");
2 26 Feb 07 jari 2987          return -2.0;
2 26 Feb 07 jari 2988       } else if (degreeOfFreedom <= 30) {
2 26 Feb 07 jari 2989          row = degreeOfFreedom-1;
2 26 Feb 07 jari 2990       } else if ((degreeOfFreedom>30) && (degreeOfFreedom<40)) {
2 26 Feb 07 jari 2991          row = 29;
2 26 Feb 07 jari 2992       } else if ((degreeOfFreedom>=40) && (degreeOfFreedom<60)) {
2 26 Feb 07 jari 2993          row = 30;
2 26 Feb 07 jari 2994       } else if ((degreeOfFreedom>=60) && (degreeOfFreedom<120)) {
2 26 Feb 07 jari 2995          row = 31;
2 26 Feb 07 jari 2996       } else if ((degreeOfFreedom>=120) && (degreeOfFreedom<1000)) {
2 26 Feb 07 jari 2997          row = 32;
2 26 Feb 07 jari 2998       } else if (degreeOfFreedom>=1000) {
2 26 Feb 07 jari 2999          row = 33;
2 26 Feb 07 jari 3000       } else if (degreeOfFreedom <= 0) {
2 26 Feb 07 jari 3001       }
2 26 Feb 07 jari 3002
2 26 Feb 07 jari 3003 /*
2 26 Feb 07 jari 3004       System.out.println("");
2 26 Feb 07 jari 3005       System.out.println("getTValue() - row = " + row);
2 26 Feb 07 jari 3006       System.out.println("getTValue() - column = " + column);
2 26 Feb 07 jari 3007 */
2 26 Feb 07 jari 3008
2 26 Feb 07 jari 3009       return ttbl[row][column];
2 26 Feb 07 jari 3010
2 26 Feb 07 jari 3011     }
2 26 Feb 07 jari 3012
2 26 Feb 07 jari 3013
2 26 Feb 07 jari 3014     /**
2 26 Feb 07 jari 3015      * getSampleMeans(): obtain the mean value for each sample
2 26 Feb 07 jari 3016      *
2 26 Feb 07 jari 3017      */
2 26 Feb 07 jari 3018     private Matrix getSampleMeans(Matrix expMatrix) {
2 26 Feb 07 jari 3019
2 26 Feb 07 jari 3020         int numOfColumns = expMatrix.getColumnDimension();
2 26 Feb 07 jari 3021         int numOfRows = expMatrix.getRowDimension();
2 26 Feb 07 jari 3022
2 26 Feb 07 jari 3023   Matrix means = new Matrix(1, numOfColumns);
2 26 Feb 07 jari 3024
2 26 Feb 07 jari 3025   double currentSum = 0;
2 26 Feb 07 jari 3026         double currentMean =0;
2 26 Feb 07 jari 3027   int denominator = 0;
2 26 Feb 07 jari 3028   double value;
2 26 Feb 07 jari 3029
2 26 Feb 07 jari 3030 //        System.out.println("getSampleMeans() - numOfColumns = " + numOfColumns);
2 26 Feb 07 jari 3031
2 26 Feb 07 jari 3032   for (int column=0; column<numOfColumns; column++) {
2 26 Feb 07 jari 3033       currentSum = 0.0;
2 26 Feb 07 jari 3034       denominator = 0;
2 26 Feb 07 jari 3035       for (int row=0; row<numOfRows; row++) {
2 26 Feb 07 jari 3036     value = expMatrix.get(row, column);
2 26 Feb 07 jari 3037     if (!Double.isNaN(value)) {
2 26 Feb 07 jari 3038         currentSum += value;
2 26 Feb 07 jari 3039         denominator++;
2 26 Feb 07 jari 3040     }
2 26 Feb 07 jari 3041       }
2 26 Feb 07 jari 3042
2 26 Feb 07 jari 3043 //            System.out.println("getSampleMeans() - currentSum = " + currentSum + ", denominator = " + denominator);
2 26 Feb 07 jari 3044
2 26 Feb 07 jari 3045             currentMean = (currentSum/denominator);
2 26 Feb 07 jari 3046       means.set(0, column, currentMean);
2 26 Feb 07 jari 3047   }
2 26 Feb 07 jari 3048
2 26 Feb 07 jari 3049 /*
2 26 Feb 07 jari 3050         for (int i=0; i<numOfColumns; i++) { 
2 26 Feb 07 jari 3051             System.out.println("getSampleMeans() - mean[" + i + "] = " + means.get(0,i));
2 26 Feb 07 jari 3052         }
2 26 Feb 07 jari 3053         System.out.println(" ");
2 26 Feb 07 jari 3054 */
2 26 Feb 07 jari 3055
2 26 Feb 07 jari 3056   return means;
2 26 Feb 07 jari 3057     }
2 26 Feb 07 jari 3058    
2 26 Feb 07 jari 3059
2 26 Feb 07 jari 3060     private Matrix getSampleVariances(Matrix expMatrix, Matrix means) {
2 26 Feb 07 jari 3061   final int rows = means.getRowDimension();
2 26 Feb 07 jari 3062   final int columns = means.getColumnDimension();
2 26 Feb 07 jari 3063   Matrix variances = new Matrix(rows, columns);
2 26 Feb 07 jari 3064   for (int row=0; row<rows; row++) {
2 26 Feb 07 jari 3065       for (int column=0; column<columns; column++) {
2 26 Feb 07 jari 3066     variances.set(row, column, getSampleVariance(expMatrix, means));
2 26 Feb 07 jari 3067       }
2 26 Feb 07 jari 3068   }
2 26 Feb 07 jari 3069   return variances;
2 26 Feb 07 jari 3070     }
2 26 Feb 07 jari 3071   
2 26 Feb 07 jari 3072
2 26 Feb 07 jari 3073     int denominator = 0;
2 26 Feb 07 jari 3074     /**
2 26 Feb 07 jari 3075      * getSampleNormalizedSum(): obtain the normalized sum value for all samples
2 26 Feb 07 jari 3076      *
2 26 Feb 07 jari 3077      */
2 26 Feb 07 jari 3078     private double getSampleNormalizedSum(Matrix expMatrix, Matrix means) {
2 26 Feb 07 jari 3079
2 26 Feb 07 jari 3080   double sum = 0;
2 26 Feb 07 jari 3081   double value = 0;
2 26 Feb 07 jari 3082
2 26 Feb 07 jari 3083         int numOfGenes = expMatrix.getRowDimension();
2 26 Feb 07 jari 3084         int numOfSamples = expMatrix.getColumnDimension();
2 26 Feb 07 jari 3085
2 26 Feb 07 jari 3086         for (int sample=0; sample<numOfSamples; sample++) {
2 26 Feb 07 jari 3087
2 26 Feb 07 jari 3088             double mean = means.get(0,sample);
2 26 Feb 07 jari 3089
2 26 Feb 07 jari 3090       for (int row=0; row<numOfGenes; row++) {
2 26 Feb 07 jari 3091           value = expMatrix.get(row, sample);
2 26 Feb 07 jari 3092           if (!Double.isNaN(value)) {
2 26 Feb 07 jari 3093         sum += Math.pow(value - mean, 2);
2 26 Feb 07 jari 3094         denominator++;
2 26 Feb 07 jari 3095           }
2 26 Feb 07 jari 3096             }
2 26 Feb 07 jari 3097   }
2 26 Feb 07 jari 3098   return sum;
2 26 Feb 07 jari 3099     }
2 26 Feb 07 jari 3100
2 26 Feb 07 jari 3101     /**
2 26 Feb 07 jari 3102      * getSampleVariance(): obtain the sample variance for all samples
2 26 Feb 07 jari 3103      *
2 26 Feb 07 jari 3104      */
2 26 Feb 07 jari 3105     private double getSampleVariance(Matrix expMatrix, Matrix means) {
2 26 Feb 07 jari 3106
2 26 Feb 07 jari 3107         int denominator = 0;
2 26 Feb 07 jari 3108         double sum = 0;
2 26 Feb 07 jari 3109   double value = 0;
2 26 Feb 07 jari 3110
2 26 Feb 07 jari 3111         int numOfSamples = expMatrix.getColumnDimension();
2 26 Feb 07 jari 3112         int numOfGenes = expMatrix.getRowDimension();
2 26 Feb 07 jari 3113
2 26 Feb 07 jari 3114         for (int sample=0; sample<numOfSamples; sample++) {
2 26 Feb 07 jari 3115
2 26 Feb 07 jari 3116             double mean = means.get(0,sample);
2 26 Feb 07 jari 3117
2 26 Feb 07 jari 3118       for (int row=0; row<numOfGenes; row++) {
2 26 Feb 07 jari 3119           value = expMatrix.get(row, sample);
2 26 Feb 07 jari 3120           if (!Double.isNaN(value)) {
2 26 Feb 07 jari 3121         sum += Math.pow(value - mean, 2);
2 26 Feb 07 jari 3122         denominator++;
2 26 Feb 07 jari 3123           }
2 26 Feb 07 jari 3124             }
2 26 Feb 07 jari 3125   }
2 26 Feb 07 jari 3126
2 26 Feb 07 jari 3127 /*
2 26 Feb 07 jari 3128         System.out.println("getSampleVariance() - sum = " + sum);
2 26 Feb 07 jari 3129         System.out.println("getSampleVariance() - denominator = " + denominator);
2 26 Feb 07 jari 3130         System.out.println("getSampleVariance() - sampleVariance = " + Math.sqrt(sum/(denominator-1)));
2 26 Feb 07 jari 3131         System.out.println(" ");
2 26 Feb 07 jari 3132 */
2 26 Feb 07 jari 3133
2 26 Feb 07 jari 3134   return Math.sqrt(sum/(denominator-1));
2 26 Feb 07 jari 3135     }    
2 26 Feb 07 jari 3136  
2 26 Feb 07 jari 3137     /** 
2 26 Feb 07 jari 3138      * getCovarianceMatrix(): obtain the covarianceMatrix
2 26 Feb 07 jari 3139      *
2 26 Feb 07 jari 3140      */
2 26 Feb 07 jari 3141     private Matrix getCovarianceMatrix(Matrix expMatrix, Matrix means) {
2 26 Feb 07 jari 3142
2 26 Feb 07 jari 3143   final int rows = expMatrix.getRowDimension();
2 26 Feb 07 jari 3144   final int columns = expMatrix.getColumnDimension();
2 26 Feb 07 jari 3145
2 26 Feb 07 jari 3146         Matrix [] X = new Matrix[columns];
2 26 Feb 07 jari 3147
2 26 Feb 07 jari 3148         for (int column=0; column<columns; column++) {
2 26 Feb 07 jari 3149             X[column] = new Matrix(1, rows);
2 26 Feb 07 jari 3150             X[column] = (expMatrix.getMatrix(0, rows-1, column, column)).transpose(); 
2 26 Feb 07 jari 3151
2 26 Feb 07 jari 3152 /*
2 26 Feb 07 jari 3153       System.out.println("getCovarianceMatrix() -- X[" + column + "] Matrix : ");
2 26 Feb 07 jari 3154       for (int i=0; i<X[column].getColumnDimension(); i++) {
2 26 Feb 07 jari 3155     System.out.print(X[column].get(0, i) + ", ");
2 26 Feb 07 jari 3156       } 
2 26 Feb 07 jari 3157       System.out.println(" ");
2 26 Feb 07 jari 3158       System.out.println(" ");
2 26 Feb 07 jari 3159 */
2 26 Feb 07 jari 3160         }
2 26 Feb 07 jari 3161
2 26 Feb 07 jari 3162         Matrix covariance = new Matrix(rows, rows);
2 26 Feb 07 jari 3163
2 26 Feb 07 jari 3164 /*
2 26 Feb 07 jari 3165   System.out.println("getCovarianceMatrix() -- mean matrix: ");
2 26 Feb 07 jari 3166   for (int i=0; i<means.getColumnDimension(); i++) {
2 26 Feb 07 jari 3167       System.out.print(means.get(0, i) + ", ");
2 26 Feb 07 jari 3168   } 
2 26 Feb 07 jari 3169   System.out.println(" ");
2 26 Feb 07 jari 3170   System.out.println(" ");
2 26 Feb 07 jari 3171 */
2 26 Feb 07 jari 3172
2 26 Feb 07 jari 3173
2 26 Feb 07 jari 3174         for (int column=0; column < columns; column++) {
2 26 Feb 07 jari 3175             covariance = covariance.plus(((X[column].minus(means)).transpose()).times((X[column].minus(means)))); 
2 26 Feb 07 jari 3176         }
2 26 Feb 07 jari 3177
2 26 Feb 07 jari 3178         covariance = covariance.times(1.0/(double)(columns-1));
2 26 Feb 07 jari 3179
2 26 Feb 07 jari 3180 /*
2 26 Feb 07 jari 3181         System.out.println("getCovarianceMatrix() -- covariance matrix: ");
2 26 Feb 07 jari 3182         for (int row=0; row<covariance.getRowDimension(); row++) {
2 26 Feb 07 jari 3183             for (int column=0; column<covariance.getColumnDimension(); column++) {
2 26 Feb 07 jari 3184                 System.out.print(covariance.get(row, column) + ", ");
2 26 Feb 07 jari 3185             }
2 26 Feb 07 jari 3186             System.out.println(" ");
2 26 Feb 07 jari 3187         } 
2 26 Feb 07 jari 3188         System.out.println(" ");
2 26 Feb 07 jari 3189 */
2 26 Feb 07 jari 3190
2 26 Feb 07 jari 3191   return covariance;
2 26 Feb 07 jari 3192     }
2 26 Feb 07 jari 3193
2 26 Feb 07 jari 3194
2 26 Feb 07 jari 3195     /**
2 26 Feb 07 jari 3196      * getGroupMean(): obtain the mean value for a given group
2 26 Feb 07 jari 3197      *
2 26 Feb 07 jari 3198      */
2 26 Feb 07 jari 3199     private double getGroupMean(double[] group) {
2 26 Feb 07 jari 3200   double sum = 0;
2 26 Feb 07 jari 3201   int n = 0;
2 26 Feb 07 jari 3202   for (int i = 0; i < group.length; i++) {
2 26 Feb 07 jari 3203       if (!Double.isNaN(group[i])) {
2 26 Feb 07 jari 3204     sum = sum + group[i];
2 26 Feb 07 jari 3205     n++;
2 26 Feb 07 jari 3206       }
2 26 Feb 07 jari 3207   }
2 26 Feb 07 jari 3208   
2 26 Feb 07 jari 3209   if (n == 0) {
2 26 Feb 07 jari 3210             return Double.NaN;
2 26 Feb 07 jari 3211         }
2 26 Feb 07 jari 3212   double mean =  sum / n;
2 26 Feb 07 jari 3213         
2 26 Feb 07 jari 3214         if (Double.isInfinite(mean)) {
2 26 Feb 07 jari 3215             return Double.NaN;
2 26 Feb 07 jari 3216         }
2 26 Feb 07 jari 3217         
2 26 Feb 07 jari 3218   return mean;
2 26 Feb 07 jari 3219     }
2 26 Feb 07 jari 3220
2 26 Feb 07 jari 3221
2 26 Feb 07 jari 3222     private boolean isFoundInVector(int element, Vector vect) {
2 26 Feb 07 jari 3223         boolean found = false;
2 26 Feb 07 jari 3224         for (int i = 0; i < vect.size(); i++) {
2 26 Feb 07 jari 3225             if (element == ((Integer)(vect.get(i))).intValue()) {
2 26 Feb 07 jari 3226                 found = true;
2 26 Feb 07 jari 3227                 break;
2 26 Feb 07 jari 3228             }
2 26 Feb 07 jari 3229         }
2 26 Feb 07 jari 3230         return found;
2 26 Feb 07 jari 3231     }
2 26 Feb 07 jari 3232
2 26 Feb 07 jari 3233
2 26 Feb 07 jari 3234     public Matrix getJamaMatrix(FloatMatrix floatMatrix) {
2 26 Feb 07 jari 3235
2 26 Feb 07 jari 3236         int rows = floatMatrix.getRowDimension();
2 26 Feb 07 jari 3237         int columns = floatMatrix.getColumnDimension();
2 26 Feb 07 jari 3238
2 26 Feb 07 jari 3239         double[][] values = new double[rows][columns];
2 26 Feb 07 jari 3240
2 26 Feb 07 jari 3241         for(int i=0; i<rows; i++) {
2 26 Feb 07 jari 3242             for(int j=0; j<columns; j++) {
2 26 Feb 07 jari 3243                 values[i][j] = (double) floatMatrix.A[i][j];
2 26 Feb 07 jari 3244             }
2 26 Feb 07 jari 3245         }
2 26 Feb 07 jari 3246
2 26 Feb 07 jari 3247         Matrix jamaMatrix = new Matrix(values);
2 26 Feb 07 jari 3248
2 26 Feb 07 jari 3249         return jamaMatrix;
2 26 Feb 07 jari 3250     }
2 26 Feb 07 jari 3251
2 26 Feb 07 jari 3252
2 26 Feb 07 jari 3253     public FloatMatrix getFloatMatrix(Matrix doubleMatrix) {
2 26 Feb 07 jari 3254         int rows = doubleMatrix.getRowDimension();
2 26 Feb 07 jari 3255         int columns = doubleMatrix.getColumnDimension();
2 26 Feb 07 jari 3256
2 26 Feb 07 jari 3257         float[][] values = new float[rows][columns];
2 26 Feb 07 jari 3258
2 26 Feb 07 jari 3259         for(int i=0; i<rows; i++) {
2 26 Feb 07 jari 3260             for(int j=0; j<columns; j++) {
2 26 Feb 07 jari 3261                 values[i][j] = (float) ((doubleMatrix.getArray())[i][j]);
2 26 Feb 07 jari 3262             }
2 26 Feb 07 jari 3263         }
2 26 Feb 07 jari 3264
2 26 Feb 07 jari 3265         FloatMatrix floatMatrix = new FloatMatrix(values);
2 26 Feb 07 jari 3266
2 26 Feb 07 jari 3267         return floatMatrix;
2 26 Feb 07 jari 3268     }
2 26 Feb 07 jari 3269
2 26 Feb 07 jari 3270
2 26 Feb 07 jari 3271     private Matrix getMeansForGenes(Vector[] clusters) {
2 26 Feb 07 jari 3272   Matrix means = new Matrix(clusters.length, numberOfSamples);
2 26 Feb 07 jari 3273   Matrix mean;
2 26 Feb 07 jari 3274   for (int i=0; i<clusters.length; i++) {
2 26 Feb 07 jari 3275       mean = getMeanForGenes(clusters[i]);
2 26 Feb 07 jari 3276             means.setMatrix(i, i, 0, numberOfSamples-1, mean);
2 26 Feb 07 jari 3277   }
2 26 Feb 07 jari 3278   return means;
2 26 Feb 07 jari 3279     }
2 26 Feb 07 jari 3280     
2 26 Feb 07 jari 3281     private Matrix getMeanForGenes(Vector cluster) {
2 26 Feb 07 jari 3282   Matrix mean = new Matrix(1, numberOfSamples);
2 26 Feb 07 jari 3283   double currentMean;
2 26 Feb 07 jari 3284   int n = cluster.size();
2 26 Feb 07 jari 3285   int denom = 0;
2 26 Feb 07 jari 3286   double value;
2 26 Feb 07 jari 3287   for (int i=0; i<numberOfSamples; i++) {
2 26 Feb 07 jari 3288       currentMean = 0.0d;
2 26 Feb 07 jari 3289       denom  = 0;
2 26 Feb 07 jari 3290       for (int j=0; j<n; j++) {
2 26 Feb 07 jari 3291     value = expMatrix.get(((Integer) cluster.get(j)).intValue(), i);
2 26 Feb 07 jari 3292     if (!Double.isNaN(value)) {
2 26 Feb 07 jari 3293         currentMean += value;
2 26 Feb 07 jari 3294         denom++;
2 26 Feb 07 jari 3295     }
2 26 Feb 07 jari 3296       }
2 26 Feb 07 jari 3297       mean.set(0, i, currentMean/(double)denom);
2 26 Feb 07 jari 3298   }
2 26 Feb 07 jari 3299   return mean;
2 26 Feb 07 jari 3300     }
2 26 Feb 07 jari 3301
2 26 Feb 07 jari 3302     private Matrix getVariancesForGenes(Vector[] clusters, Matrix means) {
2 26 Feb 07 jari 3303   final int rows = means.getRowDimension();
2 26 Feb 07 jari 3304   final int columns = means.getColumnDimension();
2 26 Feb 07 jari 3305   Matrix variances = new Matrix(rows, columns);
2 26 Feb 07 jari 3306   for (int row=0; row<rows; row++) {
2 26 Feb 07 jari 3307       for (int column=0; column<columns; column++) {
2 26 Feb 07 jari 3308     variances.set(row, column, getClusterVarianceForGenes(clusters[row], column, means.get(row, column)));
2 26 Feb 07 jari 3309       }
2 26 Feb 07 jari 3310   }
2 26 Feb 07 jari 3311   return variances;
2 26 Feb 07 jari 3312     }
2 26 Feb 07 jari 3313
2 26 Feb 07 jari 3314     int validN;
2 26 Feb 07 jari 3315     private double getSampleNormalizedSumForGenes(Vector cluster, int column, double mean) {
2 26 Feb 07 jari 3316   final int size = cluster.size();
2 26 Feb 07 jari 3317   double sum = 0f;
2 26 Feb 07 jari 3318   validN = 0;
2 26 Feb 07 jari 3319   double value;
2 26 Feb 07 jari 3320   for (int i=0; i<size; i++) {
2 26 Feb 07 jari 3321       value = expMatrix.get(((Integer) cluster.get(i)).intValue(), column);
2 26 Feb 07 jari 3322       if (!Double.isNaN(value)) {
2 26 Feb 07 jari 3323     sum += Math.pow(value-mean, 2);
2 26 Feb 07 jari 3324     validN++;
2 26 Feb 07 jari 3325       }
2 26 Feb 07 jari 3326   }
2 26 Feb 07 jari 3327   return sum;
2 26 Feb 07 jari 3328     }
2 26 Feb 07 jari 3329  
2 26 Feb 07 jari 3330     private double getClusterVarianceForGenes(Vector cluster, int column, double mean) {
2 26 Feb 07 jari 3331   return Math.sqrt(getSampleNormalizedSumForGenes(cluster, column, mean)/(validN-1));
2 26 Feb 07 jari 3332     }    
2 26 Feb 07 jari 3333
2 26 Feb 07 jari 3334
2 26 Feb 07 jari 3335     private Matrix getMeans(Vector[] clusters) {
2 26 Feb 07 jari 3336   Matrix means = new Matrix(clusters.length, numberOfGenes);
2 26 Feb 07 jari 3337   Matrix mean;
2 26 Feb 07 jari 3338   for (int i=0; i<clusters.length; i++) {
2 26 Feb 07 jari 3339       mean = getMean(clusters[i]);
2 26 Feb 07 jari 3340             means.setMatrix(i, i, 0, numberOfGenes-1, mean);
2 26 Feb 07 jari 3341   }
2 26 Feb 07 jari 3342   return means;
2 26 Feb 07 jari 3343     }
2 26 Feb 07 jari 3344     
2 26 Feb 07 jari 3345     private Matrix getMean(Vector cluster) {
2 26 Feb 07 jari 3346   Matrix mean = new Matrix(1, numberOfGenes);
2 26 Feb 07 jari 3347   double currentMean;
2 26 Feb 07 jari 3348   int n = cluster.size();
2 26 Feb 07 jari 3349   int denom = 0;
2 26 Feb 07 jari 3350   double value;
2 26 Feb 07 jari 3351   for (int i=0; i<numberOfGenes; i++) {
2 26 Feb 07 jari 3352       currentMean = 0.0d;
2 26 Feb 07 jari 3353       denom  = 0;
2 26 Feb 07 jari 3354       for (int j=0; j<n; j++) {
2 26 Feb 07 jari 3355     value = expMatrixTranspose.get(((Integer) cluster.get(j)).intValue(), i);
2 26 Feb 07 jari 3356     if (!Double.isNaN(value)) {
2 26 Feb 07 jari 3357         currentMean += value;
2 26 Feb 07 jari 3358         denom++;
2 26 Feb 07 jari 3359     }
2 26 Feb 07 jari 3360       }
2 26 Feb 07 jari 3361       mean.set(0, i, currentMean/(double)denom);
2 26 Feb 07 jari 3362   }
2 26 Feb 07 jari 3363   return mean;
2 26 Feb 07 jari 3364     }
2 26 Feb 07 jari 3365
2 26 Feb 07 jari 3366     private Matrix getVariances(Vector[] clusters, Matrix means) {
2 26 Feb 07 jari 3367   final int rows = means.getRowDimension();
2 26 Feb 07 jari 3368   final int columns = means.getColumnDimension();
2 26 Feb 07 jari 3369   Matrix variances = new Matrix(rows, columns);
2 26 Feb 07 jari 3370   for (int row=0; row<rows; row++) {
2 26 Feb 07 jari 3371       for (int column=0; column<columns; column++) {
2 26 Feb 07 jari 3372     variances.set(row, column, getClusterVariance(clusters[row], column, means.get(row, column)));
2 26 Feb 07 jari 3373       }
2 26 Feb 07 jari 3374   }
2 26 Feb 07 jari 3375   return variances;
2 26 Feb 07 jari 3376     }
2 26 Feb 07 jari 3377
2 26 Feb 07 jari 3378  
2 26 Feb 07 jari 3379     private double getSampleNormalizedSum(Vector cluster, int column, double mean) {
2 26 Feb 07 jari 3380   final int size = cluster.size();
2 26 Feb 07 jari 3381   double sum = 0f;
2 26 Feb 07 jari 3382   validN = 0;
2 26 Feb 07 jari 3383   double value;
2 26 Feb 07 jari 3384   for (int i=0; i<size; i++) {
2 26 Feb 07 jari 3385       value = expMatrixTranspose.get(((Integer) cluster.get(i)).intValue(), column);
2 26 Feb 07 jari 3386       if (!Double.isNaN(value)) {
2 26 Feb 07 jari 3387     sum += Math.pow(value-mean, 2);
2 26 Feb 07 jari 3388     validN++;
2 26 Feb 07 jari 3389       }
2 26 Feb 07 jari 3390   }
2 26 Feb 07 jari 3391   return sum;
2 26 Feb 07 jari 3392     }
2 26 Feb 07 jari 3393  
2 26 Feb 07 jari 3394     private double getClusterVariance(Vector cluster, int column, double mean) {
2 26 Feb 07 jari 3395   return Math.sqrt(getSampleNormalizedSum(cluster, column, mean)/(validN-1));
2 26 Feb 07 jari 3396     }    
2 26 Feb 07 jari 3397    
2 26 Feb 07 jari 3398
2 26 Feb 07 jari 3399     private int[] convert2int(Vector source) {
2 26 Feb 07 jari 3400   int[] int_matrix = new int[source.size()];
2 26 Feb 07 jari 3401   for (int i=0; i<int_matrix.length; i++) {
2 26 Feb 07 jari 3402       int_matrix[i] = ((Integer) source.get(i)).intValue();
2 26 Feb 07 jari 3403   }
2 26 Feb 07 jari 3404   return int_matrix;
2 26 Feb 07 jari 3405     }    
2 26 Feb 07 jari 3406
2 26 Feb 07 jari 3407
2 26 Feb 07 jari 3408     private void printMatrix(Matrix matrix) {
2 26 Feb 07 jari 3409
2 26 Feb 07 jari 3410             DecimalFormat df = new DecimalFormat("####.000000");
2 26 Feb 07 jari 3411
2 26 Feb 07 jari 3412       for (int row=0; row<matrix.getRowDimension(); row++) {
2 26 Feb 07 jari 3413     for (int column=0; column<matrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 3414       System.out.print(df.format(matrix.get(row, column)) + "  ");
2 26 Feb 07 jari 3415     }
2 26 Feb 07 jari 3416     System.out.println(" ");
2 26 Feb 07 jari 3417       } 
2 26 Feb 07 jari 3418       System.out.println(" ");
2 26 Feb 07 jari 3419       System.out.println(" ");
2 26 Feb 07 jari 3420     }
2 26 Feb 07 jari 3421
2 26 Feb 07 jari 3422     private void printMatrix(FloatMatrix matrix) {
2 26 Feb 07 jari 3423
2 26 Feb 07 jari 3424             DecimalFormat df = new DecimalFormat("####.000000");
2 26 Feb 07 jari 3425
2 26 Feb 07 jari 3426       for (int row=0; row<matrix.getRowDimension(); row++) {
2 26 Feb 07 jari 3427     for (int column=0; column<matrix.getColumnDimension(); column++) {
2 26 Feb 07 jari 3428       System.out.print(df.format(matrix.get(row, column)) + "  ");
2 26 Feb 07 jari 3429     }
2 26 Feb 07 jari 3430     System.out.println(" ");
2 26 Feb 07 jari 3431       } 
2 26 Feb 07 jari 3432       System.out.println(" ");
2 26 Feb 07 jari 3433       System.out.println(" ");
2 26 Feb 07 jari 3434     }
2 26 Feb 07 jari 3435
2 26 Feb 07 jari 3436     public void abort() {
2 26 Feb 07 jari 3437   stop = true;
2 26 Feb 07 jari 3438     }
2 26 Feb 07 jari 3439 }

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