mev-4.0.01/source/org/tigr/microarray/mev/cluster/algorithm/impl/ExperimentUtil.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: ExperimentUtil.java,v $
2 26 Feb 07 jari 7  * $Revision: 1.5 $
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
2 26 Feb 07 jari 13 package org.tigr.microarray.mev.cluster.algorithm.impl;
2 26 Feb 07 jari 14
2 26 Feb 07 jari 15 import org.tigr.microarray.mev.cluster.algorithm.Algorithm;
2 26 Feb 07 jari 16 import org.tigr.util.FloatMatrix;
2 26 Feb 07 jari 17
2 26 Feb 07 jari 18 public class ExperimentUtil {
2 26 Feb 07 jari 19     
2 26 Feb 07 jari 20     public static float distance(FloatMatrix matrix, int e1, int e2, int function, float factor, boolean absolute) {
2 26 Feb 07 jari 21   float result = Float.NaN;
2 26 Feb 07 jari 22   switch (function) {
2 26 Feb 07 jari 23       case Algorithm.PEARSON:
2 26 Feb 07 jari 24     result = pearson(matrix, e1, e2, factor);
2 26 Feb 07 jari 25     factor *= -1;
2 26 Feb 07 jari 26     break;
2 26 Feb 07 jari 27       case Algorithm.COSINE:
2 26 Feb 07 jari 28     result = cosine(matrix, e1, e2, factor);
2 26 Feb 07 jari 29     factor *= -1;
2 26 Feb 07 jari 30     break;
2 26 Feb 07 jari 31       case Algorithm.COVARIANCE:
2 26 Feb 07 jari 32     result = covariance(matrix, e1, e2, factor);
2 26 Feb 07 jari 33     factor *= -1;
2 26 Feb 07 jari 34     break;
2 26 Feb 07 jari 35       case Algorithm.EUCLIDEAN:
2 26 Feb 07 jari 36     result = euclidian(matrix, e1, e2, factor);
2 26 Feb 07 jari 37     break;
2 26 Feb 07 jari 38       case Algorithm.DOTPRODUCT:
2 26 Feb 07 jari 39     result = dotProduct(matrix, e1, e2, factor);
2 26 Feb 07 jari 40     factor *= -1;
2 26 Feb 07 jari 41     break;
2 26 Feb 07 jari 42       case Algorithm.PEARSONUNCENTERED:
2 26 Feb 07 jari 43     result = pearsonUncentered(matrix, e1, e2, factor);
2 26 Feb 07 jari 44     factor *= -1;
2 26 Feb 07 jari 45     break;
2 26 Feb 07 jari 46       case Algorithm.PEARSONSQARED:
2 26 Feb 07 jari 47     result = (float)Math.pow(pearsonUncentered(matrix, e1, e2, factor), 2);
2 26 Feb 07 jari 48     factor *= -1;
2 26 Feb 07 jari 49     break;
2 26 Feb 07 jari 50       case Algorithm.MANHATTAN:
2 26 Feb 07 jari 51     result = manhattan(matrix, e1, e2, factor);
2 26 Feb 07 jari 52     break;
2 26 Feb 07 jari 53       case Algorithm.SPEARMANRANK:
2 26 Feb 07 jari 54     result = spearmanRank(matrix, e1, e2, factor);
2 26 Feb 07 jari 55     factor *= -1;
2 26 Feb 07 jari 56     break;
2 26 Feb 07 jari 57       case Algorithm.KENDALLSTAU:
2 26 Feb 07 jari 58     result = kendallsTau(matrix, e1, e2, factor);
2 26 Feb 07 jari 59     factor *= -1;
2 26 Feb 07 jari 60     break;
2 26 Feb 07 jari 61       case Algorithm.MUTUALINFORMATION:
2 26 Feb 07 jari 62     result = mutualInformation(matrix, e1, e2, factor);
2 26 Feb 07 jari 63     break;
2 26 Feb 07 jari 64       default: {}
2 26 Feb 07 jari 65   }
2 26 Feb 07 jari 66   if (absolute) {
2 26 Feb 07 jari 67       result = Math.abs(result);
2 26 Feb 07 jari 68   }
2 26 Feb 07 jari 69   
2 26 Feb 07 jari 70   return result * factor;
2 26 Feb 07 jari 71     }
2 26 Feb 07 jari 72     
2 26 Feb 07 jari 73     public static float geneDistance(FloatMatrix matrix, FloatMatrix M, int g1, int g2, int function, float factor, boolean absolute) {
2 26 Feb 07 jari 74   float result = Float.NaN;
2 26 Feb 07 jari 75   switch (function) {
2 26 Feb 07 jari 76       case Algorithm.PEARSON:
2 26 Feb 07 jari 77     result = genePearson(matrix, M, g1, g2, factor);
2 26 Feb 07 jari 78     factor *= -1;
2 26 Feb 07 jari 79     break;
2 26 Feb 07 jari 80       case Algorithm.COSINE:
2 26 Feb 07 jari 81     result = geneCosine(matrix, M, g1, g2, factor);
2 26 Feb 07 jari 82     factor *= -1;
2 26 Feb 07 jari 83     break;
2 26 Feb 07 jari 84       case Algorithm.COVARIANCE:
2 26 Feb 07 jari 85     result = geneCovariance(matrix, M, g1, g2, factor);
2 26 Feb 07 jari 86     factor *= -1;
2 26 Feb 07 jari 87     break;
2 26 Feb 07 jari 88       case Algorithm.EUCLIDEAN:
2 26 Feb 07 jari 89     result = geneEuclidianDistance(matrix, M, g1, g2, factor);
2 26 Feb 07 jari 90     break;
2 26 Feb 07 jari 91       case Algorithm.DOTPRODUCT:
2 26 Feb 07 jari 92     result = geneDotProduct(matrix, M, g1, g2, factor);
2 26 Feb 07 jari 93     factor *= -1;
2 26 Feb 07 jari 94     break;
2 26 Feb 07 jari 95       case Algorithm.PEARSONUNCENTERED:
2 26 Feb 07 jari 96     result = genePearsonUncentered(matrix, M, g1, g2, factor);
2 26 Feb 07 jari 97     factor *= -1;
2 26 Feb 07 jari 98     break;
2 26 Feb 07 jari 99       case Algorithm.PEARSONSQARED:
2 26 Feb 07 jari 100     result = (float)Math.pow(genePearsonUncentered(matrix, M, g1, g2, factor), 2)*factor;
2 26 Feb 07 jari 101     factor *= -1;
2 26 Feb 07 jari 102     break;
2 26 Feb 07 jari 103       case Algorithm.MANHATTAN:
2 26 Feb 07 jari 104     result = geneManhattan(matrix, M, g1, g2, factor);
2 26 Feb 07 jari 105     break;
2 26 Feb 07 jari 106       case Algorithm.SPEARMANRANK:
2 26 Feb 07 jari 107     result = geneSpearmanRank(matrix, M, g1, g2, factor);
2 26 Feb 07 jari 108     factor *= -1;
2 26 Feb 07 jari 109     break;
2 26 Feb 07 jari 110       case Algorithm.KENDALLSTAU:
2 26 Feb 07 jari 111     result = geneKendallsTau(matrix, M, g1, g2, factor);
2 26 Feb 07 jari 112     factor *= -1;
2 26 Feb 07 jari 113     break;
2 26 Feb 07 jari 114       case Algorithm.MUTUALINFORMATION:
2 26 Feb 07 jari 115     result = geneMutualInformation(matrix, M, g1, g2, factor);
2 26 Feb 07 jari 116     break;
2 26 Feb 07 jari 117       default: {}
2 26 Feb 07 jari 118   }
2 26 Feb 07 jari 119   if (absolute) {
2 26 Feb 07 jari 120       result = Math.abs(result);
2 26 Feb 07 jari 121   }
2 26 Feb 07 jari 122   return result * factor;
2 26 Feb 07 jari 123     }
2 26 Feb 07 jari 124     
2 26 Feb 07 jari 125     //=================
2 26 Feb 07 jari 126     public static float genePearsonOld(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 127   if (M == null) {
2 26 Feb 07 jari 128       M = matrix;
2 26 Feb 07 jari 129   }
2 26 Feb 07 jari 130   float TINY = Float.MIN_VALUE;
2 26 Feb 07 jari 131   double xt,yt;
2 26 Feb 07 jari 132   double sxx=0.0;
2 26 Feb 07 jari 133   double syy=0.0;
2 26 Feb 07 jari 134   double sxy=0.0;
2 26 Feb 07 jari 135   double ax =0.0;
2 26 Feb 07 jari 136   double ay =0.0;
2 26 Feb 07 jari 137   int k = matrix.getColumnDimension();
2 26 Feb 07 jari 138   int n = 0;
2 26 Feb 07 jari 139   int j;
2 26 Feb 07 jari 140   for (j=0; j<k; j++) {
2 26 Feb 07 jari 141       if ((!Float.isNaN(matrix.get(g1,j))) && (!Float.isNaN(M.get(g2,j)))) {
2 26 Feb 07 jari 142     ax += matrix.get(g1,j);
2 26 Feb 07 jari 143     ay += M.get(g2,j);
2 26 Feb 07 jari 144     n++;
2 26 Feb 07 jari 145       }
2 26 Feb 07 jari 146   }
2 26 Feb 07 jari 147   ax /= n;
2 26 Feb 07 jari 148   ay /= n;
2 26 Feb 07 jari 149   for (j=0; j<k; j++) {
2 26 Feb 07 jari 150       if ((!Float.isNaN(matrix.get(g1,j))) && (!Float.isNaN(M.get(g2,j)))) {
2 26 Feb 07 jari 151     xt=matrix.get(g1,j)-ax;
2 26 Feb 07 jari 152     yt=M.get(g2,j)-ay;
2 26 Feb 07 jari 153     sxx+=xt*xt;
2 26 Feb 07 jari 154     syy+=yt*yt;
2 26 Feb 07 jari 155     sxy+=xt*yt;
2 26 Feb 07 jari 156       }
2 26 Feb 07 jari 157   }
2 26 Feb 07 jari 158   return(float)(sxy/(Math.sqrt(sxx*syy)+TINY)*factor);
2 26 Feb 07 jari 159     }
2 26 Feb 07 jari 160     
2 26 Feb 07 jari 161     //=================
2 26 Feb 07 jari 162     
2 26 Feb 07 jari 163     public static float genePearson(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 164   if (M == null) {
2 26 Feb 07 jari 165       M = matrix;
2 26 Feb 07 jari 166   }
2 26 Feb 07 jari 167   
2 26 Feb 07 jari 168   float[] arrX = matrix.A[g1];
2 26 Feb 07 jari 169   float[] arrY = M.A[g2];
2 26 Feb 07 jari 170   int nArrSize = matrix.getColumnDimension();
2 26 Feb 07 jari 171   
2 26 Feb 07 jari 172   double dblXY = 0f;
2 26 Feb 07 jari 173   double dblX  = 0f;
2 26 Feb 07 jari 174   double dblXX = 0f;
2 26 Feb 07 jari 175   double dblY  = 0f;
2 26 Feb 07 jari 176   double dblYY = 0f;
2 26 Feb 07 jari 177   
2 26 Feb 07 jari 178   double v_1, v_2;
2 26 Feb 07 jari 179   int iValidValCount = 0;
2 26 Feb 07 jari 180   for (int i=0; i<nArrSize; i++) {
2 26 Feb 07 jari 181       v_1 = arrX[i];
2 26 Feb 07 jari 182       v_2 = arrY[i];
2 26 Feb 07 jari 183       if (Double.isNaN(v_1) || Double.isNaN(v_2)) {
2 26 Feb 07 jari 184     continue;
2 26 Feb 07 jari 185       }
2 26 Feb 07 jari 186       iValidValCount++;
2 26 Feb 07 jari 187       dblXY += v_1*v_2;
2 26 Feb 07 jari 188       dblXX += v_1*v_1;
2 26 Feb 07 jari 189       dblYY += v_2*v_2;
2 26 Feb 07 jari 190       dblX  += v_1;
2 26 Feb 07 jari 191       dblY  += v_2;
2 26 Feb 07 jari 192   }
2 26 Feb 07 jari 193   if (iValidValCount == 0)
2 26 Feb 07 jari 194       return 0f;
2 26 Feb 07 jari 195   
2 26 Feb 07 jari 196   //Allows for a comparison of two 'flat' genes (genes with no variability in their
2 26 Feb 07 jari 197   // expression values), ie. 0, 0, 0, 0, 0
2 26 Feb 07 jari 198   boolean nonFlat = false;
2 26 Feb 07 jari 199   NON_FLAT_CHECK: for (int j = 1; j < nArrSize; j++) {
2 26 Feb 07 jari 200       if ((!Float.isNaN(arrX[j])) && (!Float.isNaN(arrY[j]))) {
2 26 Feb 07 jari 201     if (arrX[j] != arrX[j-1]) {
2 26 Feb 07 jari 202         nonFlat = true;
2 26 Feb 07 jari 203         break NON_FLAT_CHECK;
2 26 Feb 07 jari 204     }
2 26 Feb 07 jari 205     if (arrY[j] != arrY[j-1]) {
2 26 Feb 07 jari 206         nonFlat = true;
2 26 Feb 07 jari 207         break NON_FLAT_CHECK;
2 26 Feb 07 jari 208     }
2 26 Feb 07 jari 209       }
2 26 Feb 07 jari 210   }
2 26 Feb 07 jari 211   
2 26 Feb 07 jari 212   if (nonFlat == false) {
2 26 Feb 07 jari 213       return 1.0f;
2 26 Feb 07 jari 214   }
2 26 Feb 07 jari 215   
2 26 Feb 07 jari 216   
2 26 Feb 07 jari 217   double dblAvgX = dblX/iValidValCount;
2 26 Feb 07 jari 218   double dblAvgY = dblY/iValidValCount;
2 26 Feb 07 jari 219   double dblUpper = dblXY-dblX*dblAvgY-dblAvgX*dblY+dblAvgX*dblAvgY*((double)iValidValCount);
2 26 Feb 07 jari 220   double p1 = (dblXX-dblAvgX*dblX*2d+dblAvgX*dblAvgX*((double)iValidValCount));
2 26 Feb 07 jari 221   double p2 = (dblYY-dblAvgY*dblY*2d+dblAvgY*dblAvgY*((double)iValidValCount));
2 26 Feb 07 jari 222   double dblLower = p1*p2;
2 26 Feb 07 jari 223   return(float)(dblUpper/(Math.sqrt(dblLower)+Double.MIN_VALUE)*(double)factor);
2 26 Feb 07 jari 224     }
2 26 Feb 07 jari 225     
2 26 Feb 07 jari 226     public static float geneCosine(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 227   if (M == null) {
2 26 Feb 07 jari 228       M = matrix;
2 26 Feb 07 jari 229   }
2 26 Feb 07 jari 230   double xt,yt;
2 26 Feb 07 jari 231   double sxy=0.0;
2 26 Feb 07 jari 232   double sxx=0.0;
2 26 Feb 07 jari 233   double syy=0.0;
2 26 Feb 07 jari 234   double tx=0.0;
2 26 Feb 07 jari 235   double ty=0.0;
2 26 Feb 07 jari 236   int k = matrix.getColumnDimension();
2 26 Feb 07 jari 237   int n = 0;
2 26 Feb 07 jari 238   int j;
2 26 Feb 07 jari 239   for (j=0; j<k; j++) {
2 26 Feb 07 jari 240       if ((!Float.isNaN(matrix.get(g1, j))) && (!Float.isNaN(M.get(g2, j)))) {
2 26 Feb 07 jari 241     tx = matrix.get(g1, j);
2 26 Feb 07 jari 242     ty = M.get(g2, j);
2 26 Feb 07 jari 243     sxy += tx*ty;
2 26 Feb 07 jari 244     sxx += tx*tx;
2 26 Feb 07 jari 245     syy += ty*ty;
2 26 Feb 07 jari 246     n++;
2 26 Feb 07 jari 247       }
2 26 Feb 07 jari 248   }
2 26 Feb 07 jari 249   return(float)(sxy/(Math.sqrt(sxx)*Math.sqrt(syy))*factor);
2 26 Feb 07 jari 250     }
2 26 Feb 07 jari 251     
2 26 Feb 07 jari 252     public static float geneCovariance(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 253   if (M == null) {
2 26 Feb 07 jari 254       M=matrix;
2 26 Feb 07 jari 255   }
2 26 Feb 07 jari 256   double xt,yt;
2 26 Feb 07 jari 257   double sxy=0.0;
2 26 Feb 07 jari 258   double ax=0.0;
2 26 Feb 07 jari 259   double ay=0.0;
2 26 Feb 07 jari 260   int k = matrix.getColumnDimension();
2 26 Feb 07 jari 261   int n = 0;
2 26 Feb 07 jari 262   int j;
2 26 Feb 07 jari 263   for (j=0; j<k; j++) {
2 26 Feb 07 jari 264       if ((!Float.isNaN(matrix.get(g1, j))) && (!Float.isNaN(M.get(g2, j)))) {
2 26 Feb 07 jari 265     ax += matrix.get(g1, j);
2 26 Feb 07 jari 266     ay += M.get(g2, j);
2 26 Feb 07 jari 267     n++;
2 26 Feb 07 jari 268       }
2 26 Feb 07 jari 269   }
2 26 Feb 07 jari 270   ax /= n;
2 26 Feb 07 jari 271   ay /= n;
2 26 Feb 07 jari 272   for (j=0; j<k; j++) {
2 26 Feb 07 jari 273       if ((!Float.isNaN(matrix.get(g1, j))) && (!Float.isNaN(M.get(g2, j)))) {
2 26 Feb 07 jari 274     xt=matrix.get(g1, j)-ax;
2 26 Feb 07 jari 275     yt=M.get(g2, j)-ay;
2 26 Feb 07 jari 276     sxy+=xt*yt;
2 26 Feb 07 jari 277       }
2 26 Feb 07 jari 278   }
2 26 Feb 07 jari 279   return(float)(sxy/((n-1)*1.0)*factor);
2 26 Feb 07 jari 280     }
2 26 Feb 07 jari 281     
2 26 Feb 07 jari 282     public static float geneEuclidianDistance(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 283   if (M == null) {
2 26 Feb 07 jari 284       M = matrix;
2 26 Feb 07 jari 285   }
2 26 Feb 07 jari 286   int k = matrix.getColumnDimension();
2 26 Feb 07 jari 287   int n = 0;
2 26 Feb 07 jari 288   double sum = 0.0;
2 26 Feb 07 jari 289   for (int i=0; i<k; i++) {
2 26 Feb 07 jari 290       if ((!Float.isNaN(matrix.get(g1,i))) && (!Float.isNaN(M.get(g2,i)))) {
2 26 Feb 07 jari 291     sum+=Math.pow((matrix.get(g1,i)-M.get(g2,i)),2);
2 26 Feb 07 jari 292     n++;
2 26 Feb 07 jari 293       }
2 26 Feb 07 jari 294   }
2 26 Feb 07 jari 295   return(float)(Math.sqrt(sum)*factor);
2 26 Feb 07 jari 296     }
2 26 Feb 07 jari 297     
2 26 Feb 07 jari 298     public static float geneDotProduct(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 299   if (M == null) {
2 26 Feb 07 jari 300       M = matrix;
2 26 Feb 07 jari 301   }
2 26 Feb 07 jari 302   int k=matrix.getColumnDimension();
2 26 Feb 07 jari 303   int n=0;
2 26 Feb 07 jari 304   double sum=0.0;
2 26 Feb 07 jari 305   for (int i=0; i<k; i++) {
2 26 Feb 07 jari 306       if ((!Float.isNaN(matrix.get(g1,i))) && (!Float.isNaN(M.get(g2,i)))) {
2 26 Feb 07 jari 307     sum+=matrix.get(g1,i)*M.get(g2,i);
2 26 Feb 07 jari 308     n++;
2 26 Feb 07 jari 309       }
2 26 Feb 07 jari 310   }
2 26 Feb 07 jari 311   return(float)(sum/((double)n)*factor);
2 26 Feb 07 jari 312     }
2 26 Feb 07 jari 313     
2 26 Feb 07 jari 314     public static float genePearsonUncentered(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 315   if (M == null) {
2 26 Feb 07 jari 316       M = matrix;
2 26 Feb 07 jari 317   }
2 26 Feb 07 jari 318   float TINY = Float.MIN_VALUE;
2 26 Feb 07 jari 319   double xt,yt;
2 26 Feb 07 jari 320   double sxx=0.0;
2 26 Feb 07 jari 321   double syy=0.0;
2 26 Feb 07 jari 322   double sxy=0.0;
2 26 Feb 07 jari 323   double ax =0.0;
2 26 Feb 07 jari 324   double ay =0.0;
2 26 Feb 07 jari 325   int k = matrix.getColumnDimension();
2 26 Feb 07 jari 326   int n = 0;
2 26 Feb 07 jari 327   int j;
2 26 Feb 07 jari 328   for (j=0; j<k; j++) {
2 26 Feb 07 jari 329       if ((!Float.isNaN(matrix.get(g1,j))) && (!Float.isNaN(M.get(g2,j)))) {
2 26 Feb 07 jari 330     ax += matrix.get(g1,j);
2 26 Feb 07 jari 331     ay += M.get(g2,j);
2 26 Feb 07 jari 332     n++;
2 26 Feb 07 jari 333       }
2 26 Feb 07 jari 334   }
2 26 Feb 07 jari 335   ax /= n;
2 26 Feb 07 jari 336   ay /= n;
2 26 Feb 07 jari 337   for (j=0; j<k; j++) {
2 26 Feb 07 jari 338       if ((!Float.isNaN(matrix.get(g1,j))) && (!Float.isNaN(M.get(g2,j)))) {
2 26 Feb 07 jari 339     xt=matrix.get(g1,j);
2 26 Feb 07 jari 340     yt=M.get(g2,j);
2 26 Feb 07 jari 341     sxx+=xt*xt;
2 26 Feb 07 jari 342     syy+=yt*yt;
2 26 Feb 07 jari 343     sxy+=xt*yt;
2 26 Feb 07 jari 344       }
2 26 Feb 07 jari 345   }
2 26 Feb 07 jari 346   return(float)(sxy/(Math.sqrt(sxx*syy)+TINY)*factor);
2 26 Feb 07 jari 347     }
2 26 Feb 07 jari 348     
2 26 Feb 07 jari 349     public static float geneManhattan(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 350   if (M == null) {
2 26 Feb 07 jari 351       M = matrix;
2 26 Feb 07 jari 352   }
2 26 Feb 07 jari 353   int j;
2 26 Feb 07 jari 354   double sum = 0.0;
2 26 Feb 07 jari 355   int n = matrix.getColumnDimension();
2 26 Feb 07 jari 356   for (j=0; j<n; j++) {
2 26 Feb 07 jari 357       if ((!Float.isNaN(matrix.get(g1, j))) && (!Float.isNaN(M.get(g2, j)))) {
2 26 Feb 07 jari 358     sum += Math.abs(matrix.get(g1,j)-matrix.get(g2, j));
2 26 Feb 07 jari 359       }
2 26 Feb 07 jari 360   }
2 26 Feb 07 jari 361   return(float)(sum*factor);
2 26 Feb 07 jari 362     }
2 26 Feb 07 jari 363     
2 26 Feb 07 jari 364     public static float geneSpearmanRank(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 365   if (M == null) {
2 26 Feb 07 jari 366       M = matrix;
2 26 Feb 07 jari 367   }
2 26 Feb 07 jari 368   int j;
2 26 Feb 07 jari 369   double vard,t,sg,sf,fac,en3n,en,df,aved;
2 26 Feb 07 jari 370   double[] wksp1;
2 26 Feb 07 jari 371   double[] wksp2;
2 26 Feb 07 jari 372   double d;
2 26 Feb 07 jari 373   
2 26 Feb 07 jari 374   int n = matrix.getColumnDimension();
2 26 Feb 07 jari 375   wksp1=new double[n];
2 26 Feb 07 jari 376   wksp2=new double[n];
2 26 Feb 07 jari 377   for (j=0;j<n;j++) {
2 26 Feb 07 jari 378       wksp1[j]=matrix.get(g1,j);
2 26 Feb 07 jari 379       wksp2[j]=M.get(g2,j);
2 26 Feb 07 jari 380   }
2 26 Feb 07 jari 381   sort2(wksp1,wksp2); // Sort each of the data arrays, and convert the entries to ranks.
2 26 Feb 07 jari 382   sf=crank(wksp1);
2 26 Feb 07 jari 383   sort2(wksp2,wksp1);
2 26 Feb 07 jari 384   sg=crank(wksp2);
2 26 Feb 07 jari 385   d=0.0;
2 26 Feb 07 jari 386   for (j=0;j<n;j++)
2 26 Feb 07 jari 387       d += Math.pow((wksp1[j]-wksp2[j]),2); // Sum the squared diference of ranks.
2 26 Feb 07 jari 388   en=n;
2 26 Feb 07 jari 389   en3n=en*en*en-en;
2 26 Feb 07 jari 390   aved=en3n/6.0-(sf+sg)/12.0; // Expectation value of D,
2 26 Feb 07 jari 391   fac=(1.0-sf/en3n)*(1.0-sg/en3n);
2 26 Feb 07 jari 392   vard=((en-1.0)*en*en*Math.pow((en+1.0),2)/36.0)*fac; // and variance of D give
2 26 Feb 07 jari 393   
2 26 Feb 07 jari 394   return(float)((1.0-(6.0/en3n)*(d+(sf+sg)/12.0))/Math.sqrt(fac)*factor); // Rank correlation coe cient,
2 26 Feb 07 jari 395     }
2 26 Feb 07 jari 396     
2 26 Feb 07 jari 397     public static float geneKendallsTau(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 398   if (M == null) {
2 26 Feb 07 jari 399       M = matrix;
2 26 Feb 07 jari 400   }
2 26 Feb 07 jari 401   float TINY = Float.MIN_VALUE;
2 26 Feb 07 jari 402   int n = matrix.getColumnDimension();
2 26 Feb 07 jari 403   int n2 = 0;
2 26 Feb 07 jari 404   int n1 = 0;
2 26 Feb 07 jari 405   int is = 0;
2 26 Feb 07 jari 406   double svar,aa,a2,a1;
2 26 Feb 07 jari 407   for (int j=0; j<n-1; j++) {               //Loop over rst member of pair,
2 26 Feb 07 jari 408       for (int k=(j+1); k<n; k++) {         //and second member.
2 26 Feb 07 jari 409     a1=matrix.get(g1,j)-matrix.get(g1,k);
2 26 Feb 07 jari 410     a2=M.get(g2,j)-M.get(g2,k);
2 26 Feb 07 jari 411     aa=a1*a2;
2 26 Feb 07 jari 412     if (aa!=0.0) {                    // Neither array has a tie.
2 26 Feb 07 jari 413         ++n1;
2 26 Feb 07 jari 414         ++n2;
2 26 Feb 07 jari 415         if (aa > 0.0) ++is;
2 26 Feb 07 jari 416         else --is;
2 26 Feb 07 jari 417     } else {                          // One or both arrays have ties.
2 26 Feb 07 jari 418         if (a1!=0.0) ++n1;             // An \extra x" event.
2 26 Feb 07 jari 419         if (a2!=0.0) ++n2;             // An \extra y" event.
2 26 Feb 07 jari 420     }
2 26 Feb 07 jari 421       }
2 26 Feb 07 jari 422   }
2 26 Feb 07 jari 423   return(float)(is/(Math.sqrt((double) n1)*Math.sqrt((double)n2)+TINY)*factor);
2 26 Feb 07 jari 424     }
2 26 Feb 07 jari 425     
2 26 Feb 07 jari 426    
2 26 Feb 07 jari 427     public static float geneMutualInformation(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 428   if (M == null) {
2 26 Feb 07 jari 429       M = matrix;
2 26 Feb 07 jari 430   }
2 26 Feb 07 jari 431   int n=M.getColumnDimension();
2 26 Feb 07 jari 432   int NumberOfBins=(int)Math.floor(Math.log(n)/Math.log(2));
2 26 Feb 07 jari 433   int Values=0;
2 26 Feb 07 jari 434   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 435       if ((!Float.isNaN(matrix.get(g1, i))) && (!Float.isNaN(M.get(g2,i))))
2 26 Feb 07 jari 436     Values++;
2 26 Feb 07 jari 437   }
2 26 Feb 07 jari 438   // lexa: position of arguments are different from mutualInformation method ???
2 26 Feb 07 jari 439   FloatMatrix Gene1Array = new FloatMatrix(Values, 1);
2 26 Feb 07 jari 440   FloatMatrix Gene2Array = new FloatMatrix(Values, 1);
2 26 Feb 07 jari 441   int k=0;
2 26 Feb 07 jari 442   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 443       if ((!Float.isNaN(matrix.get(g1, i))) && (!Float.isNaN(matrix.get(g2,i)))) {
2 26 Feb 07 jari 444     //Gene1Array.set(0,k,matrix.get(g1,i));
2 26 Feb 07 jari 445     //Gene2Array.set(0,k,M.get(g2,i));
2 26 Feb 07 jari 446                 Gene1Array.set(k,0,matrix.get(g1,i));
2 26 Feb 07 jari 447     Gene2Array.set(k,0,M.get(g2,i));
2 26 Feb 07 jari 448                 
2 26 Feb 07 jari 449     k++;
2 26 Feb 07 jari 450       }
2 26 Feb 07 jari 451   }
2 26 Feb 07 jari 452   n=Values;
2 26 Feb 07 jari 453   makeDigitalExperiment(Gene1Array, 0);
2 26 Feb 07 jari 454   makeDigitalExperiment(Gene2Array, 0);
2 26 Feb 07 jari 455   double[] P1=new double[NumberOfBins];
2 26 Feb 07 jari 456   double[] P2=new double[NumberOfBins];
2 26 Feb 07 jari 457   double[][] P12=new double[NumberOfBins][NumberOfBins];
2 26 Feb 07 jari 458   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 459             P1[(int)Gene1Array.get(0,i)-1]++;
2 26 Feb 07 jari 460       P2[(int)Gene2Array.get(0,i)-1]++;
2 26 Feb 07 jari 461       P12[(int)Gene1Array.get(0,i)-1][(int)Gene2Array.get(0,i)-1]++;
2 26 Feb 07 jari 462   }
2 26 Feb 07 jari 463   for (int i=0; i<P1.length; i++) {
2 26 Feb 07 jari 464       P1[i]/=n;
2 26 Feb 07 jari 465       P2[i]/=n;
2 26 Feb 07 jari 466       for (int j=0; j<P1.length; j++) {
2 26 Feb 07 jari 467     P12[i][j]/=n;
2 26 Feb 07 jari 468       }
2 26 Feb 07 jari 469   }
2 26 Feb 07 jari 470   double H1=0;
2 26 Feb 07 jari 471   double H2=0;
2 26 Feb 07 jari 472   double H12=0;
2 26 Feb 07 jari 473   double MI=0;
2 26 Feb 07 jari 474   for (int i=0; i<P1.length; i++) {
2 26 Feb 07 jari 475       if (P1[i]!=0) H1 +=P1[i]*Math.log(P1[i])/Math.log(2);
2 26 Feb 07 jari 476       if (P2[i]!=0) H2 +=P2[i]*Math.log(P2[i])/Math.log(2);
2 26 Feb 07 jari 477       for (int j=0; j<P2.length; j++) {
2 26 Feb 07 jari 478     if (P12[i][j]!=0) H12+=P12[i][j]*Math.log(P12[i][j])/Math.log(2);
2 26 Feb 07 jari 479       }
2 26 Feb 07 jari 480   }
2 26 Feb 07 jari 481   H1=-H1;
2 26 Feb 07 jari 482   H2=-H2;
2 26 Feb 07 jari 483   H12=-H12;
2 26 Feb 07 jari 484   MI=(H1+H2-H12)/Math.max(H1,H2);
2 26 Feb 07 jari 485   return(float)((1-MI)*factor);
2 26 Feb 07 jari 486     }
2 26 Feb 07 jari 487     
2 26 Feb 07 jari 488     
2 26 Feb 07 jari 489     public static float pearson(FloatMatrix matrix, int e1, int e2, float factor) {
2 26 Feb 07 jari 490   float TINY = Float.MIN_VALUE;
2 26 Feb 07 jari 491   int n, j, k;
2 26 Feb 07 jari 492   double xt,yt;
2 26 Feb 07 jari 493   double sxx=0.0;
2 26 Feb 07 jari 494   double syy=0.0;
2 26 Feb 07 jari 495   double sxy=0.0;
2 26 Feb 07 jari 496   double ax =0.0;
2 26 Feb 07 jari 497   double ay =0.0;
2 26 Feb 07 jari 498   k=matrix.getRowDimension();
2 26 Feb 07 jari 499   n=0;
2 26 Feb 07 jari 500   for (j=0; j<k; j++) {
2 26 Feb 07 jari 501       if ((!Float.isNaN(matrix.get(j,e1))) && (!Float.isNaN(matrix.get(j,e2)))) {
2 26 Feb 07 jari 502     ax += matrix.get(j,e1);
2 26 Feb 07 jari 503     ay += matrix.get(j,e2);
2 26 Feb 07 jari 504     n++;
2 26 Feb 07 jari 505       }
2 26 Feb 07 jari 506   }
2 26 Feb 07 jari 507   ax /= n;
2 26 Feb 07 jari 508   ay /= n;
2 26 Feb 07 jari 509   for (j=0; j<k; j++) {
2 26 Feb 07 jari 510       if ((!Float.isNaN(matrix.get(j,e1))) && (!Float.isNaN(matrix.get(j,e2)))) {
2 26 Feb 07 jari 511     xt=matrix.get(j,e1)-ax;
2 26 Feb 07 jari 512     yt=matrix.get(j,e2)-ay;
2 26 Feb 07 jari 513     sxx+=xt*xt;
2 26 Feb 07 jari 514     syy+=yt*yt;
2 26 Feb 07 jari 515     sxy+=xt*yt;
2 26 Feb 07 jari 516       }
2 26 Feb 07 jari 517   }
2 26 Feb 07 jari 518   return(float)(sxy/(Math.sqrt(sxx*syy)+TINY)*factor);
2 26 Feb 07 jari 519     }
2 26 Feb 07 jari 520     
2 26 Feb 07 jari 521     public static float cosine(FloatMatrix matrix, int e1, int e2, float factor) {
2 26 Feb 07 jari 522   int n, j, k;
2 26 Feb 07 jari 523   double xt,yt;
2 26 Feb 07 jari 524   double sxy=0.0;
2 26 Feb 07 jari 525   double sxx=0.0;
2 26 Feb 07 jari 526   double syy=0.0;
2 26 Feb 07 jari 527   double tx=0.0;
2 26 Feb 07 jari 528   double ty=0.0;
2 26 Feb 07 jari 529   k=matrix.getRowDimension();
2 26 Feb 07 jari 530   for (j=0; j<k; j++) {
2 26 Feb 07 jari 531       if ((!Float.isNaN(matrix.get(j,e1))) && (!Float.isNaN(matrix.get(j,e2)))) {
2 26 Feb 07 jari 532     tx=matrix.get(j,e1);
2 26 Feb 07 jari 533     ty=matrix.get(j,e2);
2 26 Feb 07 jari 534     sxy+=tx*ty;
2 26 Feb 07 jari 535     sxx+=tx*tx;
2 26 Feb 07 jari 536     syy+=ty*ty;
2 26 Feb 07 jari 537       }
2 26 Feb 07 jari 538   }
2 26 Feb 07 jari 539   return(float)(sxy/(Math.sqrt(sxx)*Math.sqrt(syy))*factor);
2 26 Feb 07 jari 540     }
2 26 Feb 07 jari 541     
2 26 Feb 07 jari 542     public static float covariance(FloatMatrix matrix, int e1, int e2, float factor) {
2 26 Feb 07 jari 543   int n, j, k;
2 26 Feb 07 jari 544   double xt,yt;
2 26 Feb 07 jari 545   double sxy=0.0;
2 26 Feb 07 jari 546   double ax=0.0;
2 26 Feb 07 jari 547   double ay=0.0;
2 26 Feb 07 jari 548   k=matrix.getRowDimension();
2 26 Feb 07 jari 549   n=0;
2 26 Feb 07 jari 550   for (j=0; j<k; j++) {
2 26 Feb 07 jari 551       if ((!Float.isNaN(matrix.get(j,e1))) && (!Float.isNaN(matrix.get(j,e2)))) {
2 26 Feb 07 jari 552     ax += matrix.get(j,e1);
2 26 Feb 07 jari 553     ay += matrix.get(j,e2);
2 26 Feb 07 jari 554     n++;
2 26 Feb 07 jari 555       }
2 26 Feb 07 jari 556   }
2 26 Feb 07 jari 557   ax /= n;
2 26 Feb 07 jari 558   ay /= n;
2 26 Feb 07 jari 559   for (j=0; j<k; j++) {
2 26 Feb 07 jari 560       if ((!Float.isNaN(matrix.get(j,e1))) && (!Float.isNaN(matrix.get(j,e2)))) {
2 26 Feb 07 jari 561     xt=matrix.get(j,e1)-ax;
2 26 Feb 07 jari 562     yt=matrix.get(j,e2)-ay;
2 26 Feb 07 jari 563     sxy+=xt*yt;
2 26 Feb 07 jari 564       }
2 26 Feb 07 jari 565   }
2 26 Feb 07 jari 566   return(float)(sxy/((n-1)*1.0f)*factor);
2 26 Feb 07 jari 567     }
2 26 Feb 07 jari 568     
2 26 Feb 07 jari 569     public static float euclidian(FloatMatrix matrix, int e1, int e2, float factor) {
2 26 Feb 07 jari 570   int n=matrix.getRowDimension();
2 26 Feb 07 jari 571   double sum=0.0;
2 26 Feb 07 jari 572   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 573       if ((!Float.isNaN(matrix.get(i,e1))) && (!Float.isNaN(matrix.get(i,e2)))) {
2 26 Feb 07 jari 574     sum+=Math.pow((matrix.get(i,e1)-matrix.get(i,e2)),2);
2 26 Feb 07 jari 575       }
2 26 Feb 07 jari 576   }
2 26 Feb 07 jari 577   return(float)(Math.sqrt(sum)*factor);
2 26 Feb 07 jari 578     }
2 26 Feb 07 jari 579     
2 26 Feb 07 jari 580     public static float dotProduct(FloatMatrix matrix, int e1, int e2, float factor) {
2 26 Feb 07 jari 581   int k=matrix.getRowDimension();
2 26 Feb 07 jari 582   double sum=0.0;
2 26 Feb 07 jari 583   int n=0;
2 26 Feb 07 jari 584   for (int i=0; i<k; i++) {
2 26 Feb 07 jari 585       if ((!Float.isNaN(matrix.get(i,e1))) && (!Float.isNaN(matrix.get(i,e2)))) {
2 26 Feb 07 jari 586     sum+=matrix.get(i,e1)*matrix.get(i,e2);
2 26 Feb 07 jari 587     n++;
2 26 Feb 07 jari 588       }
2 26 Feb 07 jari 589   }
2 26 Feb 07 jari 590   return(float)(sum/((double)n)*factor);
2 26 Feb 07 jari 591     }
2 26 Feb 07 jari 592     
2 26 Feb 07 jari 593     public static float pearsonUncentered(FloatMatrix matrix, int e1, int e2, float factor) {
2 26 Feb 07 jari 594   float TINY = Float.MIN_VALUE;
2 26 Feb 07 jari 595   int n, j, k;
2 26 Feb 07 jari 596   double xt,yt;
2 26 Feb 07 jari 597   double sxx=0.0;
2 26 Feb 07 jari 598   double syy=0.0;
2 26 Feb 07 jari 599   double sxy=0.0;
2 26 Feb 07 jari 600   double ax =0.0;
2 26 Feb 07 jari 601   double ay =0.0;
2 26 Feb 07 jari 602   k=matrix.getRowDimension();
2 26 Feb 07 jari 603   n=0;
2 26 Feb 07 jari 604   for (j=0; j<k; j++) {
2 26 Feb 07 jari 605       if ((!Float.isNaN(matrix.get(j,e1))) && (!Float.isNaN(matrix.get(j,e2)))) {
2 26 Feb 07 jari 606     ax += matrix.get(j,e1);
2 26 Feb 07 jari 607     ay += matrix.get(j,e2);
2 26 Feb 07 jari 608     n++;
2 26 Feb 07 jari 609       }
2 26 Feb 07 jari 610   }
2 26 Feb 07 jari 611   ax /= n;
2 26 Feb 07 jari 612   ay /= n;
2 26 Feb 07 jari 613   for (j=0; j<k; j++) {
2 26 Feb 07 jari 614       if ((!Float.isNaN(matrix.get(j,e1))) && (!Float.isNaN(matrix.get(j,e2)))) {
2 26 Feb 07 jari 615     xt=matrix.get(j,e1);
2 26 Feb 07 jari 616     yt=matrix.get(j,e2);
2 26 Feb 07 jari 617     sxx+=xt*xt;
2 26 Feb 07 jari 618     syy+=yt*yt;
2 26 Feb 07 jari 619     sxy+=xt*yt;
2 26 Feb 07 jari 620       }
2 26 Feb 07 jari 621   }
2 26 Feb 07 jari 622   return(float)(sxy/(Math.sqrt(sxx*syy)+TINY)*factor);
2 26 Feb 07 jari 623     }
2 26 Feb 07 jari 624     
2 26 Feb 07 jari 625     public static float manhattan(FloatMatrix matrix, int e1, int e2, float factor) {
2 26 Feb 07 jari 626   int n, j;
2 26 Feb 07 jari 627   double sum=0.0;
2 26 Feb 07 jari 628   n=matrix.getRowDimension();
2 26 Feb 07 jari 629   for (j=0; j<n; j++) {
2 26 Feb 07 jari 630       if ((!Float.isNaN(matrix.get(j,e1))) && (!Float.isNaN(matrix.get(j,e2)))) {
2 26 Feb 07 jari 631     sum+=Math.abs(matrix.get(j,e1)-matrix.get(j,e2));
2 26 Feb 07 jari 632       }
2 26 Feb 07 jari 633   }
2 26 Feb 07 jari 634   return(float)(sum*factor);
2 26 Feb 07 jari 635     }
2 26 Feb 07 jari 636     
2 26 Feb 07 jari 637     public static float spearmanRank(FloatMatrix matrix, int e1, int e2, float factor) {
2 26 Feb 07 jari 638   int j;
2 26 Feb 07 jari 639   double vard,t,sg,sf,fac,en3n,en,df,aved;
2 26 Feb 07 jari 640   double[] wksp1;
2 26 Feb 07 jari 641   double[] wksp2;
2 26 Feb 07 jari 642   double d;
2 26 Feb 07 jari 643   
2 26 Feb 07 jari 644   int n=matrix.getRowDimension();
2 26 Feb 07 jari 645   wksp1=new double[n];
2 26 Feb 07 jari 646   wksp2=new double[n];
2 26 Feb 07 jari 647   for (j=0;j<n;j++) {
2 26 Feb 07 jari 648       wksp1[j]=matrix.get(j,e1);
2 26 Feb 07 jari 649       wksp2[j]=matrix.get(j,e2);
2 26 Feb 07 jari 650   }
2 26 Feb 07 jari 651   sort2(wksp1,wksp2); // Sort each of the data arrays, and convert the entries to ranks.
2 26 Feb 07 jari 652   sf=crank(wksp1);
2 26 Feb 07 jari 653   sort2(wksp2,wksp1);
2 26 Feb 07 jari 654   sg=crank(wksp2);
2 26 Feb 07 jari 655   d=0.0;
2 26 Feb 07 jari 656   for (j=0;j<n;j++) d += Math.pow((wksp1[j]-wksp2[j]),2); // Sum the squared diference of ranks.
2 26 Feb 07 jari 657   en=n;
2 26 Feb 07 jari 658   en3n=en*en*en-en;
2 26 Feb 07 jari 659   aved=en3n/6.0-(sf+sg)/12.0; // Expectation value of D,
2 26 Feb 07 jari 660   fac=(1.0-sf/en3n)*(1.0-sg/en3n);
2 26 Feb 07 jari 661   vard=((en-1.0)*en*en*Math.pow((en+1.0),2)/36.0)*fac; // and variance of D give
2 26 Feb 07 jari 662   
2 26 Feb 07 jari 663   return(float)((1.0-(6.0/en3n)*(d+(sf+sg)/12.0))/Math.sqrt(fac)*factor); // Rank correlation coecient,
2 26 Feb 07 jari 664     }
2 26 Feb 07 jari 665     
2 26 Feb 07 jari 666     public static float kendallsTau(FloatMatrix matrix, int e1, int e2, float factor) {
2 26 Feb 07 jari 667   float TINY=Float.MIN_VALUE;
2 26 Feb 07 jari 668   int n=matrix.getRowDimension();
2 26 Feb 07 jari 669   int n2=0;
2 26 Feb 07 jari 670   int n1=0;
2 26 Feb 07 jari 671   int is=0;
2 26 Feb 07 jari 672   double svar,aa,a2,a1;
2 26 Feb 07 jari 673   for (int j=0; j<n-1; j++) {               //Loop over rst member of pair,
2 26 Feb 07 jari 674       for (int k=(j+1); k<n; k++) {         //and second member.
2 26 Feb 07 jari 675     a1=matrix.get(j,e1)-matrix.get(k,e1);
2 26 Feb 07 jari 676     a2=matrix.get(j,e2)-matrix.get(k,e2);
2 26 Feb 07 jari 677     aa=a1*a2;
2 26 Feb 07 jari 678     if (aa!=0.0) {                    // Neither array has a tie.
2 26 Feb 07 jari 679         ++n1;
2 26 Feb 07 jari 680         ++n2;
2 26 Feb 07 jari 681         if (aa > 0.0) ++is;
2 26 Feb 07 jari 682         else --is;
2 26 Feb 07 jari 683     } else {                          // One or both arrays have ties.
2 26 Feb 07 jari 684         if (a1!=0.0) ++n1;             // An \extra x" event.
2 26 Feb 07 jari 685         if (a2!=0.0) ++n2;             // An \extra y" event.
2 26 Feb 07 jari 686     }
2 26 Feb 07 jari 687       }
2 26 Feb 07 jari 688   }
2 26 Feb 07 jari 689   return(float)(is/(Math.sqrt((double) n1)*Math.sqrt((double)n2)+TINY)*factor);
2 26 Feb 07 jari 690     }
2 26 Feb 07 jari 691     
2 26 Feb 07 jari 692     public static float mutualInformation(FloatMatrix matrix, int e1, int e2, float factor) {
2 26 Feb 07 jari 693   int n=matrix.getRowDimension();
2 26 Feb 07 jari 694   int NumberOfBins=(int)Math.floor(Math.log(n)/Math.log(2));
2 26 Feb 07 jari 695   int Values=0;
2 26 Feb 07 jari 696   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 697       if ((!Float.isNaN(matrix.get(i,e1))) && (!Float.isNaN(matrix.get(i,e2)))) Values++;
2 26 Feb 07 jari 698   }
2 26 Feb 07 jari 699   
2 26 Feb 07 jari 700   FloatMatrix Experiment1Array = new FloatMatrix(1, Values);
2 26 Feb 07 jari 701   FloatMatrix Experiment2Array = new FloatMatrix(1, Values);
2 26 Feb 07 jari 702   int k=0;
2 26 Feb 07 jari 703   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 704       if ((!Float.isNaN(matrix.get(i,e1))) && (!Float.isNaN(matrix.get(i,e2)))) {
2 26 Feb 07 jari 705     Experiment1Array.set(k,0,matrix.get(i,e1));
2 26 Feb 07 jari 706     Experiment2Array.set(k,0,matrix.get(i,e2));
2 26 Feb 07 jari 707     k++;
2 26 Feb 07 jari 708       }
2 26 Feb 07 jari 709   }
2 26 Feb 07 jari 710   n=Values;
2 26 Feb 07 jari 711   makeDigitalExperiment(Experiment1Array, 0);
2 26 Feb 07 jari 712   makeDigitalExperiment(Experiment2Array, 0);
2 26 Feb 07 jari 713   double[] P1=new double[NumberOfBins];
2 26 Feb 07 jari 714   double[] P2=new double[NumberOfBins];
2 26 Feb 07 jari 715   double[][] P12=new double[NumberOfBins][NumberOfBins];
2 26 Feb 07 jari 716   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 717       P1[(int)Experiment1Array.get(i,0)-1]++;
2 26 Feb 07 jari 718       P2[(int)Experiment2Array.get(i,0)-1]++;
2 26 Feb 07 jari 719       P12[(int)Experiment1Array.get(i,0)-1][(int)Experiment2Array.get(i,0)-1]++;
2 26 Feb 07 jari 720   }
2 26 Feb 07 jari 721   for (int i=0; i<P1.length; i++) {
2 26 Feb 07 jari 722       P1[i]/=n;
2 26 Feb 07 jari 723       P2[i]/=n;
2 26 Feb 07 jari 724       for (int j=0; j<P1.length; j++) {
2 26 Feb 07 jari 725     P12[i][j]/=n;
2 26 Feb 07 jari 726       }
2 26 Feb 07 jari 727   }
2 26 Feb 07 jari 728   double H1=0;
2 26 Feb 07 jari 729   double H2=0;
2 26 Feb 07 jari 730   double H12=0;
2 26 Feb 07 jari 731   double MI=0;
2 26 Feb 07 jari 732   for (int i=0; i<P1.length; i++) {
2 26 Feb 07 jari 733       if (P1[i]!=0) H1 +=P1[i]*Math.log(P1[i])/Math.log(2);
2 26 Feb 07 jari 734       if (P2[i]!=0) H2 +=P2[i]*Math.log(P2[i])/Math.log(2);
2 26 Feb 07 jari 735       for (int j=0; j<P2.length; j++) {
2 26 Feb 07 jari 736     if (P12[i][j]!=0) H12+=P12[i][j]*Math.log(P12[i][j])/Math.log(2);
2 26 Feb 07 jari 737       }
2 26 Feb 07 jari 738   }
2 26 Feb 07 jari 739   H1=-H1;
2 26 Feb 07 jari 740   H2=-H2;
2 26 Feb 07 jari 741   H12=-H12;
2 26 Feb 07 jari 742   if (Math.max(H1,H2)!=0) {
2 26 Feb 07 jari 743       MI=(H1+H2-H12)/Math.max(H1,H2);
2 26 Feb 07 jari 744   } else {
2 26 Feb 07 jari 745       MI=(H1+H2-H12);
2 26 Feb 07 jari 746   }
2 26 Feb 07 jari 747   return(float)((1-MI)*factor);
2 26 Feb 07 jari 748     }
2 26 Feb 07 jari 749     
2 26 Feb 07 jari 750     public static void makeDigitalExperiment(FloatMatrix matrix, int e) {
2 26 Feb 07 jari 751   int n=matrix.getRowDimension();
2 26 Feb 07 jari 752   int NumberOfBins=(int)Math.floor(Math.log(n)/Math.log(2));
2 26 Feb 07 jari 753   int Step=1000000/NumberOfBins;
2 26 Feb 07 jari 754   float Minimum=Float.MAX_VALUE;
2 26 Feb 07 jari 755   float Maximum=0;
2 26 Feb 07 jari 756   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 757       if (matrix.get(i,e)<Minimum) Minimum=matrix.get(i,e);
2 26 Feb 07 jari 758   }
2 26 Feb 07 jari 759   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 760       matrix.set(i,e,matrix.get(i,e)-Minimum);
2 26 Feb 07 jari 761   }
2 26 Feb 07 jari 762   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 763       if (matrix.get(i,e)>Maximum) Maximum=matrix.get(i,e);
2 26 Feb 07 jari 764   }
2 26 Feb 07 jari 765   if (Maximum!=0) {
2 26 Feb 07 jari 766       for (int i=0; i<n; i++) {
2 26 Feb 07 jari 767     matrix.set(i,e,matrix.get(i,e)/Maximum);
2 26 Feb 07 jari 768       }
2 26 Feb 07 jari 769   }
2 26 Feb 07 jari 770   for (int i=0; i<n; i++) {
2 26 Feb 07 jari 771       if (matrix.get(i,e)==1.0) {
2 26 Feb 07 jari 772     matrix.set(i,e,(float)NumberOfBins);
2 26 Feb 07 jari 773       } else {
2 26 Feb 07 jari 774     matrix.set(i,e,(float)(Math.floor(matrix.get(i,e)*1000000/Step)+1));
2 26 Feb 07 jari 775       }
2 26 Feb 07 jari 776   }
2 26 Feb 07 jari 777     }
2 26 Feb 07 jari 778     
2 26 Feb 07 jari 779     /** Given a sorted array w[0..n-1], replaces the elements by their rank,
2 26 Feb 07 jari 780      *  including midranking of ties, and returns as s the sum of f3 ? f,
2 26 Feb 07 jari 781      *  where fis the number of elements in each tie.
2 26 Feb 07 jari 782      */
2 26 Feb 07 jari 783     public static double crank(double w[]) {
2 26 Feb 07 jari 784   int j=0,ji,jt;
2 26 Feb 07 jari 785   double t,rank;
2 26 Feb 07 jari 786   double s=0.0;
2 26 Feb 07 jari 787   int n=w.length;
2 26 Feb 07 jari 788   while (j < n-1) {
2 26 Feb 07 jari 789       if (w[j+1] != w[j]) { // Not a tie.
2 26 Feb 07 jari 790     w[j]=j;
2 26 Feb 07 jari 791     ++j;
2 26 Feb 07 jari 792       } else { // A tie:
2 26 Feb 07 jari 793     for (jt=j+1; jt<n && w[jt]==w[j]; jt++); // How far does it go?
2 26 Feb 07 jari 794     rank=0.5*(j+jt-1);                       // This is the mean rank of the tie,
2 26 Feb 07 jari 795     for (ji=j;ji<=(jt-1);ji++) w[ji]=rank;   // so enter it into all the tied entries,
2 26 Feb 07 jari 796     t=jt-j;
2 26 Feb 07 jari 797     s += t*t*t-t;                            // and update s.
2 26 Feb 07 jari 798     j=jt;
2 26 Feb 07 jari 799       }
2 26 Feb 07 jari 800   }
2 26 Feb 07 jari 801   if (j == n-1) w[n-1]=n-1; // If the last element was not tied, this is its rank.
2 26 Feb 07 jari 802   return s;
2 26 Feb 07 jari 803     }
2 26 Feb 07 jari 804     
2 26 Feb 07 jari 805     /**
2 26 Feb 07 jari 806      * Sorts an array arr[1..n] into ascending order using Quicksort,
2 26 Feb 07 jari 807      * while making the corresponding rearrangement of the array brr[1..n].
2 26 Feb 07 jari 808      */
2 26 Feb 07 jari 809     public static void sort2(double[] arr, double[] brr) {
2 26 Feb 07 jari 810   int n=arr.length;
2 26 Feb 07 jari 811   int i,ir=n-1,j,k,l=0;
2 26 Feb 07 jari 812   int[] istack;
2 26 Feb 07 jari 813   int jstack=0;
2 26 Feb 07 jari 814   double a,b,temp;
2 26 Feb 07 jari 815   double dummy;
2 26 Feb 07 jari 816   istack=new int[50];
2 26 Feb 07 jari 817   for (;;) {
2 26 Feb 07 jari 818       if (ir-l < 7) {
2 26 Feb 07 jari 819     for (j=l+1;j<=ir;j++) {
2 26 Feb 07 jari 820         a=arr[j];
2 26 Feb 07 jari 821         b=brr[j];
2 26 Feb 07 jari 822         for (i=j-1;i>=l;i--) {
2 26 Feb 07 jari 823       if (arr[i] <= a) break;
2 26 Feb 07 jari 824       arr[i+1]=arr[i];
2 26 Feb 07 jari 825       brr[i+1]=brr[i];
2 26 Feb 07 jari 826         }
2 26 Feb 07 jari 827         arr[i+1]=a;
2 26 Feb 07 jari 828         brr[i+1]=b;
2 26 Feb 07 jari 829     }
2 26 Feb 07 jari 830     if (jstack==0) {
2 26 Feb 07 jari 831         istack=null;
2 26 Feb 07 jari 832         return;
2 26 Feb 07 jari 833     }
2 26 Feb 07 jari 834     ir=istack[jstack];
2 26 Feb 07 jari 835     l=istack[jstack-1];
2 26 Feb 07 jari 836     jstack -= 2;
2 26 Feb 07 jari 837       } else {
2 26 Feb 07 jari 838     k=(l+ir) >> 1;
2 26 Feb 07 jari 839     dummy=arr[k];
2 26 Feb 07 jari 840     arr[k]=arr[l+1];
2 26 Feb 07 jari 841     arr[l+1]=arr[k];
2 26 Feb 07 jari 842     
2 26 Feb 07 jari 843     dummy=brr[k];
2 26 Feb 07 jari 844     brr[k]=brr[l+1];
2 26 Feb 07 jari 845     brr[l+1]=brr[k];
2 26 Feb 07 jari 846     
2 26 Feb 07 jari 847     if (arr[l] > arr[ir]) {
2 26 Feb 07 jari 848         dummy=arr[l];
2 26 Feb 07 jari 849         arr[l]=arr[ir];
2 26 Feb 07 jari 850         arr[ir]=dummy;
2 26 Feb 07 jari 851         
2 26 Feb 07 jari 852         dummy=brr[l];
2 26 Feb 07 jari 853         brr[l]=brr[ir];
2 26 Feb 07 jari 854         brr[ir]=dummy;
2 26 Feb 07 jari 855         
2 26 Feb 07 jari 856     }
2 26 Feb 07 jari 857     if (arr[l+1] > arr[ir]) {
2 26 Feb 07 jari 858         dummy=arr[l+1];
2 26 Feb 07 jari 859         arr[l+1]=arr[ir];
2 26 Feb 07 jari 860         arr[ir]=dummy;
2 26 Feb 07 jari 861         
2 26 Feb 07 jari 862         dummy=brr[l+1];
2 26 Feb 07 jari 863         brr[l+1]=brr[ir];
2 26 Feb 07 jari 864         brr[ir]=dummy;
2 26 Feb 07 jari 865     }
2 26 Feb 07 jari 866     if (arr[l] > arr[l+1]) {
2 26 Feb 07 jari 867         dummy=arr[l];
2 26 Feb 07 jari 868         arr[l]=arr[l+1];
2 26 Feb 07 jari 869         arr[l+1]=dummy;
2 26 Feb 07 jari 870         
2 26 Feb 07 jari 871         dummy=brr[l];
2 26 Feb 07 jari 872         brr[l]=brr[l+1];
2 26 Feb 07 jari 873         brr[l+1]=dummy;
2 26 Feb 07 jari 874     }
2 26 Feb 07 jari 875     i=l+1; // Initialize pointers for partitioning.
2 26 Feb 07 jari 876     j=ir;
2 26 Feb 07 jari 877     a=arr[l+1]; //Partitioning element.
2 26 Feb 07 jari 878     b=brr[l+1];
2 26 Feb 07 jari 879     for (;;) { // Beginning of innermost loop.
2 26 Feb 07 jari 880         do i++; while (arr[i] < a);
2 26 Feb 07 jari 881         do j--; while (arr[j] > a);
2 26 Feb 07 jari 882         if (j < i) break; // Pointers crossed. Partitioning complete.
2 26 Feb 07 jari 883         dummy=arr[i];
2 26 Feb 07 jari 884         arr[i]=arr[j];
2 26 Feb 07 jari 885         arr[j]=dummy;
2 26 Feb 07 jari 886         dummy=brr[i];
2 26 Feb 07 jari 887         brr[i]=brr[j];
2 26 Feb 07 jari 888         brr[j]=dummy;
2 26 Feb 07 jari 889     } // End of innermost loop.
2 26 Feb 07 jari 890     arr[l+1]=arr[j]; // Insert partitioning element in both arrays.
2 26 Feb 07 jari 891     arr[j]=a;
2 26 Feb 07 jari 892     brr[l+1]=brr[j];
2 26 Feb 07 jari 893     brr[j]=b;
2 26 Feb 07 jari 894     jstack += 2;
2 26 Feb 07 jari 895     if (jstack > 50) System.out.println("NSTACK too small in sort2.");
2 26 Feb 07 jari 896     if (ir-i+1 >= j-l) {
2 26 Feb 07 jari 897         istack[jstack]=ir;
2 26 Feb 07 jari 898         istack[jstack-1]=i;
2 26 Feb 07 jari 899         ir=j-1;
2 26 Feb 07 jari 900     } else {
2 26 Feb 07 jari 901         istack[jstack]=j-1;
2 26 Feb 07 jari 902         istack[jstack-1]=l;
2 26 Feb 07 jari 903         l=i;
2 26 Feb 07 jari 904     }
2 26 Feb 07 jari 905       }
2 26 Feb 07 jari 906   }
2 26 Feb 07 jari 907     }
2 26 Feb 07 jari 908     
2 26 Feb 07 jari 909
2 26 Feb 07 jari 910     
2 26 Feb 07 jari 911     /**
2 26 Feb 07 jari 912     public static float geneMutualInformation(FloatMatrix matrix, FloatMatrix M, int g1, int g2, float factor) {
2 26 Feb 07 jari 913         if (M == null) {
2 26 Feb 07 jari 914             M = matrix;
2 26 Feb 07 jari 915         }
2 26 Feb 07 jari 916         
2 26 Feb 07 jari 917         //System.out.println("Mutual info: Current genes: gene " + g1 + " gene " + g2);
2 26 Feb 07 jari 918         
2 26 Feb 07 jari 919         int numExps = matrix.getColumnDimension();
2 26 Feb 07 jari 920         float[] gene1 = new float[numExps];
2 26 Feb 07 jari 921         float[] gene2 = new float[numExps];
2 26 Feb 07 jari 922         
2 26 Feb 07 jari 923         for (int i = 0 ; i < numExps; i++) {
2 26 Feb 07 jari 924             gene1[i] = matrix.A[g1][i];
2 26 Feb 07 jari 925             gene2[i] = matrix.A[g2][i];
2 26 Feb 07 jari 926         }
2 26 Feb 07 jari 927         
2 26 Feb 07 jari 928         //Vector naNValuedElements = new Vector();
2 26 Feb 07 jari 929         boolean[] naNFlags = new boolean[numExps];
2 26 Feb 07 jari 930         int naNCounter = 0;
2 26 Feb 07 jari 931         
2 26 Feb 07 jari 932         for (int i = 0; i < numExps; i++) {
2 26 Feb 07 jari 933             naNFlags[i] = false;
2 26 Feb 07 jari 934         }
2 26 Feb 07 jari 935         
2 26 Feb 07 jari 936         for (int i = 0; i < numExps; i++) {
2 26 Feb 07 jari 937             if ((Float.isNaN(gene1[i]))||(Float.isNaN(gene2[i]))) {
2 26 Feb 07 jari 938                 //naNValuedElements.add(new Integer(i));
2 26 Feb 07 jari 939                 naNFlags[i] = true;
2 26 Feb 07 jari 940                 naNCounter++;
2 26 Feb 07 jari 941             }
2 26 Feb 07 jari 942         }
2 26 Feb 07 jari 943         
2 26 Feb 07 jari 944         int reducedGeneLength =  numExps - naNCounter; 
2 26 Feb 07 jari 945         int reducedGeneIndex = 0;
2 26 Feb 07 jari 946         float[] reducedGene1 = new float[numExps - naNCounter];
2 26 Feb 07 jari 947         float[] reducedGene2 = new float[numExps - naNCounter]; 
2 26 Feb 07 jari 948         
2 26 Feb 07 jari 949         for (int i = 0; i < numExps; i++) {
2 26 Feb 07 jari 950             if (naNFlags[i] == false) {
2 26 Feb 07 jari 951                 reducedGene1[reducedGeneIndex] = gene1[i];
2 26 Feb 07 jari 952                 reducedGene2[reducedGeneIndex] = gene2[i];
2 26 Feb 07 jari 953                 reducedGeneIndex++;
2 26 Feb 07 jari 954             }
2 26 Feb 07 jari 955         }
2 26 Feb 07 jari 956         //System.out.println("Entropy of gene " + g1);
2 26 Feb 07 jari 957         double h1 = calculateSingleEntropy(reducedGene1);
2 26 Feb 07 jari 958         //System.out.println();
2 26 Feb 07 jari 959         //System.out.println("Entropy of gene " + g2);
2 26 Feb 07 jari 960         double h2 = calculateSingleEntropy(reducedGene2);
2 26 Feb 07 jari 961         //System.out.println();
2 26 Feb 07 jari 962         double h12 = calculateJointEntropy(reducedGene1, reducedGene2);
2 26 Feb 07 jari 963         //System.out.println("Entropy of gene " + g1 + " = " + h1);
2 26 Feb 07 jari 964         //System.out.println("Entropy of gene " + g2 + " = " + h2);
2 26 Feb 07 jari 965         //System.out.println("Joint entropy of gene " + g1 + " and gene " + g2 + " = " + h12);
2 26 Feb 07 jari 966         
2 26 Feb 07 jari 967         
2 26 Feb 07 jari 968         //double mutualInfo=(h1+h2-h12)/Math.max(h1,h2);
2 26 Feb 07 jari 969         //float miDistance = (float)((1-mutualInfo)*factor);
2 26 Feb 07 jari 970          
2 26 Feb 07 jari 971         double mutualInfo = h1 + h2 - h12;
2 26 Feb 07 jari 972         float miDistance = (float)(h12 - mutualInfo)*factor;
2 26 Feb 07 jari 973         //System.out.println("Mutual info distance between gene " + g1 + " and gene " + g2 + " = " + miDistance);
2 26 Feb 07 jari 974         return miDistance;
2 26 Feb 07 jari 975         
2 26 Feb 07 jari 976         //float max1 = getMax(gene1);
2 26 Feb 07 jari 977         //loat min1 = getMin(gene1);
2 26 Feb 07 jari 978         
2 26 Feb 07 jari 979         //float max2 = getMax(gene2);
2 26 Feb 07 jari 980         //float min2 = getMin(gene2);          
2 26 Feb 07 jari 981     }
2 26 Feb 07 jari 982
2 26 Feb 07 jari 983     private static float getMax(float[] gene) {
2 26 Feb 07 jari 984         float max = Float.NEGATIVE_INFINITY;
2 26 Feb 07 jari 985         
2 26 Feb 07 jari 986         for (int i = 0; i < gene.length; i++) {
2 26 Feb 07 jari 987             if (gene[i] > max) {
2 26 Feb 07 jari 988                 max = gene[i];
2 26 Feb 07 jari 989             }
2 26 Feb 07 jari 990         }
2 26 Feb 07 jari 991         
2 26 Feb 07 jari 992         return max;
2 26 Feb 07 jari 993     }
2 26 Feb 07 jari 994
2 26 Feb 07 jari 995     private static float getMin(float[] gene) {
2 26 Feb 07 jari 996         float min = Float.POSITIVE_INFINITY;
2 26 Feb 07 jari 997         
2 26 Feb 07 jari 998         for (int i = 0; i < gene.length; i++) {
2 26 Feb 07 jari 999             if (gene[i] < min) {
2 26 Feb 07 jari 1000                 min = gene[i];
2 26 Feb 07 jari 1001             }
2 26 Feb 07 jari 1002         }
2 26 Feb 07 jari 1003         
2 26 Feb 07 jari 1004         return min;
2 26 Feb 07 jari 1005     }  
2 26 Feb 07 jari 1006     
2 26 Feb 07 jari 1007     
2 26 Feb 07 jari 1008     private static Vector[] makeBins(float[] gene) {
2 26 Feb 07 jari 1009   
2 26 Feb 07 jari 1010         
2 26 Feb 07 jari 1011         if (getMax(gene) == getMin(gene)) {
2 26 Feb 07 jari 1012             for (int i = 0; i < gene.length; i++) {//this is needed to introduce some variability into "flat" genes, otherwise entropy calculations cause subsequent problems
2 26 Feb 07 jari 1013                 gene[i] = 0;
2 26 Feb 07 jari 1014             }
2 26 Feb 07 jari 1015             gene[0] = gene[0] + Float.MIN_VALUE;
2 26 Feb 07 jari 1016             
2 26 Feb 07 jari 1017         }
2 26 Feb 07 jari 1018         
2 26 Feb 07 jari 1019         float max = getMax(gene);
2 26 Feb 07 jari 1020         float min = getMin(gene);
2 26 Feb 07 jari 1021         
2 26 Feb 07 jari 1022         //System.out.println("Max = " + max);
2 26 Feb 07 jari 1023         //System.out.println("Min = " + min);
2 26 Feb 07 jari 1024         
2 26 Feb 07 jari 1025         Vector[] bins = new Vector[10];
2 26 Feb 07 jari 1026         for (int i = 0; i < bins.length; i++) {
2 26 Feb 07 jari 1027             bins[i] = new Vector();
2 26 Feb 07 jari 1028         }
2 26 Feb 07 jari 1029         
2 26 Feb 07 jari 1030         int counter = 0;
2 26 Feb 07 jari 1031         
2 26 Feb 07 jari 1032         for (int i = 0; i < gene.length; i++) {
2 26 Feb 07 jari 1033             if (gene[i] == min) {
2 26 Feb 07 jari 1034                 bins[0].add(new Integer(i));//otherwise, the min value(s) would never be binned, as the "for" loop below only checks to see
2 26 Feb 07 jari 1035                 //if the value is > currentLowerBound. bin[0] includes the currentLowerBound (i.e., min), all other bins include values > currentLowerBound
2 26 Feb 07 jari 1036                 counter++;
2 26 Feb 07 jari 1037             }
2 26 Feb 07 jari 1038         }
2 26 Feb 07 jari 1039         
2 26 Feb 07 jari 1040         
2 26 Feb 07 jari 1041         
2 26 Feb 07 jari 1042         float currentLowerBound = min;
2 26 Feb 07 jari 1043         float interval = (max - min) / 10;
2 26 Feb 07 jari 1044         float currentUpperBound = currentLowerBound + interval;
2 26 Feb 07 jari 1045         
2 26 Feb 07 jari 1046         for (int i = 0; i < bins.length; i++) {
2 26 Feb 07 jari 1047             currentUpperBound = currentLowerBound + interval;
2 26 Feb 07 jari 1048             for (int j = 0; j < gene.length; j++) {
2 26 Feb 07 jari 1049                 if ((gene[j] > currentLowerBound)&&(gene[j] <= currentUpperBound)) {
2 26 Feb 07 jari 1050                     bins[i].add(new Integer(j));
2 26 Feb 07 jari 1051                     counter++;
2 26 Feb 07 jari 1052                 }
2 26 Feb 07 jari 1053             }
2 26 Feb 07 jari 1054             
2 26 Feb 07 jari 1055             currentLowerBound = currentUpperBound;
2 26 Feb 07 jari 1056         }
2 26 Feb 07 jari 1057         
2 26 Feb 07 jari 1058         if (currentUpperBound < max) { // might happen because of rounding off, in which case the max value(s) won't get binned
2 26 Feb 07 jari 1059             for (int i = 0; i < gene.length; i++) {
2 26 Feb 07 jari 1060                 if (gene[i] == max) {
2 26 Feb 07 jari 1061                     bins[9].add(new Integer(i));
2 26 Feb 07 jari 1062                     counter++;
2 26 Feb 07 jari 1063                 }
2 26 Feb 07 jari 1064             }
2 26 Feb 07 jari 1065         }
2 26 Feb 07 jari 1066         
2 26 Feb 07 jari 1067         if (counter > gene.length) {
2 26 Feb 07 jari 1068             System.out.println("Warning: Mutual info: too many elements added to bins");
2 26 Feb 07 jari 1069         } else if (counter < gene.length) {
2 26 Feb 07 jari 1070             System.out.println("Warning: Mutual info: too few elements added to bins");
2 26 Feb 07 jari 1071         }      
2 26 Feb 07 jari 1072         return bins;
2 26 Feb 07 jari 1073     }
2 26 Feb 07 jari 1074     
2 26 Feb 07 jari 1075     private static double calculateSingleEntropy(float[] gene) {
2 26 Feb 07 jari 1076         Vector[] bins = makeBins(gene);
2 26 Feb 07 jari 1077         double[] probArray = new double[bins.length];
2 26 Feb 07 jari 1078         
2 26 Feb 07 jari 1079         for (int i = 0; i < probArray.length; i++) {
2 26 Feb 07 jari 1080             //System.out.println("bins[" + i + "].size() = " + bins[i].size());
2 26 Feb 07 jari 1081             //System.out.println("gene.length = " + gene.length);
2 26 Feb 07 jari 1082             int num = bins[i].size();
2 26 Feb 07 jari 1083             int denom = gene.length;
2 26 Feb 07 jari 1084             probArray[i] = (double)num/(double)denom;
2 26 Feb 07 jari 1085             //System.out.println("probArray[" + i + "] = " + probArray[i]);
2 26 Feb 07 jari 1086         }
2 26 Feb 07 jari 1087         
2 26 Feb 07 jari 1088         double entropy = 0;
2 26 Feb 07 jari 1089         
2 26 Feb 07 jari 1090         for (int i = 0; i < probArray.length; i++) {
2 26 Feb 07 jari 1091             if (probArray[i] != 0) {
2 26 Feb 07 jari 1092                 entropy = entropy + probArray[i]*(Math.log(probArray[i])/Math.log(2));
2 26 Feb 07 jari 1093             }
2 26 Feb 07 jari 1094         }
2 26 Feb 07 jari 1095         
2 26 Feb 07 jari 1096         entropy = (-1)*entropy;
2 26 Feb 07 jari 1097         return entropy;
2 26 Feb 07 jari 1098         
2 26 Feb 07 jari 1099     }
2 26 Feb 07 jari 1100     
2 26 Feb 07 jari 1101     private static double calculateJointEntropy(float[] gene1, float[] gene2) {
2 26 Feb 07 jari 1102         Vector[] bins1 = makeBins(gene1);
2 26 Feb 07 jari 1103         Vector[] bins2 = makeBins(gene2);
2 26 Feb 07 jari 1104         
2 26 Feb 07 jari 1105         float[][] jointProbMatrix = new float[bins1.length][];
2 26 Feb 07 jari 1106         
2 26 Feb 07 jari 1107         for (int i = 0; i < jointProbMatrix.length; i++) {
2 26 Feb 07 jari 1108             jointProbMatrix[i] = new float[bins1.length];
2 26 Feb 07 jari 1109         }
2 26 Feb 07 jari 1110         
2 26 Feb 07 jari 1111         for (int i = 0; i < jointProbMatrix.length; i++) {
2 26 Feb 07 jari 1112             for (int j = 0; j < jointProbMatrix[0].length; j++) {
2 26 Feb 07 jari 1113                 jointProbMatrix[i][j] = 0;
2 26 Feb 07 jari 1114             }
2 26 Feb 07 jari 1115         }
2 26 Feb 07 jari 1116         
2 26 Feb 07 jari 1117         int[] gene1States = new int[gene1.length];
2 26 Feb 07 jari 1118         int[] gene2States = new int[gene2.length];
2 26 Feb 07 jari 1119         
2 26 Feb 07 jari 1120         for (int i = 0; i < gene1States.length; i++) {//initializing
2 26 Feb 07 jari 1121             gene1States[i] = -1;
2 26 Feb 07 jari 1122             gene2States[i] = -1;
2 26 Feb 07 jari 1123         }
2 26 Feb 07 jari 1124         
2 26 Feb 07 jari 1125         for (int i = 0; i < bins1.length; i++) {
2 26 Feb 07 jari 1126             Vector currentBin = bins1[i];
2 26 Feb 07 jari 1127             if (currentBin.size() != 0) {
2 26 Feb 07 jari 1128                 for (int j = 0 ; j < currentBin.size(); j++) {
2 26 Feb 07 jari 1129                     int currentElement = ((Integer)(currentBin.get(j))).intValue();
2 26 Feb 07 jari 1130                     //System.out.println("currentElement from bins1 = " + currentElement);                   
2 26 Feb 07 jari 1131                     gene1States[currentElement] = i;
2 26 Feb 07 jari 1132                     //System.out.println("gene1States[" + currentElement + "] = " + gene1States[currentElement]);
2 26 Feb 07 jari 1133                 }
2 26 Feb 07 jari 1134             }
2 26 Feb 07 jari 1135         }
2 26 Feb 07 jari 1136
2 26 Feb 07 jari 1137         //System.out.println();
2 26 Feb 07 jari 1138  
2 26 Feb 07 jari 1139         for (int i = 0; i < bins2.length; i++) {
2 26 Feb 07 jari 1140             Vector currentBin = bins2[i];
2 26 Feb 07 jari 1141             if (currentBin.size() != 0) {
2 26 Feb 07 jari 1142                 for (int j = 0 ; j < currentBin.size(); j++) {
2 26 Feb 07 jari 1143                     int currentElement = ((Integer)(currentBin.get(j))).intValue();
2 26 Feb 07 jari 1144                     //System.out.println("currentElement from bins2 = " + currentElement);
2 26 Feb 07 jari 1145                     gene2States[currentElement] = i;
2 26 Feb 07 jari 1146                     //System.out.println("gene2States[" + currentElement + "] = " + gene2States[currentElement]);                    
2 26 Feb 07 jari 1147                 }
2 26 Feb 07 jari 1148             }
2 26 Feb 07 jari 1149         }  
2 26 Feb 07 jari 1150         
2 26 Feb 07 jari 1151         for (int i = 0; i < gene1States.length; i++) {
2 26 Feb 07 jari 1152             if (gene1States.length == -1) {
2 26 Feb 07 jari 1153                 System.out.println("Warning in mutual info: gene1States[" + i + "] has not been initialized");
2 26 Feb 07 jari 1154             }
2 26 Feb 07 jari 1155
2 26 Feb 07 jari 1156             if (gene2States.length == -1) {
2 26 Feb 07 jari 1157                 System.out.println("Warning in mutual info: gene2States[" + i + "] has not been initialized");
2 26 Feb 07 jari 1158             }            
2 26 Feb 07 jari 1159             
2 26 Feb 07 jari 1160         }
2 26 Feb 07 jari 1161
2 26 Feb 07 jari 1162         for (int i = 0 ; i < gene1States.length; i++) {
2 26 Feb 07 jari 1163             int state1 = gene1States[i];
2 26 Feb 07 jari 1164             int state2 = gene2States[i];
2 26 Feb 07 jari 1165             jointProbMatrix[state1][state2] = jointProbMatrix[state1][state2] + 1;
2 26 Feb 07 jari 1166         }
2 26 Feb 07 jari 1167         
2 26 Feb 07 jari 1168         for (int i = 0; i < jointProbMatrix.length; i++) {
2 26 Feb 07 jari 1169             for (int j = 0; j < jointProbMatrix[0].length; j++) {
2 26 Feb 07 jari 1170                 jointProbMatrix[i][j] = (jointProbMatrix[i][j])/(gene1States.length);
2 26 Feb 07 jari 1171             }
2 26 Feb 07 jari 1172            
2 26 Feb 07 jari 1173         }
2 26 Feb 07 jari 1174
2 26 Feb 07 jari 1175         
2 26 Feb 07 jari 1176         double jointEntropy = 0;
2 26 Feb 07 jari 1177         
2 26 Feb 07 jari 1178         
2 26 Feb 07 jari 1179         for (int i = 0; i < jointProbMatrix.length; i++) {
2 26 Feb 07 jari 1180             for (int j = 0; j < jointProbMatrix[0].length; j++) {
2 26 Feb 07 jari 1181                 if (jointProbMatrix[i][j] != 0) {
2 26 Feb 07 jari 1182                     jointEntropy = jointEntropy + (jointProbMatrix[i][j])*(Math.log(jointProbMatrix[i][j])/Math.log(2));
2 26 Feb 07 jari 1183                 }
2 26 Feb 07 jari 1184             }
2 26 Feb 07 jari 1185         }
2 26 Feb 07 jari 1186          
2 26 Feb 07 jari 1187         
2 26 Feb 07 jari 1188        jointEntropy = (-1)*jointEntropy;
2 26 Feb 07 jari 1189        
2 26 Feb 07 jari 1190        return jointEntropy;
2 26 Feb 07 jari 1191         
2 26 Feb 07 jari 1192     }
2 26 Feb 07 jari 1193   
2 26 Feb 07 jari 1194     */
2 26 Feb 07 jari 1195     
2 26 Feb 07 jari 1196 }

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