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