binaries/src/muscle/intmath.cpp

   1 #include "muscle.h"\r
   2 #include <math.h>\r
   3 \r
   4 PROB ScoreToProb(SCORE Score)\r
   5         {\r
   6         if (MINUS_INFINITY >= Score)\r
   7                 return 0.0;\r
   8         return (PROB) pow(2.0, (double) Score/INTSCALE);\r
   9         }\r
  10 \r
  11 static const double log2e = log2(exp(1.0));\r
  12 \r
  13 double lnTolog2(double ln)\r
  14         {\r
  15         return ln*log2e;\r
  16         }\r
  17 \r
  18 double log2(double x)\r
  19         {\r
  20         if (0 == x)\r
  21                 return MINUS_INFINITY;\r
  22 \r
  23         static const double dInvLn2 = 1.0/log(2.0);\r
  24 // Multiply by inverse of log(2) just in case multiplication\r
  25 // is faster than division.\r
  26         return log(x)*dInvLn2;\r
  27         }\r
  28 \r
  29 SCORE ProbToScore(PROB Prob)\r
  30         {\r
  31         if (0.0 == Prob)\r
  32                 return MINUS_INFINITY;\r
  33 //      return (SCORE) floor(INTSCALE*log2(Prob));\r
  34         return (SCORE) log2(Prob);\r
  35         }\r
  36 \r
  37 WEIGHT DoubleToWeight(double d)\r
  38         {\r
  39         assert(d >= 0);\r
  40         return (WEIGHT) (INTSCALE*d);\r
  41         }\r
  42 \r
  43 double WeightToDouble(WEIGHT w)\r
  44         {\r
  45         return (double) w / (double) INTSCALE;\r
  46         }\r
  47 \r
  48 SCORE DoubleToScore(double d)\r
  49         {\r
  50         return (SCORE)(d*(double) INTSCALE);\r
  51         }\r
  52 \r
  53 bool ScoreEq(SCORE s1, SCORE s2)\r
  54         {\r
  55         return BTEq(s1, s2);\r
  56         }\r
  57 \r
  58 static bool BTEq2(BASETYPE b1, BASETYPE b2)\r
  59         {\r
  60         double diff = fabs(b1 - b2);\r
  61         if (diff < 0.0001)\r
  62                 return true;\r
  63         double sum = fabs(b1) + fabs(b2);\r
  64         return diff/sum < 0.005;\r
  65         }\r
  66 \r
  67 bool BTEq(double b1, double b2)\r
  68         {\r
  69         return BTEq2((BASETYPE) b1, (BASETYPE) b2);\r
  70         }\r
  71 \r
  72 const double dLn2 = log(2.0);\r
  73 \r
  74 // pow2(x)=2^x\r
  75 double pow2(double x)\r
  76         {\r
  77         if (MINUS_INFINITY == x)\r
  78                 return 0;\r
  79         return exp(x*dLn2);\r
  80         }\r
  81 \r
  82 // lp2(x) = log2(1 + 2^-x), x >= 0\r
  83 double lp2(double x)\r
  84         {\r
  85         return log2(1 + pow2(-x));\r
  86         }\r
  87 \r
  88 // SumLog(x, y) = log2(2^x + 2^y)\r
  89 SCORE SumLog(SCORE x, SCORE y)\r
  90         {\r
  91         return (SCORE) log2(pow2(x) + pow2(y));\r
  92         }\r
  93 \r
  94 // SumLog(x, y, z) = log2(2^x + 2^y + 2^z)\r
  95 SCORE SumLog(SCORE x, SCORE y, SCORE z)\r
  96         {\r
  97         return (SCORE) log2(pow2(x) + pow2(y) + pow2(z));\r
  98         }\r
  99 \r
 100 // SumLog(w, x, y, z) = log2(2^w + 2^x + 2^y + 2^z)\r
 101 SCORE SumLog(SCORE w, SCORE x, SCORE y, SCORE z)\r
 102         {\r
 103         return (SCORE) log2(pow2(w) + pow2(x) + pow2(y) + pow2(z));\r
 104         }\r
 105 \r
 106 SCORE lp2Fast(SCORE x)\r
 107         {\r
 108         assert(x >= 0);\r
 109         const int iTableSize = 1000;\r
 110         const double dRange = 20.0;\r
 111         const double dScale = dRange/iTableSize;\r
 112         static SCORE dValue[iTableSize];\r
 113         static bool bInit = false;\r
 114         if (!bInit)\r
 115                 {\r
 116                 for (int i = 0; i < iTableSize; ++i)\r
 117                         dValue[i] = (SCORE) lp2(i*dScale);\r
 118                 bInit = true;\r
 119                 }\r
 120         if (x >= dRange)\r
 121                 return 0.0;\r
 122         int i = (int) (x/dScale);\r
 123         assert(i >= 0 && i < iTableSize);\r
 124         SCORE dResult = dValue[i];\r
 125         assert(BTEq(dResult, lp2(x)));\r
 126         return dResult;\r
 127         }\r
 128 \r
 129 // SumLog(x, y) = log2(2^x + 2^y)\r
 130 SCORE SumLogFast(SCORE x, SCORE y)\r
 131         {\r
 132         if (MINUS_INFINITY == x)\r
 133                 {\r
 134                 if (MINUS_INFINITY == y)\r
 135                         return MINUS_INFINITY;\r
 136                 return y;\r
 137                 }\r
 138         else if (MINUS_INFINITY == y)\r
 139                 return x;\r
 140 \r
 141         SCORE dResult;\r
 142         if (x > y)\r
 143                 dResult = x + lp2Fast(x-y);\r
 144         else\r
 145                 dResult = y + lp2Fast(y-x);\r
 146         assert(SumLog(x, y) == dResult);\r
 147         return dResult;\r
 148         }\r
 149 \r
 150 SCORE SumLogFast(SCORE x, SCORE y, SCORE z)\r
 151         {\r
 152         SCORE dResult = SumLogFast(x, SumLogFast(y, z));\r
 153         assert(SumLog(x, y, z) == dResult);\r
 154         return dResult;\r
 155         }\r
 156 \r
 157 SCORE SumLogFast(SCORE w, SCORE x, SCORE y, SCORE z)\r
 158         {\r
 159         SCORE dResult = SumLogFast(SumLogFast(w, x), SumLogFast(y, z));\r
 160         assert(SumLog(w, x, y, z) == dResult);\r
 161         return dResult;\r
 162         }\r
 163 \r
 164 double VecSum(const double v[], unsigned n)\r
 165         {\r
 166         double dSum = 0.0;\r
 167         for (unsigned i = 0; i < n; ++i)\r
 168                 dSum += v[i];\r
 169         return dSum;\r
 170         }\r
 171 \r
 172 void Normalize(PROB p[], unsigned n)\r
 173         {\r
 174         unsigned i;\r
 175         PROB dSum = 0.0;\r
 176         for (i = 0; i < n; ++i)\r
 177                 dSum += p[i];\r
 178         if (0.0 == dSum)\r
 179                 Quit("Normalize, sum=0");\r
 180         for (i = 0; i < n; ++i)\r
 181                 p[i] /= dSum;\r
 182         }\r
 183 \r
 184 void NormalizeUnlessZero(PROB p[], unsigned n)\r
 185         {\r
 186         unsigned i;\r
 187         PROB dSum = 0.0;\r
 188         for (i = 0; i < n; ++i)\r
 189                 dSum += p[i];\r
 190         if (0.0 == dSum)\r
 191                 return;\r
 192         for (i = 0; i < n; ++i)\r
 193                 p[i] /= dSum;\r
 194         }\r
 195 \r
 196 void Normalize(PROB p[], unsigned n, double dRequiredTotal)\r
 197         {\r
 198         unsigned i;\r
 199         double dSum = 0.0;\r
 200         for (i = 0; i < n; ++i)\r
 201                 dSum += p[i];\r
 202         if (0.0 == dSum)\r
 203                 Quit("Normalize, sum=0");\r
 204         double dFactor = dRequiredTotal / dSum;\r
 205         for (i = 0; i < n; ++i)\r
 206                 p[i] *= (PROB) dFactor;\r
 207         }\r
 208 \r
 209 bool VectorIsZero(const double dValues[], unsigned n)\r
 210         {\r
 211         for (unsigned i = 0; i < n; ++i)\r
 212                 if (dValues[i] != 0.0)\r
 213                         return false;\r
 214         return true;\r
 215         }\r
 216 \r
 217 void VectorSet(double dValues[], unsigned n, double d)\r
 218         {\r
 219         for (unsigned i = 0; i < n; ++i)\r
 220                 dValues[i] = d;\r
 221         }\r
 222 \r
 223 bool VectorIsZero(const float dValues[], unsigned n)\r
 224         {\r
 225         for (unsigned i = 0; i < n; ++i)\r
 226                 if (dValues[i] != 0.0)\r
 227                         return false;\r
 228         return true;\r
 229         }\r
 230 \r
 231 void VectorSet(float dValues[], unsigned n, float d)\r
 232         {\r
 233         for (unsigned i = 0; i < n; ++i)\r
 234                 dValues[i] = d;\r
 235         }\r
 236 \r
 237 double Correl(const double P[], const double Q[], unsigned uCount)\r
 238         {\r
 239         double dSumP = 0.0;\r
 240         double dSumQ = 0.0;\r
 241         for (unsigned n = 0; n < uCount; ++n)\r
 242                 {\r
 243                 dSumP += P[n];\r
 244                 dSumQ += Q[n];\r
 245                 }\r
 246         const double dMeanP = dSumP/uCount;\r
 247         const double dMeanQ = dSumQ/uCount;\r
 248 \r
 249         double dSum1 = 0.0;\r
 250         double dSum2 = 0.0;\r
 251         double dSum3 = 0.0;\r
 252         for (unsigned n = 0; n < uCount; ++n)\r
 253                 {\r
 254                 const double dDiffP = P[n] - dMeanP;\r
 255                 const double dDiffQ = Q[n] - dMeanQ;\r
 256                 dSum1 += dDiffP*dDiffQ;\r
 257                 dSum2 += dDiffP*dDiffP;\r
 258                 dSum3 += dDiffQ*dDiffQ;\r
 259                 }\r
 260         if (0 == dSum1)\r
 261                 return 0;\r
 262         const double dCorrel = dSum1 / sqrt(dSum2*dSum3);\r
 263         return dCorrel;\r
 264         }\r
 265 \r
 266 float Correl(const float P[], const float Q[], unsigned uCount)\r
 267         {\r
 268         float dSumP = 0.0;\r
 269         float dSumQ = 0.0;\r
 270         for (unsigned n = 0; n < uCount; ++n)\r
 271                 {\r
 272                 dSumP += P[n];\r
 273                 dSumQ += Q[n];\r
 274                 }\r
 275         const float dMeanP = dSumP/uCount;\r
 276         const float dMeanQ = dSumQ/uCount;\r
 277 \r
 278         float dSum1 = 0.0;\r
 279         float dSum2 = 0.0;\r
 280         float dSum3 = 0.0;\r
 281         for (unsigned n = 0; n < uCount; ++n)\r
 282                 {\r
 283                 const float dDiffP = P[n] - dMeanP;\r
 284                 const float dDiffQ = Q[n] - dMeanQ;\r
 285                 dSum1 += dDiffP*dDiffQ;\r
 286                 dSum2 += dDiffP*dDiffP;\r
 287                 dSum3 += dDiffQ*dDiffQ;\r
 288                 }\r
 289         if (0 == dSum1)\r
 290                 return 0;\r
 291         const float dCorrel = dSum1 / (float) sqrt(dSum2*dSum3);\r
 292         return dCorrel;\r
 293         }\r
 294 \r
 295 // Simple (but slow) function to compute Pearson ranks\r
 296 // that allows for ties. Correctness and simplicity\r
 297 // are priorities over speed here.\r
 298 void Rank(const float P[], float Ranks[], unsigned uCount)\r
 299         {\r
 300         for (unsigned n = 0; n < uCount; ++n)\r
 301                 {\r
 302                 unsigned uNumberGreater = 0;\r
 303                 unsigned uNumberEqual = 0;\r
 304                 unsigned uNumberLess = 0;\r
 305                 double dValue = P[n];\r
 306                 for (unsigned i = 0; i < uCount; ++i)\r
 307                         {\r
 308                         double v = P[i];\r
 309                         if (v == dValue)\r
 310                                 ++uNumberEqual;\r
 311                         else if (v < dValue)\r
 312                                 ++uNumberLess;\r
 313                         else\r
 314                                 ++uNumberGreater;\r
 315                         }\r
 316                 assert(uNumberEqual >= 1);\r
 317                 assert(uNumberEqual + uNumberLess + uNumberGreater == uCount);\r
 318                 Ranks[n] = (float) (1 + uNumberLess + (uNumberEqual - 1)/2.0);\r
 319                 }\r
 320         }\r
 321 \r
 322 void Rank(const double P[], double Ranks[], unsigned uCount)\r
 323         {\r
 324         for (unsigned n = 0; n < uCount; ++n)\r
 325                 {\r
 326                 unsigned uNumberGreater = 0;\r
 327                 unsigned uNumberEqual = 0;\r
 328                 unsigned uNumberLess = 0;\r
 329                 double dValue = P[n];\r
 330                 for (unsigned i = 0; i < uCount; ++i)\r
 331                         {\r
 332                         double v = P[i];\r
 333                         if (v == dValue)\r
 334                                 ++uNumberEqual;\r
 335                         else if (v < dValue)\r
 336                                 ++uNumberLess;\r
 337                         else\r
 338                                 ++uNumberGreater;\r
 339                         }\r
 340                 assert(uNumberEqual >= 1);\r
 341                 assert(uNumberEqual + uNumberLess + uNumberGreater == uCount);\r
 342                 Ranks[n] = (double) (1 + uNumberLess + (uNumberEqual - 1)/2.0);\r
 343                 }\r
 344         }\r
 345 \r
 346 FCOUNT SumCounts(const FCOUNT Counts[])\r
 347         {\r
 348         FCOUNT Sum = 0;\r
 349         for (int i = 0; i < 20; ++i)\r
 350                 Sum += Counts[i];\r
 351         return Sum;\r
 352         }\r