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

Code
Comments
Other
Rev Date Author Line
2 26 Feb 07 jari 1 /*
2 26 Feb 07 jari 2 Copyright @ 1999-2003, The Institute for Genomic Research (TIGR).
2 26 Feb 07 jari 3 All rights reserved.
2 26 Feb 07 jari 4 */
2 26 Feb 07 jari 5 /*
2 26 Feb 07 jari 6  * $RCSfile: PermutationTest.java,v $
2 26 Feb 07 jari 7  * $Revision: 1.4 $
2 26 Feb 07 jari 8  * $Date: 2006/02/23 20:59:45 $
2 26 Feb 07 jari 9  * $Author: caliente $
2 26 Feb 07 jari 10  * $State: Exp $
2 26 Feb 07 jari 11  */
2 26 Feb 07 jari 12 package org.tigr.microarray.mev.cluster.algorithm.impl;
2 26 Feb 07 jari 13
2 26 Feb 07 jari 14 import java.util.Arrays;
2 26 Feb 07 jari 15
2 26 Feb 07 jari 16 import org.tigr.microarray.mev.cluster.algorithm.AbortException;
2 26 Feb 07 jari 17 import org.tigr.microarray.mev.cluster.algorithm.AbstractAlgorithm;
2 26 Feb 07 jari 18 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmData;
2 26 Feb 07 jari 19 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmEvent;
2 26 Feb 07 jari 20 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmException;
2 26 Feb 07 jari 21 import org.tigr.microarray.mev.cluster.algorithm.AlgorithmParameters;
2 26 Feb 07 jari 22 import org.tigr.microarray.mev.cluster.algorithm.impl.util.IntSorter;
2 26 Feb 07 jari 23 import org.tigr.util.FloatMatrix;
2 26 Feb 07 jari 24
2 26 Feb 07 jari 25 public class PermutationTest extends AbstractAlgorithm {
2 26 Feb 07 jari 26
2 26 Feb 07 jari 27     private static final int c_DecileCount = 10;
2 26 Feb 07 jari 28     public static final double LOG2 = Math.log(2.0);
2 26 Feb 07 jari 29     public static final float MINVAL = 0f;
2 26 Feb 07 jari 30     public static final float MAXVAL = 1f;
2 26 Feb 07 jari 31
2 26 Feb 07 jari 32     private double getEntropy(float[] pVector) {
2 26 Feb 07 jari 33         double fltMin = Double.MAX_VALUE;
2 26 Feb 07 jari 34         double fltMax = -Double.MAX_VALUE;
2 26 Feb 07 jari 35         int i=0;
2 26 Feb 07 jari 36         int[] arrDeciles = new int[c_DecileCount];
2 26 Feb 07 jari 37
2 26 Feb 07 jari 38         final int iSize = pVector.length;
2 26 Feb 07 jari 39         int iValCount = 0;
2 26 Feb 07 jari 40         for (i=0; i<iSize; i++) {
2 26 Feb 07 jari 41             if (Double.isNaN(pVector[i]))
2 26 Feb 07 jari 42                 continue;
2 26 Feb 07 jari 43             fltMin = Math.min(fltMin, pVector[i]);
2 26 Feb 07 jari 44             fltMax = Math.max(fltMax, pVector[i]);
2 26 Feb 07 jari 45             iValCount++;
2 26 Feb 07 jari 46         }
2 26 Feb 07 jari 47
2 26 Feb 07 jari 48         double fltStep = (fltMax-fltMin)/(c_DecileCount);
2 26 Feb 07 jari 49         if (fltStep == 0d) {
2 26 Feb 07 jari 50             return -1.0*Math.log(1.0)/LOG2;
2 26 Feb 07 jari 51         }
2 26 Feb 07 jari 52
2 26 Feb 07 jari 53         if (fltMin == Double.MAX_VALUE)
2 26 Feb 07 jari 54             return 0d;
2 26 Feb 07 jari 55
2 26 Feb 07 jari 56         Arrays.fill(arrDeciles, 0);
2 26 Feb 07 jari 57         for (i=0; i<iSize; i++) {
2 26 Feb 07 jari 58             if (Double.isNaN(pVector[i]))
2 26 Feb 07 jari 59                 continue;
2 26 Feb 07 jari 60             int iDecileInd = (int)Math.ceil((pVector[i]-fltMin)/fltStep)-1;
2 26 Feb 07 jari 61             if (iDecileInd < 0) {
2 26 Feb 07 jari 62                 iDecileInd = 0;
2 26 Feb 07 jari 63             }
2 26 Feb 07 jari 64             arrDeciles[iDecileInd]++;
2 26 Feb 07 jari 65         }
2 26 Feb 07 jari 66         if (iValCount == 0)
2 26 Feb 07 jari 67             return 0d;
2 26 Feb 07 jari 68
2 26 Feb 07 jari 69         // finally, calculate entropy
2 26 Feb 07 jari 70         double dblEntropy=0;
2 26 Feb 07 jari 71
2 26 Feb 07 jari 72         for (i=0; i<c_DecileCount; i++) {
2 26 Feb 07 jari 73             if (arrDeciles[i] == 0) {
2 26 Feb 07 jari 74                 continue;
2 26 Feb 07 jari 75             }
2 26 Feb 07 jari 76             double dblPx=((double)arrDeciles[i])/iValCount;
2 26 Feb 07 jari 77             dblEntropy += dblPx*Math.log(dblPx)/LOG2; // log2(x)==log(x)/log(2)
2 26 Feb 07 jari 78         }
2 26 Feb 07 jari 79         return -dblEntropy;
2 26 Feb 07 jari 80     }
2 26 Feb 07 jari 81
2 26 Feb 07 jari 82     private boolean stop = false;
2 26 Feb 07 jari 83
2 26 Feb 07 jari 84     private int[]   m_arrMainHisto = null;  // Real histo stores here
2 26 Feb 07 jari 85     //private int[][] m_arrPermHistos = null; // for every i-th permutaion m_PermHistos[i] stores appropriate histogram
2 26 Feb 07 jari 86     private double[] m_arrAvgHisto=null;
2 26 Feb 07 jari 87     private int     m_iHistoSize = 0;
2 26 Feb 07 jari 88     private int     m_iPermutationSize = 0;
2 26 Feb 07 jari 89     private float   m_fltHistoStep = 0.5f;
2 26 Feb 07 jari 90
2 26 Feb 07 jari 91     private float   m_fltMaxStatSignificant = 0;
2 26 Feb 07 jari 92
2 26 Feb 07 jari 93     private void Assert(int iFrom, int iTo, float fltVal) {
2 26 Feb 07 jari 94         m_arrMainHisto[GetIndByVal(fltVal)]++;
2 26 Feb 07 jari 95     }
2 26 Feb 07 jari 96
2 26 Feb 07 jari 97     private void AssertPermutted(int iPermutationNum, int iFrom, int iTo, float fltVal) {
2 26 Feb 07 jari 98         //m_arrPermHistos[iPermutationNum][GetIndByVal(fltVal)]++;
2 26 Feb 07 jari 99         m_arrAvgHisto[GetIndByVal(fltVal)]++;
2 26 Feb 07 jari 100         m_fltMaxStatSignificant=Math.max(m_fltMaxStatSignificant, fltVal);
2 26 Feb 07 jari 101     }
2 26 Feb 07 jari 102
2 26 Feb 07 jari 103     private int GetIndByVal(float fltVal) {
2 26 Feb 07 jari 104         int iDecileInd = (int)Math.ceil((fltVal-(MINVAL))/m_fltHistoStep)-1;
2 26 Feb 07 jari 105         if (iDecileInd < 0)
2 26 Feb 07 jari 106             return 0;
2 26 Feb 07 jari 107         return iDecileInd;
2 26 Feb 07 jari 108     }
2 26 Feb 07 jari 109
2 26 Feb 07 jari 110     private void RandomPermute(FloatMatrix matr, int iRow) {
2 26 Feb 07 jari 111         int iFirstInd = 0;
2 26 Feb 07 jari 112         int iSecondInd =0;
2 26 Feb 07 jari 113         int iPermCount=(int)(Math.random()*matr.getColumnDimension()+1);
2 26 Feb 07 jari 114         for(int i=0; i<iPermCount;i++){
2 26 Feb 07 jari 115           do{
2 26 Feb 07 jari 116             iFirstInd = (int)(Math.random()*(matr.getColumnDimension()-1));
2 26 Feb 07 jari 117             iSecondInd = (int)(Math.random()*(matr.getColumnDimension()-1));
2 26 Feb 07 jari 118           }while(iFirstInd==iSecondInd);
2 26 Feb 07 jari 119
2 26 Feb 07 jari 120       //System.out.println(iFirstInd);
2 26 Feb 07 jari 121           float temp=matr.A[iRow][iFirstInd];
2 26 Feb 07 jari 122           matr.A[iRow][iFirstInd]=matr.A[iRow][iSecondInd];
2 26 Feb 07 jari 123           matr.A[iRow][iSecondInd]=temp;
2 26 Feb 07 jari 124         }
2 26 Feb 07 jari 125     }
2 26 Feb 07 jari 126
2 26 Feb 07 jari 127     private void Run(int[] entropyIndices, FloatMatrix expMatrix, int function, float factor, boolean absolute) throws AlgorithmException {
2 26 Feb 07 jari 128         int filteredSize = entropyIndices.length;
2 26 Feb 07 jari 129         int progress = 0;
2 26 Feb 07 jari 130         int links = 0;
2 26 Feb 07 jari 131         int sum = filteredSize*(filteredSize+1)/2;
2 26 Feb 07 jari 132         int step = sum/100+1;
2 26 Feb 07 jari 133         float value;
2 26 Feb 07 jari 134         AlgorithmEvent event = new AlgorithmEvent(this, AlgorithmEvent.SET_UNITS, 100, "Calculating Histogram");
2 26 Feb 07 jari 135         fireValueChanged(event);
2 26 Feb 07 jari 136
2 26 Feb 07 jari 137         event.setId(AlgorithmEvent.PROGRESS_VALUE);
2 26 Feb 07 jari 138
2 26 Feb 07 jari 139         FloatMatrix temp= new FloatMatrix (1,expMatrix.getColumnDimension());
2 26 Feb 07 jari 140         for (int nCurIndOuter=0; nCurIndOuter<filteredSize; nCurIndOuter++) {
2 26 Feb 07 jari 141             if (this.stop) {
2 26 Feb 07 jari 142                 throw new AbortException();
2 26 Feb 07 jari 143             }
2 26 Feb 07 jari 144             progress++;
2 26 Feb 07 jari 145             // calculation
2 26 Feb 07 jari 146            { float[] arrX = expMatrix.A[entropyIndices[nCurIndOuter]];
2 26 Feb 07 jari 147                     float[] arrY = temp.A[0];
2 26 Feb 07 jari 148                     for (int k=0; k<arrY.length; k++)// copy outer vectr
2 26 Feb 07 jari 149                       arrY[k]=arrX[k];
2 26 Feb 07 jari 150            }
2 26 Feb 07 jari 151
2 26 Feb 07 jari 152             for (int nCurIndInner=nCurIndOuter+1; nCurIndInner<filteredSize; nCurIndInner++) {
2 26 Feb 07 jari 153                 //value = AlgorithmUtil.geneDistance(expMatrix, null, entropyIndices[nCurIndOuter], entropyIndices[nCurIndInner], function, factor, absolute);
2 26 Feb 07 jari 154                 //value = value*value; // = abs(r^2)
2 26 Feb 07 jari 155                  //value = AlgorithmUtil.genePearson(temp, expMatrix, 0, entropyIndices[nCurIndInner],1);
2 26 Feb 07 jari 156                  //value*=value;
2 26 Feb 07 jari 157
2 26 Feb 07 jari 158                 //Assert(entropyIndices[nCurIndOuter], entropyIndices[nCurIndInner], value);
2 26 Feb 07 jari 159
2 26 Feb 07 jari 160                 for (int i=0; i<m_iPermutationSize; i++) {
2 26 Feb 07 jari 161                     RandomPermute(temp, 0);
2 26 Feb 07 jari 162                     //value = AlgorithmUtil.geneDistance(temp, expMatrix, 0, entropyIndices[nCurIndInner], function, factor, absolute);
2 26 Feb 07 jari 163                     //value = value*value; // = abs(r^2)
2 26 Feb 07 jari 164                     value = ExperimentUtil.genePearson(temp, expMatrix, 0, entropyIndices[nCurIndInner],1);
2 26 Feb 07 jari 165                     value*=value;
2 26 Feb 07 jari 166                     AssertPermutted(i, entropyIndices[nCurIndOuter], entropyIndices[nCurIndInner], value);
2 26 Feb 07 jari 167                 }
2 26 Feb 07 jari 168
2 26 Feb 07 jari 169                 // progress events handling
2 26 Feb 07 jari 170                 progress++;
2 26 Feb 07 jari 171                 if (progress%step == 0) {
2 26 Feb 07 jari 172                     event.setIntValue(progress/step);
2 26 Feb 07 jari 173                     event.setDescription("Calculating permutation histogram ");
2 26 Feb 07 jari 174                     fireValueChanged(event);
2 26 Feb 07 jari 175                 }
2 26 Feb 07 jari 176             }
2 26 Feb 07 jari 177         }
2 26 Feb 07 jari 178     }
2 26 Feb 07 jari 179
2 26 Feb 07 jari 180     public AlgorithmData execute(AlgorithmData data) throws AlgorithmException {
2 26 Feb 07 jari 181
2 26 Feb 07 jari 182         FloatMatrix expMatrix = data.getMatrix("experiment");
2 26 Feb 07 jari 183         if (expMatrix == null)
2 26 Feb 07 jari 184             return null;
2 26 Feb 07 jari 185
2 26 Feb 07 jari 186         AlgorithmParameters map = data.getParams();
2 26 Feb 07 jari 187         m_iHistoSize = map.getInt("decile-count", 1000);
2 26 Feb 07 jari 188         m_iPermutationSize = map.getInt("permutation-count", 20);
2 26 Feb 07 jari 189         m_fltHistoStep = (MAXVAL-MINVAL)/(float)m_iHistoSize;
2 26 Feb 07 jari 190         m_arrMainHisto = new int[m_iHistoSize];
2 26 Feb 07 jari 191         m_arrAvgHisto=new double[m_iHistoSize];
2 26 Feb 07 jari 192         //m_arrPermHistos = new int[m_iPermutationSize][];
2 26 Feb 07 jari 193         //for (int i=0; i<m_iPermutationSize; i++)
2 26 Feb 07 jari 194         //   m_arrPermHistos[i] = new int[m_iHistoSize];
2 26 Feb 07 jari 195
2 26 Feb 07 jari 196         int function = map.getInt("distance-function", PEARSON);
2 26 Feb 07 jari 197         float factor = map.getFloat("distance-factor", 1.0f);
2 26 Feb 07 jari 198         boolean absolute = map.getBoolean("distance-absolute", true);
2 26 Feb 07 jari 199
2 26 Feb 07 jari 200         int nGenes = expMatrix.getRowDimension();
2 26 Feb 07 jari 201         int filteredSize = nGenes;
2 26 Feb 07 jari 202
2 26 Feb 07 jari 203         // For Entropy Filter
2 26 Feb 07 jari 204         boolean bFilterByEntropy = map.getBoolean("filter-by-entropy");
2 26 Feb 07 jari 205         float fltTopNPercent=map.getFloat("top-n-percent", 100f);
2 26 Feb 07 jari 206         if (fltTopNPercent < 0f || fltTopNPercent>100f) {
2 26 Feb 07 jari 207             throw new AlgorithmException("Filter value is out of range (0, 100)%");
2 26 Feb 07 jari 208         }
2 26 Feb 07 jari 209
2 26 Feb 07 jari 210         int[] entropyIndices = new int[nGenes];
2 26 Feb 07 jari 211         for (int i=0; i<entropyIndices.length; i++) {
2 26 Feb 07 jari 212             entropyIndices[i] = i;
2 26 Feb 07 jari 213         }
2 26 Feb 07 jari 214
2 26 Feb 07 jari 215         if (bFilterByEntropy) {
2 26 Feb 07 jari 216             double[] entropyValues = new double[nGenes];
2 26 Feb 07 jari 217             for (int i=0; i<entropyValues.length; i++) {
2 26 Feb 07 jari 218                 entropyValues[i] = getEntropy(expMatrix.A[i]);
2 26 Feb 07 jari 219             }
2 26 Feb 07 jari 220             IntSorter.sort(entropyIndices, new RelNetComparator(entropyValues));
2 26 Feb 07 jari 221             filteredSize = (int)((float)nGenes*fltTopNPercent/100f);
2 26 Feb 07 jari 222         }
2 26 Feb 07 jari 223
2 26 Feb 07 jari 224         Run(entropyIndices, expMatrix, function, factor, absolute);
2 26 Feb 07 jari 225
2 26 Feb 07 jari 226         // Retrieve threshold from multidimensional histogramms
2 26 Feb 07 jari 227
2 26 Feb 07 jari 228         // make average histo
2 26 Feb 07 jari 229         // TODO: add min histo,max histo
2 26 Feb 07 jari 230
2 26 Feb 07 jari 231         //for (int j=0; j<m_arrAvgHisto.length; j++) {
2 26 Feb 07 jari 232           //  double fltAvg=0;
2 26 Feb 07 jari 233             //double fltMin=Double.MAX_VALUE; // never used in the future, but required for stdDev output in case of displayint a histogramm
2 26 Feb 07 jari 234             //double fltMax=-Double.MAX_VALUE;// never used in the future, but required for stdDev output in case of displayint a histogramm
2 26 Feb 07 jari 235             //for (int i=0; i<m_arrPermHistos.length; i++) {
2 26 Feb 07 jari 236              //   fltAvg+=(double)(m_arrPermHistos[i][j])/(double)m_iPermutationSize;
2 26 Feb 07 jari 237               //  fltMin=Math.min(fltMin, m_arrPermHistos[i][j]);
2 26 Feb 07 jari 238                // fltMax=Math.max(fltMax, m_arrPermHistos[i][j]);
2 26 Feb 07 jari 239             //}
2 26 Feb 07 jari 240             //m_arrAvgHisto[j]=fltAvg;
2 26 Feb 07 jari 241         //}
2 26 Feb 07 jari 242
2 26 Feb 07 jari 243         float fltMinStatSignificantThreshold=m_fltMaxStatSignificant;
2 26 Feb 07 jari 244         if (fltMinStatSignificantThreshold>=1.0){//relaxed method(heuristic)
2 26 Feb 07 jari 245           for (int i=m_arrAvgHisto.length-1; i>0; i--) {
2 26 Feb 07 jari 246             //System.out.println("m_arrAvgHisto["+i+"]="+m_arrAvgHisto[i]);
2 26 Feb 07 jari 247               if (m_arrAvgHisto[i]/m_iPermutationSize>=1f ) {
2 26 Feb 07 jari 248                   fltMinStatSignificantThreshold=(i+1)*m_fltHistoStep;
2 26 Feb 07 jari 249                   if (fltMinStatSignificantThreshold>1f)// just in case
2 26 Feb 07 jari 250                       fltMinStatSignificantThreshold=1f;
2 26 Feb 07 jari 251                   break;
2 26 Feb 07 jari 252               }
2 26 Feb 07 jari 253           }
2 26 Feb 07 jari 254         }
2 26 Feb 07 jari 255
2 26 Feb 07 jari 256
2 26 Feb 07 jari 257         // return the result
2 26 Feb 07 jari 258         AlgorithmData result = new AlgorithmData();
2 26 Feb 07 jari 259         result.addParam("threshold", String.valueOf(fltMinStatSignificantThreshold));
2 26 Feb 07 jari 260         return result;
2 26 Feb 07 jari 261     }
2 26 Feb 07 jari 262
2 26 Feb 07 jari 263     public void abort() {
2 26 Feb 07 jari 264         this.stop = true;
2 26 Feb 07 jari 265     }
2 26 Feb 07 jari 266 }

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