--- /dev/null
+#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