Mac binaries
[jabaws.git] / website / archive / binaries / mac / src / muscle / intmath.cpp
diff --git a/website/archive/binaries/mac/src/muscle/intmath.cpp b/website/archive/binaries/mac/src/muscle/intmath.cpp
new file mode 100644 (file)
index 0000000..40c25bb
--- /dev/null
@@ -0,0 +1,354 @@
+#include "muscle.h"\r
+#include <math.h>\r
+\r
+PROB ScoreToProb(SCORE Score)\r
+       {\r
+       if (MINUS_INFINITY >= Score)\r
+               return 0.0;\r
+       return (PROB) pow(2.0, (double) Score/INTSCALE);\r
+       }\r
+\r
+//#if  0\r
+//static const double log2e = log2(exp(1.0));\r
+//\r
+//double lnTolog2(double ln)\r
+//     {\r
+//     return ln*log2e;\r
+//     }\r
+//\r
+//double log2(double x)\r
+//     {\r
+//     if (0 == x)\r
+//             return MINUS_INFINITY;\r
+//\r
+//     static const double dInvLn2 = 1.0/log(2.0);\r
+//// Multiply by inverse of log(2) just in case multiplication\r
+//// is faster than division.\r
+//     return log(x)*dInvLn2;\r
+//     }\r
+//#endif\r
+\r
+//SCORE ProbToScore(PROB Prob)\r
+//     {\r
+//     if (0.0 == Prob)\r
+//             return MINUS_INFINITY;\r
+////   return (SCORE) floor(INTSCALE*log2(Prob));\r
+//     return (SCORE) log2(Prob);\r
+//     }\r
+\r
+WEIGHT DoubleToWeight(double d)\r
+       {\r
+       assert(d >= 0);\r
+       return (WEIGHT) (INTSCALE*d);\r
+       }\r
+\r
+double WeightToDouble(WEIGHT w)\r
+       {\r
+       return (double) w / (double) INTSCALE;\r
+       }\r
+\r
+SCORE DoubleToScore(double d)\r
+       {\r
+       return (SCORE)(d*(double) INTSCALE);\r
+       }\r
+\r
+bool ScoreEq(SCORE s1, SCORE s2)\r
+       {\r
+       return BTEq(s1, s2);\r
+       }\r
+\r
+static bool BTEq2(BASETYPE b1, BASETYPE b2)\r
+       {\r
+       double diff = fabs(b1 - b2);\r
+       if (diff < 0.0001)\r
+               return true;\r
+       double sum = fabs(b1) + fabs(b2);\r
+       return diff/sum < 0.005;\r
+       }\r
+\r
+bool BTEq(double b1, double b2)\r
+       {\r
+       return BTEq2((BASETYPE) b1, (BASETYPE) b2);\r
+       }\r
+\r
+//const double dLn2 = log(2.0);\r
+\r
+//// pow2(x)=2^x\r
+//double pow2(double x)\r
+//     {\r
+//     if (MINUS_INFINITY == x)\r
+//             return 0;\r
+//     return exp(x*dLn2);\r
+//     }\r
+\r
+//// lp2(x) = log2(1 + 2^-x), x >= 0\r
+//double lp2(double x)\r
+//     {\r
+//     return log2(1 + pow2(-x));\r
+//     }\r
+\r
+// SumLog(x, y) = log2(2^x + 2^y)\r
+//SCORE SumLog(SCORE x, SCORE y)\r
+//     {\r
+//     return (SCORE) log2(pow2(x) + pow2(y));\r
+//     }\r
+//\r
+//// SumLog(x, y, z) = log2(2^x + 2^y + 2^z)\r
+//SCORE SumLog(SCORE x, SCORE y, SCORE z)\r
+//     {\r
+//     return (SCORE) log2(pow2(x) + pow2(y) + pow2(z));\r
+//     }\r
+//\r
+//// SumLog(w, x, y, z) = log2(2^w + 2^x + 2^y + 2^z)\r
+//SCORE SumLog(SCORE w, SCORE x, SCORE y, SCORE z)\r
+//     {\r
+//     return (SCORE) log2(pow2(w) + pow2(x) + pow2(y) + pow2(z));\r
+//     }\r
+\r
+//SCORE lp2Fast(SCORE x)\r
+//     {\r
+//     assert(x >= 0);\r
+//     const int iTableSize = 1000;\r
+//     const double dRange = 20.0;\r
+//     const double dScale = dRange/iTableSize;\r
+//     static SCORE dValue[iTableSize];\r
+//     static bool bInit = false;\r
+//     if (!bInit)\r
+//             {\r
+//             for (int i = 0; i < iTableSize; ++i)\r
+//                     dValue[i] = (SCORE) lp2(i*dScale);\r
+//             bInit = true;\r
+//             }\r
+//     if (x >= dRange)\r
+//             return 0.0;\r
+//     int i = (int) (x/dScale);\r
+//     assert(i >= 0 && i < iTableSize);\r
+//     SCORE dResult = dValue[i];\r
+//     assert(BTEq(dResult, lp2(x)));\r
+//     return dResult;\r
+//     }\r
+//\r
+//// SumLog(x, y) = log2(2^x + 2^y)\r
+//SCORE SumLogFast(SCORE x, SCORE y)\r
+//     {\r
+//     if (MINUS_INFINITY == x)\r
+//             {\r
+//             if (MINUS_INFINITY == y)\r
+//                     return MINUS_INFINITY;\r
+//             return y;\r
+//             }\r
+//     else if (MINUS_INFINITY == y)\r
+//             return x;\r
+//\r
+//     SCORE dResult;\r
+//     if (x > y)\r
+//             dResult = x + lp2Fast(x-y);\r
+//     else\r
+//             dResult = y + lp2Fast(y-x);\r
+//     assert(SumLog(x, y) == dResult);\r
+//     return dResult;\r
+//     }\r
+//\r
+//SCORE SumLogFast(SCORE x, SCORE y, SCORE z)\r
+//     {\r
+//     SCORE dResult = SumLogFast(x, SumLogFast(y, z));\r
+//     assert(SumLog(x, y, z) == dResult);\r
+//     return dResult;\r
+//     }\r
+\r
+//SCORE SumLogFast(SCORE w, SCORE x, SCORE y, SCORE z)\r
+//     {\r
+//     SCORE dResult = SumLogFast(SumLogFast(w, x), SumLogFast(y, z));\r
+//     assert(SumLog(w, x, y, z) == dResult);\r
+//     return dResult;\r
+//     }\r
+\r
+double VecSum(const double v[], unsigned n)\r
+       {\r
+       double dSum = 0.0;\r
+       for (unsigned i = 0; i < n; ++i)\r
+               dSum += v[i];\r
+       return dSum;\r
+       }\r
+\r
+void Normalize(PROB p[], unsigned n)\r
+       {\r
+       unsigned i;\r
+       PROB dSum = 0.0;\r
+       for (i = 0; i < n; ++i)\r
+               dSum += p[i];\r
+       if (0.0 == dSum)\r
+               Quit("Normalize, sum=0");\r
+       for (i = 0; i < n; ++i)\r
+               p[i] /= dSum;\r
+       }\r
+\r
+void NormalizeUnlessZero(PROB p[], unsigned n)\r
+       {\r
+       unsigned i;\r
+       PROB dSum = 0.0;\r
+       for (i = 0; i < n; ++i)\r
+               dSum += p[i];\r
+       if (0.0 == dSum)\r
+               return;\r
+       for (i = 0; i < n; ++i)\r
+               p[i] /= dSum;\r
+       }\r
+\r
+void Normalize(PROB p[], unsigned n, double dRequiredTotal)\r
+       {\r
+       unsigned i;\r
+       double dSum = 0.0;\r
+       for (i = 0; i < n; ++i)\r
+               dSum += p[i];\r
+       if (0.0 == dSum)\r
+               Quit("Normalize, sum=0");\r
+       double dFactor = dRequiredTotal / dSum;\r
+       for (i = 0; i < n; ++i)\r
+               p[i] *= (PROB) dFactor;\r
+       }\r
+\r
+bool VectorIsZero(const double dValues[], unsigned n)\r
+       {\r
+       for (unsigned i = 0; i < n; ++i)\r
+               if (dValues[i] != 0.0)\r
+                       return false;\r
+       return true;\r
+       }\r
+\r
+void VectorSet(double dValues[], unsigned n, double d)\r
+       {\r
+       for (unsigned i = 0; i < n; ++i)\r
+               dValues[i] = d;\r
+       }\r
+\r
+bool VectorIsZero(const float dValues[], unsigned n)\r
+       {\r
+       for (unsigned i = 0; i < n; ++i)\r
+               if (dValues[i] != 0.0)\r
+                       return false;\r
+       return true;\r
+       }\r
+\r
+void VectorSet(float dValues[], unsigned n, float d)\r
+       {\r
+       for (unsigned i = 0; i < n; ++i)\r
+               dValues[i] = d;\r
+       }\r
+\r
+double Correl(const double P[], const double Q[], unsigned uCount)\r
+       {\r
+       double dSumP = 0.0;\r
+       double dSumQ = 0.0;\r
+       for (unsigned n = 0; n < uCount; ++n)\r
+               {\r
+               dSumP += P[n];\r
+               dSumQ += Q[n];\r
+               }\r
+       const double dMeanP = dSumP/uCount;\r
+       const double dMeanQ = dSumQ/uCount;\r
+\r
+       double dSum1 = 0.0;\r
+       double dSum2 = 0.0;\r
+       double dSum3 = 0.0;\r
+       for (unsigned n = 0; n < uCount; ++n)\r
+               {\r
+               const double dDiffP = P[n] - dMeanP;\r
+               const double dDiffQ = Q[n] - dMeanQ;\r
+               dSum1 += dDiffP*dDiffQ;\r
+               dSum2 += dDiffP*dDiffP;\r
+               dSum3 += dDiffQ*dDiffQ;\r
+               }\r
+       if (0 == dSum1)\r
+               return 0;\r
+       const double dCorrel = dSum1 / sqrt(dSum2*dSum3);\r
+       return dCorrel;\r
+       }\r
+\r
+float Correl(const float P[], const float Q[], unsigned uCount)\r
+       {\r
+       float dSumP = 0.0;\r
+       float dSumQ = 0.0;\r
+       for (unsigned n = 0; n < uCount; ++n)\r
+               {\r
+               dSumP += P[n];\r
+               dSumQ += Q[n];\r
+               }\r
+       const float dMeanP = dSumP/uCount;\r
+       const float dMeanQ = dSumQ/uCount;\r
+\r
+       float dSum1 = 0.0;\r
+       float dSum2 = 0.0;\r
+       float dSum3 = 0.0;\r
+       for (unsigned n = 0; n < uCount; ++n)\r
+               {\r
+               const float dDiffP = P[n] - dMeanP;\r
+               const float dDiffQ = Q[n] - dMeanQ;\r
+               dSum1 += dDiffP*dDiffQ;\r
+               dSum2 += dDiffP*dDiffP;\r
+               dSum3 += dDiffQ*dDiffQ;\r
+               }\r
+       if (0 == dSum1)\r
+               return 0;\r
+       const float dCorrel = dSum1 / (float) sqrt(dSum2*dSum3);\r
+       return dCorrel;\r
+       }\r
+\r
+// Simple (but slow) function to compute Pearson ranks\r
+// that allows for ties. Correctness and simplicity\r
+// are priorities over speed here.\r
+void Rank(const float P[], float Ranks[], unsigned uCount)\r
+       {\r
+       for (unsigned n = 0; n < uCount; ++n)\r
+               {\r
+               unsigned uNumberGreater = 0;\r
+               unsigned uNumberEqual = 0;\r
+               unsigned uNumberLess = 0;\r
+               double dValue = P[n];\r
+               for (unsigned i = 0; i < uCount; ++i)\r
+                       {\r
+                       double v = P[i];\r
+                       if (v == dValue)\r
+                               ++uNumberEqual;\r
+                       else if (v < dValue)\r
+                               ++uNumberLess;\r
+                       else\r
+                               ++uNumberGreater;\r
+                       }\r
+               assert(uNumberEqual >= 1);\r
+               assert(uNumberEqual + uNumberLess + uNumberGreater == uCount);\r
+               Ranks[n] = (float) (1 + uNumberLess + (uNumberEqual - 1)/2.0);\r
+               }\r
+       }\r
+\r
+void Rank(const double P[], double Ranks[], unsigned uCount)\r
+       {\r
+       for (unsigned n = 0; n < uCount; ++n)\r
+               {\r
+               unsigned uNumberGreater = 0;\r
+               unsigned uNumberEqual = 0;\r
+               unsigned uNumberLess = 0;\r
+               double dValue = P[n];\r
+               for (unsigned i = 0; i < uCount; ++i)\r
+                       {\r
+                       double v = P[i];\r
+                       if (v == dValue)\r
+                               ++uNumberEqual;\r
+                       else if (v < dValue)\r
+                               ++uNumberLess;\r
+                       else\r
+                               ++uNumberGreater;\r
+                       }\r
+               assert(uNumberEqual >= 1);\r
+               assert(uNumberEqual + uNumberLess + uNumberGreater == uCount);\r
+               Ranks[n] = (double) (1 + uNumberLess + (uNumberEqual - 1)/2.0);\r
+               }\r
+       }\r
+\r
+FCOUNT SumCounts(const FCOUNT Counts[])\r
+       {\r
+       FCOUNT Sum = 0;\r
+       for (int i = 0; i < 20; ++i)\r
+               Sum += Counts[i];\r
+       return Sum;\r
+       }\r